text:
# Dynamical Rotational Instability at Low 𝑇/𝑊

## 1 Introduction

Rotational instabilities are potentially important in the evolution of massive stellar cores. For example, the core of a massive star that has been prevented from collapsing to neutron star densities by centrifugal forces may rotate rapidly enough for rotational instabilities to develop (Hayashi, Eriguchi, & Hashimoto, 1998, 1999). Massive cores spun up by accretion from a binary companion (Wagoner, 1984) and the remnants of compact binary mergers (Rasio & Shapiro, 1992; Zhuge, Centrella, & McMillan, 1994) may also reach fast enough rotation rates to become unstable. Rotational instabilities may produce detectable gravitational radiation (Schutz, 1989; Thorne, 1996). If enough angular momentum is shed, full collapse to a neutron star or black hole may occur.

We focus here on global rotational instabilities that arise in fluids from growing azimuthal modes $`e^{im\varphi }`$ (Tassoul, 1978). Dynamical instabilities are driven by hydrodynamics and gravity, and grow on the order of the dynamical timescales of the system. In objects with Maclaurin-like rotation laws, dynamical rotational instabilities are generally expected to arise at fairly high values of $`\beta T/|W|0.27`$ in Newtonian gravity (see New, Centrella, & Tohline and references therein), and $`\beta 0.250.26`$ in general relativity (Shibata, Baumgarte, & Shapiro, 2000). Here, $`T`$ is the rotational kinetic energy, and $`W`$ is the gravitational potential energy. Dynamical instability may set in at lower values of $`\beta `$ for thick, self-gravitating disks (Pickett, Durisen, & Davis, 1996; Woodward, Tohline, & Hachisu, 1994) and in tori (Tohline & Hachisu, 1990). Secular instabilities develop on longer timescales, and arise from dissipative processes such as viscosity and gravitational radiation reaction. They are expected to set in for relatively low values of $`\beta 0.14`$, although the rotation law, polytropic index, and general relativity can affect the value of the instability limit (Imamura, et al., 1995; Stergioulas & Friedman, 1998).

Studies of scenarios for the formation of centrifugally-hung cores have generally concluded that these objects will have $`\beta <0.27`$ (Tohline, 1984; Eriguchi & Müller, 1985; Müller & Eriguchi, 1985). Direct numerical simulations of axisymmetric stellar collapse by Zwerger & Müller (1997) also indicate that, for the cases studied, it is difficult for collapsing cores to reach and exceed $`\beta 0.27`$. Overall, it is generally assumed that centrifugally-hung cores will evolve secularly, as in the “fizzler” scenario (Hayashi, Eriguchi, & Hashimoto, 1998, 1999).

In this Letter, we report on new 3-D numerical simulations of rotating stellar cores that exhibit dynamical rotational instability at relatively low values of $`\beta 0.14`$. These models are computed using Newtonian gravity with no backreaction; this is a reasonable approximation for a stellar core of mass $`M1.4M_{}`$ hung up at a radius $`R100`$km, so that $`(GM/Rc^2)0.02`$. Our results show that compact objects can become dynamically unstable in a region of parameter space previously thought to be solely the province of secular instabilities.

## 2 Numerical Simulations

Rotating stellar cores are initially assumed to be axisymmetric equilibrium polytropes with equation of state $`P=K\rho ^\mathrm{\Gamma }=K\rho ^{1+1/N}`$, with $`\mathrm{\Gamma }=1.3`$ ($`N=3.33`$). The angular velocity distribution is given by the so-called $`j`$constant rotation law, $`\mathrm{\Omega }^2=j_o^2/(d^2+\varpi ^2)^2`$, where $`\varpi `$ is the cylindrical radius and $`d`$ is an arbitrary constant (Hachisu, 1986). As $`d0`$, the specific angular momentum approaches the constant value $`j_o`$. Here, we use $`d=0.2`$. The models are computed on a uniformly spaced $`(\varpi ,z)`$ grid using the self-consistent field method of Hachisu (1986). Varying the axis ratio yields equilibria with different values of $`\beta `$; see New, Centrella, & Tohline (2000) for details. Figure Dynamical Rotational Instability at Low $`T/W`$ shows density contours in the $`xz`$ plane for the 4 models that we evolved here, $`\beta =0.090,0.12,0.14,`$ and $`0.18`$. Notice that, with the exception of the $`\beta =0.090`$ model, the density maxima are toroidal; similar structures were used as initial data in the collapse models of Rampp, Ruffert, & Müller (1998).

These initial models were introduced into 2 different 3-D hydrocodes. The first of these is the L code discussed in New, Centrella, & Tohline (2000), with a cylindrical $`(\varpi ,z,\varphi )`$ grid of $`64\times 64\times 128`$ zones. The second is the PPM code discussed by Brown (2000), with a Cartesian $`(x,y,z)`$ grid of $`128\times 128\times 128`$ zones. Both codes have the same resolution along the $`x`$ and $`z`$axes, with the equatorial radius $`R_\mathrm{E}`$ of the initial model extending out to zone 48 in $`x`$-direction. The cylindrical code imposes equatorial plane symmetry, whereas the Cartesian code imposes no symmetries. Random perturbations of 1% were imposed on the densities, and the models were evolved forward in time.

## 3 Results

The model with initial $`\beta =0.14`$ exhibits a dynamical instability that is very similar in simulations performed with both the cylindrical and Cartesian hydrocodes. These runs were stopped when the loss of mass from the edge of the grid became significant. Figure Dynamical Rotational Instability at Low $`T/W`$ shows density contours in the equatorial plane from the run with the cylindrical code. The green (lightest greyscale) contour in the first frame delineates the toroidal density maximum seen in Fig. Dynamical Rotational Instability at Low $`T/W`$. As time proceeds, this torus “pinches off” to leave a single high density region, indicative of an $`m=1`$ mode. At late times, the appearance of higher density contours signals that the dense region is starting to collapse. Time is measured in units of $`t_\mathrm{D}=(R_\mathrm{E}^3/GM)^{1/2}`$, which is the dynamical time for a sphere of mass $`M`$ and radius $`R_\mathrm{E}`$. For $`R_\mathrm{E}200`$ km, $`t_\mathrm{D}6`$ ms. 3-D animations of the evolution of this model, using the same contour levels as in Figure Dynamical Rotational Instability at Low $`T/W`$ and looking down the $`z`$axis towards the equatorial plane, are available for both the cylindrical and Cartesian runs. Additional animations, in which the model is viewed from the side, are also available for the cylindrical and Cartesian runs. \[Note to editors: Here, the words “cylindrical” and “Cartesian” should link to the animations in each case.\]

We can quantify the instability by examining various Fourier components in the density distribution. We examine the density in a ring of fixed $`\varpi `$ and $`z`$ using a complex azimuthal Fourier decomposition, where the amplitudes $`C_m`$ of the various components $`m`$ are defined by $`C_m(\varpi ,z)=(1/2\pi )_0^{2\pi }\rho (\varpi ,z,\phi )e^{im\phi }𝑑\phi `$ (Tohline, Durisen, & McCollough, 1985). The normalized amplitude is then $`A_m=C_m/C_0`$, where $`C_0(\varpi ,z)=\overline{\rho }(\varpi ,z)`$ is the mean density in the ring, and the phase angle of the $`m^{\mathrm{th}}`$ component is $`\varphi _m(\varpi ,z)=\mathrm{tan}^1\left[\mathrm{Im}(A_m)/\mathrm{Re}(A_m)\right]`$. We write the phase angle $`\varphi _m=\sigma _mt`$, where $`\sigma _m=d\varphi /dt`$ is the eigenfrequency and $`W_m=\sigma _m/m`$ is the pattern speed of the $`m^{\mathrm{th}}`$ mode.

Figures Dynamical Rotational Instability at Low $`T/W`$ and Dynamical Rotational Instability at Low $`T/W`$ show the growth of the amplitudes $`|A_m|`$ for the first four Fourier components in the runs with $`T/|W|=0.14`$ on the cylindrical and Cartesian hydrocodes, respectively. For the cylindrical code, these amplitudes were calculated in the equatorial plane in a ring of width $`\mathrm{\Delta }\varpi =1/48`$ at radius $`\varpi =0.32`$; see New, Centrella, & Tohline (2000) for details. For the Cartesian code, the amplitudes were computed on a circle of radius $`\varpi =0.32`$ using a nonuniform discretization, which avoids grid boundaries, and a linear interpolation of density.

Similar results were found in rings at other values of $`\varpi `$. The evolution in both runs is dominated by an exponentially growing $`m=1`$ mode (thick solid line). Once this mode is well into its exponential growth phase, an $`m=2`$ mode also develops, followed at later times by $`m=3`$ and $`m=4`$ modes. (The constant amplitude $`m=4`$ signal in the Cartesian run is due to the geometry of the numerical grid.) Since these modes grow rapidly on dynamical timescales, they signal dynamical instability. Numerical values for the growth rates, eigenfrequencies, and pattern speeds are shown in Table 1 for $`m=1`$ and $`m=2`$. Notice that the cylindrical and Cartesian runs yield similar results and that the pattern speeds $`W_1W_2`$, indicating that the $`m=2`$ mode is a harmonic of the $`m=1`$ mode. Longer runs with larger grids and higher resolution are needed to obtain reliable values for the $`m=3`$ and $`m=4`$ modes.

The model with $`\beta =0.18`$ also exhibits dynamical instability. For the evolutions run with both codes, the $`m=1`$ mode grew at very close to the same rate, starting at about the same time; see Table 1. In the run on the cylindrical code, the modes $`m=1,2,3,`$ and $`4`$ appear sequentially in time and $`W_1W_2`$, as in the model with $`\beta =0.14`$. In the run on the Cartesian code, the $`m=1`$ and $`m=2`$ modes both grew at about the same rate, starting at about the same time; however, $`W_1W_2`$. At this stage, we do not understand this difference in the behavior of the $`m=2`$ mode in this case. We plan to investigate it using higher resolution simulations. The $`\beta =0.12`$ model was run for $`40t_\mathrm{D}`$, at which time the $`m=1`$ mode was just beginning to grow; longer runs with higher resolution are needed to confirm the development of instability in this case. The model with $`\beta =0.09`$ was run for $`35t_\mathrm{D}^1`$, and showed no growing modes. In the event that we can confirm instability for $`\beta =0.12`$, we will return to the $`\beta =0.09`$ case and run it further.

Previous studies of the $`m=2`$ bar-mode instability using the cylindrical hydrocode showed that significant motion of the system center of mass could develop at late times, resulting in a spurious $`m=1`$ signal (New, Centrella, & Tohline, 2000). For the simulations reported here, we monitored the position of the overall center of mass and verified that it underwent no systematic motion during the development of the $`m=1`$ mode, in both the cylindrical and Cartesian codes. When the runs with $`T/|W|=0.14`$ and $`T/|W|=0.18`$ were stopped, some small, spurious center of mass motion ($`<`$ one zone) had developed in both codes. This effect is possibly related to the loss of mass from the grid as the outer regions expand; runs with larger grids are needed to determine the late time behavior of these models. Previous studies of the bar instability using the Cartesian code (Brown, 2000) showed that angular momentum losses due to numerical inaccuracies can be signficant. For the $`\beta =0.14`$ and $`\beta =0.18`$ simulations with the Cartesian code, the artificial angular momentum losses amounted to less than $`0.4\%`$ and $`1.5\%`$, respectively.

## 4 Discussion

Our results demonstrate that dynamical instability can occur in differentially rotating polytropes with low values of $`\beta 0.14`$, and that this instability has a strong, $`m=1`$ character. Figure Dynamical Rotational Instability at Low $`T/W`$ shows that the unstable axisymmetric equilibria each have a torus of dense material centered on the rotation axis within the model. This suggests that this instability may be related to the ones found by Woodward, Tohline, & Hachisu (1994) and Tohline & Hachisu (1990) in their studies of polytropic tori with $`N=3/2`$. We note that Pickett, Durisen, & Davis (1996) also found a dominant $`m=1`$ instability that set in at a relatively low $`\beta 0.20`$, in a centrally condensed model (without an off-center density maximum). Their model was an $`N=3/2`$ polytrope with an $`n^{}=2`$ rotation law. The $`n^{}=2`$ rotation law produces configurations with a strong concentration of angular momentum in the outer regions of the model. We plan to carry out more detailed studies to investigate the character of the unstable modes seen in our simulations and their properties for various values of $`N`$ and the $`j`$-constant rotation law parameter $`d`$.

If such instability occurs in a centrifugally-hung core, collapse to neutron star densities may result. Our simulations do show that the density is increasing at the end of the unstable runs. Further studies are needed to see how dense the remnant actually becomes, and whether the $`m=1`$ mode will result in the dense region moving at a velocity comparable to those of actual neutron stars (Popov, S. et al., 2000).

Dynamical instability will also produce gravitational radiation. Although longer runs are needed to obtain the full gravitational waveforms, we can estimate the properties of the gravitational wave emission from the initial stages of the instability. For a massive stellar core with $`M1.4M_{}`$ and $`R_\mathrm{E}200`$km, the peak amplitude will occur at a frequency of roughly $`f200`$Hz. This peak amplitude will be $`h10^{24}r_{20}^1`$ for $`\beta =0.14`$, and $`h10^{23}r_{20}^1`$ for $`\beta =0.18`$. Here, $`r_{20}`$ is the distance to the source in units of 20 Mpc. Emission from such unstable cores may be detectable with advanced ground-based interferometers, such as LIGO-II. An even more optimistic scenario for detectable gravitational radiation occurs if these instabilities arise in supermassive stars111 We note that the appropriate equation of state for supermassive stars is the $`n=3`$ polytropic equation of state. (New & Shapiro, 2000). For example, if the instability occurs when a supermassive star contracts to the point that $`(GM/Rc^2)1/15`$, an approximate value for uniformly rotating stars (Baumgarte & Shapiro 1999), we estimate the frequency of the gravitational radiation to be $`f3.5\times 10^3`$Hz, with amplitude $`h10^{18}r_{20}^1`$ for $`\beta =0.14`$ and $`h10^{17}r_{20}^1`$ for $`\beta =0.18`$. Such signals would be easily detectable by the space-based LISA detector.

We thank K. Thorne and J. Tohline for stimulating discussions, and E. Mamikonyan for assistance with the animations. This work was supported in part by NSF grants PHY-9722109, PHY-0070892, NSF cooperative agreement ACI-9619020 through computing resources provided by the National Partnership for Advanced Computational Infrastructure (NPACI) at the San Diego Supercomputer Center (SDSC), and the North Carolina Supercomputing Center (NCSC). A portion of this work was performed under the auspices of the U.S. Department of Energy by Los Alamos National Laboratory under contract W-7405-ENG-36. The cylindrical simulations were run on the T3E at SDSC, and the Cartesian models on the T90 at NCSC.


format: text

text:
# References

ITEP-TH-49/00

hep-th/0010068

RG flows on the phase spaces and the $`\tau `$ functions for the generic Hamiltonian systems

A. S. Gorsky

Institute of Theoretical and Experimental Physics, Moscow 117259, B. Cheryomushkinskaya 25

## Abstract

We discuss the generic definition of the $`\tau `$ function for the arbitrary Hamiltonian system. The different approaches concerning the deformations of the curves and surfaces are compared. It is shown that the Baker-Akhiezer function for the secondary integrable system of the Toda lattice type can be identified with the coherent wave function of the initial dynamical system. The $`\tau `$ function appears to be related to the filling of the interior of the classical trajectory by coherent states. Transition from dispersionless to dispersionful Toda lattice corresponds to the quantization of the initial dynamical system.

1. Recently it was recognized that it is possible to attribute the so called $`\tau `$ function to any analytic curve in two dimensions . This can be done as follows. In the simplest case one considers the evolution of the curve when the area inside grows linearly while the moments of the curve being fixed. In more general situation all moments become new variables (”times”) too. It appears that in all cases the relevant object which encodes the information on the properties of the deformation of the curve is the quasiclassical $`\tau `$ function of the Toda lattice with the additional constraint which selects the solution to the hierarchy.

On the other hand it was shown in  that it is possible to attribute the prepotential which is the quasiclassical version of the $`\tau `$ function to any holomorphic dynamical system. It plays the role of the generating function for the analogue of the S-duality transformation in the context of Hamiltonian dynamics. Moreover it was argued that generically there exists a pair of dynamical systems such that their coordinates and actions get interchanged. Actually one can say that along this viewpoint the quasiclassical $`\tau `$ function can be attributed to the Riemann surface which is the solution to the equation of motion of the holomorphic system.

The third general setting involving prepotential is the theory of the effective actions. It appears that generically effective actions manifest a lot of universality property which is governed by the hidden integrability. The examples are the identification of the effective action of the N=2 SYM theory with $`\tau `$ function of the quasiclassical hierarchy  and similar integrable structures behind d=2 theories which have been uncovered in . Along this way times in the hierarchy are identified with the coupling constants in the corresponding theories. In more general terms integrability amounts from the deformations of the topological theories or topological sectors in nontopological theories by the proper observables .

The goal of this note is to attempt to glue these ideas altogether and formulate the place of $`\tau `$ function for the arbitrary Hamiltonian system without any appealing to the integrability.

2. Let us remind the key points from . One can consider complex coordinates $`\overline{z},z`$ and the curve determined by the equation

$$\overline{z}=S(z)$$
(1)
The Schwarz function $`S(z)`$ is analytic in a domain including the curve. One more ingredient to be defined is the map of the exterior of the curve to the exterior of the unit disk

$$\omega (z)=\frac{z}{r}+\underset{j}{}p_jz^j$$
(2)
where $`\omega `$ is defined on the unit circle. It is useful to introduce the moments of the curve

$$t_n=\frac{1}{2\pi in}z^nS(z)𝑑z$$
(3)
$$v_n=\frac{1}{2\pi i}z^nS(z)𝑑z$$
(4)
$$v_0=log|z|𝑑z$$
(5)
which yield the following expansion for the Schwarz function

$$S(z)=kt_kz^{k1}+t_0z^1+kv_kz^{k1}$$
(6)
Let us define the generating function

$$S(z)=_z\mathrm{\Omega }(z)$$
(7)
with the following expansion

$$\mathrm{\Omega }(z)=\underset{k=1}{}t_kz^k+t_0logz\underset{k=1}{}v_kz^kk^11/2v_0$$
(8)
One can easily prove the following relations

$$_{t_0}\mathrm{\Omega }(z)=\mathrm{log}\omega (z)$$
(9)
$$_{t_n}\mathrm{\Omega }(z)=(z^n(\omega ))_++1/2(z^n(\omega ))_0$$
(10)
$$_{\overline{t_n}}\mathrm{\Omega }(z)=(S^n(\omega ))_++1/2(S^n(\omega ))_0$$
(11)
The symbol $`(S(\omega ))_+`$ means the truncated Laurent series with only positive powers of $`\omega `$ kept and the $`(S(\omega ))_0`$ is the constant term in the series. The differential $`d\mathrm{\Omega }`$

$$d\mathrm{\Omega }=Sdz+log\omega dt_0+(H_kdt_k\overline{H_k}d\overline{t_k})$$
(12)
provides the Hamiltonians and $`\mathrm{\Omega }`$ itself can be immediately identified with the generating function for the canonical transformation from the pair $`(z,\overline{z})`$ to the canonical pair $`(t_0,log\omega )`$.

The dynamical equations read

$$_{t_n}S(z)=_zH_n(z)$$
(13)
$$_{\overline{t_n}}S(z)=_z\overline{H_n}(z)$$
(14)
and the consistency of (13) (14) yields the zero-curvature condition which amounts to the equations of the dispersionless Toda lattice hierarchy. The first equation of the hierarchy reads as follows

$$_{t_1\overline{t_1}}^2\varphi =_{t_0}e^{_{t_0}\varphi }$$
(15)
where $`_{t_0}\varphi =2logr`$. The Lax operator L coincides with $`z(\omega )`$ and its eigenfunction - Baker-Akhiezer function looks as follows $`\mathrm{\Psi }=e^{\frac{\mathrm{\Omega }}{h}}`$. Hamiltonians are expressed in terms of the Lax operator

$$H_k=(L^k)_++1/2(L^k)_0$$
(16)
3. In this section we shall interpret the evolution above in terms of dynamics on the phase space assuming that $`\overline{z},z`$ pair yields the phase space of some dynamical system. The curve itself corresponds to the energy level of this dynamical system. With this setup it is clear that Poisson bracket between $`\overline{z}`$ and $`z`$ is fixed by the standard symplectic form . The next step is the identification of $`logw`$ as angle variable hence the identification of the area inside the curve $`t_0`$ as the action variable is evident.

The holomorphic variables on the phase space are usually considered as the operators in the Fock space so we identify $`z=b;\overline{z}=b^+`$. Typically the creation- annihilation operators are expressed it terms of the $`(p,q)`$ variables however for our needs we assume the action-angle representation. The variable $`z=b`$ can be expressed in terms of action-angle variables in the particular dynamical system and just this expression provides the Toda lattice Lax operator

$$z(logw,t_0)=L(logw,t_0)$$
(17)
Baker-Akhiezer (BA) function is solution to the equation

$$L\mathrm{\Psi }(z,t_0)=z\mathrm{\Psi }$$
(18)
and the complex $`z`$ plane is identified with the surface where the spectral parameter of the Toda system lives on.

Now we turn to clarification of the meaning of the BA function in the generic dynamical system. The answer follows from the equation above; BA function is nothing but the coherent wave function in the action representation. Indeed the coherent wave function is the eigenfunction of the creation operator

$$\widehat{b}\mathrm{\Psi }=b\mathrm{\Psi }$$
(19)
The main properties of the coherent states can be summarized as follows . The coherent state can be represented in the form

$$|\alpha >=D(\alpha )|0>$$
(20)
where $`|0>`$ is the normalized vacuum vector annihilated by the operator b and the $`D(\alpha )=exp(\alpha b^+\overline{\alpha }b)`$ yields the action of the Heisenberg-Weyl group. The operators $`D(\alpha )`$ obey the following relation

$$D(\alpha )D(\beta )=e^{2iIm(\alpha \overline{\beta })}D(\beta )D(\alpha )$$
(21)
The coherent states $`|\alpha >`$ and $`|\beta >`$ are not orthogonal moreover the system of the coherent states $`[|\alpha >]`$ is overcompleted generically. To get the complete system one can proceed as follows. The lattice on the complex plane $`\alpha _{nm}=n\omega _1+m\omega _2`$ can be introduced with the area of the elementary cell equals to the elementary Planck cell $`\pi `$. In this case the system of the coherent states $`|\alpha _{nm}>`$ can be expressed in terms of the standard theta-functions on the torus and turns out to be complete.

Note that the system of the coherent states can be defined for the generic group as follows

$$|\psi _g>=T(g)|\psi _0>$$
(22)
If H is the stabilizer of the vector $`|\psi _0>`$ the coherent state is defined by the point x=x(g) from the homogeneous space G/H corresponding to the element g. Moreover the state corresponding to $`|x>`$ can be identified with the one dimensional projector $`P_x=|x><x|`$.

Now we are ready to recognize the meaning of the generating function $`\mathrm{\Omega }`$. From the equations above it becomes clear that it is the generating function for the canonical transformations from the $`b,b^+`$ representation to the angle-action variables. Let us also remark that the role of the Orlov-Shulman operator becomes transparent. This operator looks like $`b^+b`$ and its eigenvalue counts the ”number of particles”.

Having identified the BA function for the generic system let us consider the role of the $`\tau `$ function in the generic case. To this aim it is convenient to use the following formulae for the $`\tau `$ function

$$\tau (t,W)=<t,\overline{t}|W>$$
(23)
where the bra vector depends on times while the ket vector is fixed by the so called point of Grassmanian. One more useful representation is provided by the fermionic language

$$\tau (t,W)=\frac{<N|\mathrm{\Psi }(z_1)\mathrm{}\mathrm{\Psi }(z_N)|W>}{\mathrm{\Delta }(z)}$$
(24)
where $`\mathrm{\Delta }(z)`$ is Vandermonde determinant.

Let us briefly comment on the definition of the point of the Grassmanian W. Generally speaking the Grassmanian itself can be considered as a collection of all fermionic Bogolyubov transforms of the vacuum state. Hence we can say that W belongs to the Grassmanian if it is annihilated by some linear combination of the fermionic creation and annihilation operators. Equivalently one could consider the following definition

$$|W>=S|0>,S=exp\underset{nm}{}A_{nm}\overline{\psi }_{n1/2}\psi _{m1/2}$$
(25)
that is S can be considered as the element of $`GL(\mathrm{},C)`$.

The consideration above suggests the following picture behind the definition of the $`\tau `$ function. The peculiar classical trajectory of the dynamical system yields the curve on the phase space. Then the domain inside the trajectory is filled by the coherent states for this particular system. Since the coherent state occupies the minimal cell of the phase space the number of the coherent states packed inside the domain is finite and equals N. Since there is only one coherent state per cell for the complete set it actually behaves like a fermion implying a kind of the fermionic representation.

Therefore we can develop the second dynamical system of the Toda lattice type based on the generic dynamical system. The number of the independent time variables in the Toda system amounts from the independent parameters in the potential in the initial system plus additional time attributed to the action variable. The situation is similar to the consideration of the pair of the dual dynamical systems in . The subtle point in that paper is the choice of the Hamiltonians for the dual systems. It seems that consideration in this note provides the unified viewpoint on this issue. Indeed formulae above suggest the set of commuting Hamiltonians for the dual system. Let us emphasize that the choice of the particular initial dynamical system amounts to the choice of the particular solution to the Toda lattice hierarchy.

It is important that the energy level of the dynamical system doesn’t develop the discontinuity in some range of the energies. This happens only if the system approaches separatrices generically existing on the phase space.

Consider the example of the oscillator. Corresponding phase trajectories

$$E=p^2+q^2$$
(26)
are circles on the phase space (ellipses in the generic case) and the natural complex variables are

$$z=p+iq$$
(27)
If we calculate the expression for $`v_0`$ the following prepotential for the oscillator emerges

$$=1/2I^2logI3/4I^2$$
(28)
where I is the action variable. To get the first dispersionless Toda equation one has to consider the phase trajectories with the shifted center. The position of the center is parameterized by times $`t_1,\overline{t_1}`$. If the Toda equation with three nonvanishing times is solved the prepotential above can be reproduced. Let us also note that the dual Hamiltonian in the sense of  is

$$H=b=\sqrt{t_0}e^{i_{t_0}}$$
(29)
The wave function of the dual system coincides with the coherent wave function in the action representation.

The prepotential for the complex system is defined by its spectral curve

$$a_D=\frac{}{a}$$
(30)
where $`a_D`$ and $`a`$ are the integrals of the meromorphic differential over A and B cycles. The variable $`a`$ is just the action variable. The S duality transformation $`aa_D`$ actually maps the region of small energies to the region of the large energies. Prepotential obeys the Matone type  relation which has the meaning of the Ward identity

$$\frac{}{log\mathrm{\Lambda }}=constH$$
(31)
where H is the Hamiltonian of the complex system and $`\mathrm{\Lambda }`$ is the scale factor. The simple calculation shows that such relation doesn’t‘t hold in the real case.

Let us remind how the prepotential provides the S-duality transformation for the complex dynamical system . The action variables in dynamical system are the integrals of meromorphic differential $`\lambda `$ over the $`A`$-cycles on the spectral curve. The reason for the $`B`$-cycles to be discarded is simply the fact that the $`B`$-periods of $`\lambda `$ are not independent of the $`A`$-periods. On the other hand, one can choose as the independent periods the integrals of $`\lambda `$ over any lagrangian subspace in $`H_1(𝐓_b;Z)`$.

This leads to the following structure of the action variables in the holomorphic setting. Locally over a disc in $`B`$ one chooses a basis in $`H_1`$ of the fiber together with the set of $`A`$-cycles. This choice may differ over another disc. Over the intersection of these discs one has a $`Sp(2m,Z)`$ transformation relating the bases which has a natural interpretation as a S-duality transformation. Altogether they form an $`Sp(2m,Z)`$ bundle. It can be easily shown the two form:

$$dI^idI_i^D$$
(32)
vanishes. Therefore one can always locally find a function $``$ , such that:

$$I_i^D=\frac{}{I^i}$$
(33)
The angle variables are uniquely reconstructed once the action variables are known. Since the real system to some extend can be considered as a real section of the dynamical system one could expect the analogue of the S-duality for the real system too.

An important outcome from our identification is the natural recipe for the transition from the dispersionless to dispersionful Toda lattices. Indeed at the first step of quantization the system one imposes the Bohr-Sommerfeld quantization condition on the action variable. This is in agreement with the discretization of $`t_0`$ variable in Toda lattice. The Planck constant in the initial system plays the transparent role in the Toda lattice.

4. Let us argue now that there is clear correspondence between the integrable dynamics behind the deformation of the curves which are the phase space trajectories for the real dynamical systems and the deformation of the surfaces which represent the trajectories of the complex dynamical systems. Remind that to formulate the integrable dynamics behind the deformation of the curves one introduces the following main objects; generating function $`\mathrm{\Omega }`$, $`\tau `$ function depending on time variables $`t_k`$ which are the moments of the exterior of the classical trajectory as well as the dual variables $`v_k`$ representing the moments of interior.

It is not a difficult task to recognize their counterparts in the case of the surfaces. First let us represent the time variables in the following form

$$v_k=z^k𝑑\mathrm{\Omega }(z)$$
(34)
$$t_k=z^k𝑑\mathrm{\Omega }(z)$$
(35)
and compare it with the definition of times in the case of surfaces

$$T_n=res_{\xi =0}\xi ^ndS_{SW}$$
(36)
$$\frac{}{T^n}=\frac{1}{2\pi in}res_0\xi ^ndS_{SW}$$
(37)
where $`S_{SW}`$ in the context of the N=2 SYM theories is the so-called Seiberg-Witten meromorphic differential whose derivatives with respect of the moduli of the curve are the holomorphic differentials. We see that the variables in cases of curves and surfaces are differed just by the Legendre transform. The variable $`t_0`$ is nothing but the variable

$$a=_A𝑑S$$
(38)
where the integration is over the A cycle on the spectral surface is the complex dynamical system with only one degree of freedom is considered. The definition of conjugated variables $`v_n`$ is more delicate issue since naively there are no the direct analogue of the integrals over the B cycles

$$a_D=_B𝑑S$$
(39)
in the real case. However it can be seen that the analogous variable for the one degree of freedom can be identified with variable $`v_0`$ from . If there are several degrees of freedom the corresponding variables $`a_i`$ are defined as the integrals over the complete set of A cycles.

It is interesting to recognize the real analogue of the modulus of the spectral curve in the complex case

$$\tau =\frac{da_D}{da}$$
(40)
which plays the role of the effective coupling constant in the field theory context. In the problem under consideration the corresponding variable is

$$logr=\frac{dv_0}{dt_0}$$
(41)
and also exhibits the behaviour similar to the perturbative regime in the field theory. In what follows we shall see more indications that dynamics of the curve has a pertubative regime as a field theory counterpart.

In the context of the Seiberg-Witten solution to the N=2 SUSY YM theory the Whitham dynamics has a transparent renormalization group meaning . Indeed the derivative of the prepotential with respect to the first time in SUSY YM context looks as follows

$$\frac{}{T_1}=a\frac{}{a}2=2H(a)$$
(42)
where the term at the r.h.s. coincides with the Hamiltonian for the dynamical system providing the corresponding spectral curve. The dependence $`H(a)`$ on the action variable is the intrinsic characteristics of the dynamical system, moreover this term has an anomalous nature from the field theory viewpoint. The analogous equation has been also found for the real case  where the observed quadratic dependence on the action in the anomaly term corresponds to the perturbative limit in the field theory framework.

Given such RG interpretation in the complex case one could look for the reasonable RG structure behind the real dynamical system. The following starting point can be chosen. Assume that the RG scale is fixed by the energy corresponding to the classical trajectory or by the variable $`t_0=\mu `$. Then the phase space is naturally divided into the ”IR theory” inside the trajectory and the ”UV theory” outside. Therefore a kind of the RG problem can be formulated. To this aim introduce the classical expectation value of the observable in the theory

$$<O>=𝑑p𝑑qO(p,q)\rho (p,q)$$
(43)
where $`\rho (p,q)`$ is the phase space density (Wigner function). The vacuum expectation value at the scale $`\mu `$ is naturally introduced if the integration over the phase space up to the energy $`\mu `$ is performed.

Now we can formulate the following observation - correlators in UV theory in the outer region amount to the coupling constants in the inner region and vice versa

$$\tau (t_0,t_k)=<expt_kO_k>$$
(44)
where the matrix element is taken in the IR theory. The key point is that usually effective action $`S(\mu )`$ identified with the prepotential is defined via integrating out all modes up to the scale $`\mu `$. From  it becomes clear that is has to be done in a clever way; one has to integrate out the part of the phase space outside of the phase trajectory corresponding to the energy E. Then the Riemann-Hilbert problem is defined on the RG scale $`\mu `$ \- one has to divide the whole theory into UV and IR parts consistently with the RG flows. Remark that interpretation of the area inside the curve as the RG parameter agrees with the interpretation of the size of the matrix in the matrix model as RG scale . Let us emphasize that the picture emerged is formally identical to the consideration in  where $`\tau `$ function approach amounts to the explicit calculation of the RG dependence of the correlators. Along this viewpoint one could expect that the similar correlators which now carry the information about the spectrum of the system could be found. Therefore instead of the determination of the RG invariant characteristics of the field theories we expect the RG like evaluation of the spectral invariants which would distinguish different behaviour of the Hamiltonian systems.

Let us emphasize that the Toda lattice equations are written for the variable $`b(\omega )`$ which yields the dependence of the creation operator on the angle variable. If only zero time is involved then the evolution doesn‘t change the system however with all times switched on the dynamics proceeds in the space of the dynamical systems or can be interpreted as the RG flows. The integrals of motion in the Toda system are the ”classical correlators” in the IR theory.

5. It was shown in  that the $`\tau `$ function of the curve allows two additional complementary interpretations. First, it has the interpretation of the energy of the fermions with the unit charge inside the curve with the Coulomb interaction taken into account. Secondly it can be considered as the partition function of the normal matrix model with the potential yielding the set of times in the Toda lattice.

Let us speculate on the appearance of the matrix model in a given context. To this aim consider the phase space as a two dimensional noncommutative plane. Evidently the Planck constant $`h`$ fixes the scale of noncommutativity due to the canonical commutation relation

$$[x_1,x_2]=ih$$
(45)
The Moyal multiplication law for the functions on the phase space reads

$$f(x)g(x)=exp(ih_x_y)_{x=y}g(y)f(x)$$
(46)
Now let us remind the notion of the Morita equivalence. Qualitatively it means that for the rational noncommutativity ($`h=p/q`$) the noncommutative theory is equivalent to the commutative theory on the rescaled noncommutative manifold. The functions on the commutative manifold become matrix valued with the size of the matrix fixed by the noncommutativity.

Let us apply these circle of ideas to the generic dynamical system. Suppose that the Planck constant is a rational number $`h=\frac{1}{N}`$. If N is large we are dealing with the quasiclassical approximation. Applying the Morita equivalence, all functions become $`N\times N`$ matrix valued and the system becomes effectively classical. Therefore Morita transform maps the quantum system into the classical matrix theory. The ”classical” partition function of the matrix model is the prepotential for the dynamical system.

One more indication of the effective classicality of the matrix model follows from the matrix relation $`[M^+,M]=0`$ which corresponds to the ”classical” commutation relation. On the other hand the corresponding Planck constant in the Morita equivalent system is $`\stackrel{~}{h}=N`$ which indicates the regime of the strong noncommutativity on the phase plane. Along this viewpoint one could expect the appearance of the peculiar objects relevant for the such regime, namely the noncommutative solitons  which in the simplest case correspond to the projector operators $`|n><n|`$. One the other hand we have already remarked that the coherent states also have the description in terms of the projector operators on the homogeneous spaces. This implies the intriguing conjecture to consider the matrix model as the system of the interacting noncommutative solitons.

Let us make one more remark concerning the different interpretations of the prepotential. It has a interpretation of the energy in the matrix model, action for the c=1 string and somewhat involved interpretation in the context of particle dynamics. This could presumable have the following picture behind; it is energy for the membrane and therefore the action for string. With this picture in mind the matrix model corresponds to the M(atrix) model representation for the membrane. This is in a rough agreement with the chain of T duality transformations of the phase space. Let us also note that in the context of c=1 string ”times” correspond to the vacuum expectation values of tachyons with the different momenta which could be related to the known phenomena of an appearance of the noncommutative solitons(coherent states) as the result of the tachyon condensation .

6. To conclude we have argued that the Toda lattice governs both the deformations of the complex surfaces and analytic curves. These deformations can be considered as the motion in the space of the complex and real dynamical systems respectively or as a kind of RG flows in a peculiar system. The dynamics with respect to the time $`t_0`$ can be considered as the RG evolution with $`t_0`$ playing the role of the scale factor. One could also expect that the discrete symmetries could provide an important tool for investigation of the RG flows like in .

The Baker-Akhiezer function was definitely identified as the coherent wave function in the action-angle representation and the generating functional has the interpretation of the generating function of the canonical transformation from holomorphic to action-angle representation. The role of the quantization of the initial dynamical system in the integrability framework is clarified.

The author is grateful to V. Roubtsov, P. Wiegmann and A. Zabrodin for the interesting discussions and University of Angers where the paper has been completed for the hospitality. The work was partially supported by grants INTAS-99-1705 and CRDF-RP1- 2108 and by grant for senior scientist fellowship of CNRS 2000.


text:
# Simulations of Merging Clusters of Galaxies

## 1. Introduction

Irregular cluster morphologies in many X-ray images of cluster as well as indications from optical observations imply that many clusters are not relaxed. Hence merging is a common phenomenon in clusters of galaxies. Such mergers of subclusters are very energetic events, which affect clusters strongly, e.g. shocks emerge which are important because they are the major heating source for the intra-cluster gas and moreover they can accelerate particles.

While observations provide only snapshots of the different evolutionary stages of mergers, simulations are the only way to follow the evolution. The models reveal where shocks emerge, where they move and how strong they are. With the help of simulated X-ray images at different stages of a merger a connection between X-ray morphology and the dynamical state can be established and the effects on mass determination can be investigated. Magnetohydrodynamic simulations show how the magnetic field is enhanced during mergers and how the overall structure of the magnetic field is changed due to a merger.

## 2. Simulation Methods

In order to perform realistic simulations three-dimensional calculations are required and the different cluster components must be taken into account. It is necessary to follow the evolution of the dark matter as well as the intra-cluster gas. Dark matter (and galaxies) can be regarded as collisionless particles and can therefore be modeled by N-body simulations. The gas is best simulated by hydrodynamic calculations. For these hydrodynamic calculations two different methods have mainly been used: (1) Smoothed Particle Hydrodynamics (Lagrangian approach, e.g. Evrard 1990; Dolag et al. 1999; Takizawa 1999; Takizawa & Naito 2000) in which the gas is treated as particles and (2) grid codes (Eulerian approach, e.g. Schindler & Müller 1993; Bryan et al. 1994; Roettiger et al. 1996, 1997, 1998, 1999a,b; Ricker 1998; Quilis et al. 1998) in which the simulation volume is divided into cells. Fortunately, the choice of simulation method seems to be irrelevant. Calculations with both methods yield very similar results.

## 3. Shocks

Mergers produce shocks in the intra-cluster gas. These shocks are of particular interest for particle acceleration models. The strongest shocks emerge after the collision of subclusters, when these shock propagate outwards along the original collision axis (Schindler & Müller 1993; Roettiger et al. 1999a). However, even these shocks are relatively mild shocks with a maximum Mach number of about 3. When a dense subcluster falls into a cluster a shock emerges already before the core passage: a bow shock is visible in front of the infalling subcluster (Roettiger et al. 1997). In general, the shock structure is found to be more filamentary at early epochs and quasi-spherical at low redshifts (Quilis et al. 1998).

Observationally, the shocks are best visible in temperature maps, because they show up as steep temperature gradients. For such maps spatially resolved X-ray spectroscopy is necessary which can be performed now with high accuracy with the new X-ray observatories XMM and Chandra.

## 4. Other Effects of Mergers

Apart from shocks, mergers have many other interesting effects. For example, the X-ray luminosity increases during the collision of two subclusters (Schindler & Müller 1993). The reason is that the gas is compressed, i.e. the gas density is increased, and as the X-ray emission is proportional to the square of the density we see enhanced X-ray emission during the passage. During the core passage and at each rebounce an increase in the magnetic field is visible (Dolag et al. 1999, 2000 (in prep.); Roettiger et al. 1999b; see also Section 6). Mergers also cause a lot of turbulence. Off-centre collisions produce in addition angular momentum (Ricker 1998; Roettiger et al. 1998). The shocks heat primarily the ions as has been shown in simulations which treat ions and electrons separately (Chièze et al. 1998; Takizawa 1999).

Observationally mergers cannot only be identified by multiple X-ray maxima, but also by isophote twisting with centroid shift and elongations: the dark matter component is always elongated along the collision axis, while the gas is first elongated along the collision axis but during the core passage it is pushed out perpendicular to the collision axis, so that later an elongation perpendicular to the collision axis can be seen (Schindler & Müller 1993). Also offsets between the collisionless component and the gas have been found (Roettiger et al. 1997).

## 5. Mass Determination in Merging Clusters

Mass determination in clusters with the X-ray method can be affected strongly during mergers (Evrard et al. 1996; Roettiger et al. 1996; Schindler 1996). The reason is that during the merging process there can be quite strong deviations from two assumptions necessary for the mass determination – hydrostatic equilibrium and spherical symmetry. For example, at the positions of shocks the gas is not in hydrostatic equilibrium. Shocks cause gradients – both in the temperature and in the density – and can cause therefore an overestimation of the mass. Locally, this can lead to a mass estimate up to two times the true mass. Substructure on the other hand tends to flatten the azimuthally averaged profile and hence leads to an underestimation of the mass, in extreme cases to deviations of 50% of the true mass.

In some cases, these deviations can be corrected for, e.g. in clusters in which substructures are well distinguishable, the disturbed part can be excluded from the mass analysis and a good mass estimate can be obtained. But in general, mass determinations in non-relaxed clusters should be done very cautiously.

## 6. Simulations of Mergers with Magnetic Fields

Radio halos require the existence of magnetic fields in clusters on scales of a few Mpc. Therefore magnetohydrodynamic calculations have been performed (Dolag et al. 1999; Roettiger et al. 1999b). It has been found that the initial field distribution is irrelevant for the final structure of the magnetic field. The structure is dominated only by the cluster collapse. Faraday rotation measurements can be reproduced by the simulations for magnetic fields of the order of $`\mu `$G. Most important for the amplification of the magnetic field are shear flows, while the compression of the gas is of minor importance. Mergers change the local magnetic field strength as well as the structure of the cluster-wide field. At early stages of the merger, filamentary structures prevail, which break down later and leave a stochastically ordered magnetic field.

## References

Bryan, G.L. et al. 1994, ApJ, 437, L5
Chièze, J.-P., Alimi J.-M., Teyssier, R. 1998, ApJ, 495, 630
Dolag, K., Bartelmann, M., Lesch, H. 1999, A&A, 348, 351
Evrard, A.E. 1990, A&A, 363, 349
Evrard, A.E., Metzler, C.A., Navarro, J.N. 1996, ApJ, 469, 494
Quilis, V., Ibáñez, J.M., Sáez, D. 1998, ApJ, 502, 518
Ricker, P.M. 1998, ApJ, 496, 670
Roettiger, K., Burns, J.O., Loken, C. 1996, ApJ, 473, 651
Roettiger, K., Burns, J.O., Stone J.M. 1999a, ApJ, 518, 603
Roettiger, K., Loken, C., Burns, J.O. 1997, ApJS, 109, 307
Roettiger, K., Stone J.M., Burns, J.O. 1999b, ApJ 518, 594
Roettiger, K., Stone J.M., Mushotzky, R.F., 1998, ApJ, 493, 62
Schindler, S. 1996, A&A, 305, 756
Schindler, S., Müller, E. 1993, A&A, 272, 137
Takizawa, M. 1999, ApJ, 529, 514
Takizawa, M., Naito, T., 2000, ApJ, 535, 586


text:
# Surface molecular dynamics simulation with two orthogonal surface steps: how to beat the particle conservation problem

## I introduction

Classical Molecular Dynamics (MD) simulation of crystal surfaces is generally inadequate for the description of physical situations involving large density and structure changes including some surface reconstructions, incomplete melting, thermal preroughening and roughening, etc. due to the constraint of particle conservation. A single monatomic step can act as a source or sink of atoms: given a request for lateral density variation the step can move normally to its direction, to accommodate to it. Pairs of steps of opposite sign are generally easy to implement in a slab simulation: but this is not optimal, for it leads to a very limited size terrace, with non-negligible step-step interactions ($`1/L^2`$, if $`L`$ is the terrace size), and general hindering of the step movement. To minimize this problem, one can introduce a single step per unit cell, which will weakly interact only with its own replicas ($`L`$ being in this case the entire lateral size of the simulation cell) through periodic boundary conditions (PBC). To this end, the PBC must be modified so that a layer $`n`$ on the cell right hand side is joined at the left hand side not with itself, but with layer $`n+1`$. In this manner the single step is introduced topologically. That scheme was applied with success to the simulation of the lateral density increase of the reconstructed hexagonal overlayer of Au$`(001)`$ . The movement of the step was indeed observed, with an enhanced mobility connected in that case with the incommensurability between first and second layer .

However, a lateral density adjustment will in general take place in both directions on the surface and not just along one. The impossibility for the density to re-adjust in both directions (as implied by the unidirectional single step construction, and by rigid overall PBC parallel to the step) can also hinder interplanar sliding which might in real life require simultaneous sliding in both planar directions. The natural solution to this problem is to allow the same freedom in both orthogonal surface directions. A practical realization of this task is presented here through the implementation of two orthogonal steps. We will first of all examine the practical implementation and the applicability to surfaces of unbent as well as of bent crystals, in the light of the recently introduced variable curvature molecular dynamics method . We will then present direct examples of application, exploiting as done previously  gold surfaces as a test case.

## II The double step

A surface with regularly repeated steps separated by large terraces is a vicinal surface. A method for generating a vicinal surface for simulation is to rotate a slab by an angle $`\theta =\mathrm{arctan}(L/h)`$ with a flat surface and dimensions $`L_x\times L_y\times L_z`$ around the step axis (called $`x`$) (with $`L`$ the distance between the steps and $`h`$ the interplanar separation) and to “cleave” the crystal surface accordingly, generating a sequence of monatomic steps. If $`L`$ is chosen equal to $`L_y`$, a single step is produced. In a typical high symmetry direction of an fcc crystal, such as $`(111)`$ or $`(001)`$, the stacking of the layers is $`ABCABC\mathrm{}`$ or $`ABABAB\mathrm{}`$. If a single step is produced, the portion of terrace above the step and the one below will belong to layers of different sublattices, thus causing a mismatch in the periodic boundary conditions (PBC). This problem is solved through the introduction of an additional translation vector in the PBC, absorbing the mismatch between the two different layers . An equivalent solution is to modify the dimensions and orientation of the simulation cell, which will no longer be built with an integer number of unit cells nor with walls parallel to high-symmetry planes of the crystal . We prefer here the former approach.

If we now wish to add a single step along the orthogonal $`y`$ direction as well, a problem in the PBC will obviously appear also in the $`x`$ direction, and we will need two rotations of the slab. The first, through an angle $`\theta L_y/h`$ about the $`x`$ axis brings the $`y`$ axis into $`y^{}`$. The second, through an angle $`\varphi L_x/h`$, about $`y^{}`$.

The dimension of the cell is chosen so that a particle and its image lie at the same height $`z`$.

## III Application to a flat surface

As an example, we consider a Au$`(001)`$ surface, whose hexagonal reconstructed layer is incommensurate, with a cell approximately $`(n\times 5)`$ with $`n30`$. In the $`[110]`$ (henceforth $`y`$) direction, 6 surface rows, shifted alternatingly so as to form a quasi triangular layer, are stacked on 5 rows of the underlying square second layer. In the $`[1\overline{1}0]`$ (henceforth $`x`$) direction, one soliton-like extra atomic row is present every $`n`$ second layer rows. The surface finally appears as in Fig. 1 , where a reconstructed Au$`(001)`$ surface, 20 layer thick, is shown. Although there is in reality a single large terrace, that appears in the form of four sub-terraces of height 1,0,0, -1. The corrugation given by the $`(5\times 1)`$ reconstruction periodicity is clearly visible.

To demonstrate how the two orthogonal steps are going to act as grand-canonical sources of surface atoms, we will exploit one well established but rather peculiar property of this surface: Au$`(001)`$ increases its lateral density upon heating . Physically this is due to the competing role of the $`d`$\- and $`s`$-orbitals of gold: the consequence of this competition is a giant outwards relaxation, resulting in a consequent planar contraction of the first layer.

Since the (6 over 5) commensurability pinning in the $`y`$ direction is much stronger than the ($`n+1`$ over $`n`$) pinning in the $`x`$ direction, the lateral contractions upon heating will generally be different between x and y. Our goal here will be to discover what our approach predicts for them using the two-steps geometry. different between x and y. Our goal here will be to discover by simulation what should happen in reality by using the two-step geometry. We modeled the surface with a slab with 20 layers of 1800 atoms each. The two surfaces were reconstructed with 31 rows over 30 in the $`x`$ direction and 6 rows over 5 in the $`y`$ direction. Our interatomic potential was the glue model of Ercolessi et al.  with a many body term mimicking the metallic cohesion of valence s and d electrons. Newton’s equations were numerically integrated with standard algorithms, and velocity rescaling is used to control temperature. After an initial quench to relax the slab atomic positions, we heated slowly (rate 100 Kelvin/nsec) up to $`1000K`$. Fig. 2 shows a top view of the surface of Fig. 1 before and after it was heated in this way.

At the beginning ($`T=0`$) the four sub-terraces had equal dimensions (left panel). At $`T=900K`$ (right panel) the situation has now changed, with a visible withdrawal (even if superposed with meandering) of the x-step along the positive y direction, but no visible shift of the y-step along the x direction. As a result, the two lower subterraces 0, -1 have expanded at the expense of the two upper ones 1, 0. The relative x-step withdrawal is about 4%. Since that happens in one direction only, this measures exactly the relative lateral density increase in the top layer of our Au(001) on heating model between zero and 900 K. This value equals very closely the magnitude of the known surface lateral density increase of real Au(001) from room temperature to 900 K . Not surprisingly it is mostly the weakly pinned x-step that acted as a source of atoms for the terrace density increase, while hardly anything happened to the strongly pinned y-step. The reconstruction periodicity drift with temperature from $`(5\times 30)`$ to roughly $`(5\times 28.5)`$, which is precisely the trend seen experimentally .

## IV Double step and variable curvature

A second interesting application of the double step geometry is in conjunction with the Variable Curvature Molecular Dynamics (VCMD) method . That method allows the simulation of a bent slab, the radius of bending curvature treated as an additional degree of freedom. In the force-free case, the slab will adopt a spontaneous curvature in response to the surface stress difference of its two faces . Alternatively, a bending force or constraint can be applied to the slab, to force a given curvature. The associated fields of strain and stress in the slab will change its overall free energy and in particular that of its surfaces . In response, the equilibrium surface lateral density may also generally change, and we can gauge that change with our double step method. Special care must be taken in this case concerning PBC: if the slab is bent, a particle with velocity parallel to the bending direction which crosses the boundary of the simulation cell, must remain on the same ‘orbit’, thus the translation vector added to the usual PBC refolding vectors should be also parallel to the bending direction. With this in mind, we have considered two important cases under bending: Au(001) $`(1\times 5)`$ and the Au(111) with its $`(\sqrt{3}\times n)`$ reconstruction. Details of both calculations will be reported in a forthcoming paper . Here we just give some preliminary results concerning Au$`(111)`$. Under compressive strain, the density mismatch between the substrate and the overlayer decreases, and the $`\sqrt{3}\times n`$ periodicity should change, in particular n should increase.

We started with a single $`(\sqrt{3}\times 11)`$ reconstructed domain ($`n=11`$ is the equilibrium value in the glue model ), with 6 solitons which are parallel to the $`[11\overline{2}]`$ (henceforth $`x`$) direction and two orthogonal steps, a total of about 1200 surface atoms, for a slab thickness of 40 layers. We applied a finite curvature in the $`[1\overline{1}0]`$ (henceforth $`y`$) direction, and focused on the face under compression (see Fig. 3).

Fig. 4 shows the phase difference of the overlayer over the second layer with and without bending strain. The flat surface shows 6 solitons along the $`x`$ direction (the peak of one of the steps is visible), whereas under a compressive strain of about 4% the period of reconstruction has doubled to $`(\sqrt{3}\times 22)`$ and there are only 3 solitons. The step parallel to the soliton advances in this case, accepting atoms from the terrace. On the other hand, the other step (orthogonal to the solitons) remains pinned, most likely due to strong commensurability. Hence in comparison with the Au(001) studied in the previous chapter, we have step motion parallel rather than orthogonal to the soliton, and we also demonstrate step advancing rather than withdrawing.

## V Conclusions and perspectives

Two orthogonal single steps can be built topologically on a simulated crystal surface through periodic boundary conditions. This device is shown to permit, under favorable conditions, an effectively grand canonical study of the surface with molecular dynamics (MD) simulation, otherwise strictly canonical. Any change of surface coverage will reflect in an advancement or retraction of the two steps, that can be naturally quite anisotropic. The method is shown to work by means of two case studies. the thermal increase of lateral surface density of Au(001), and the strain-induced decrease of that of Au(111).

In both cases only one step eventually moved, but the presence of two orthogonal steps allowed to separate the lateral expansion in two directions. The conclusion that the contraction is dominant in one direction with respect to the other, althoug physically plausible, could not have been drawn with a single step geometry.

In the case of surfaces under strain, we are presently applying the same method to the curved Au $`(100)`$, and the results will appear in a future paper .

ACKNOWLEDGMENTS

We acknowledge support from MURST COFIN99, and from INFM/F. Research of D. P. in Stuttgart is supported by the Alexander Von Humboldt foundation.

## Figure captions

* Figure 1: The 20-layers slab of the reconstructed Au (100) surface with two orthogonal steps; the corrugations given by the $`5\times 1`$ periodicity is clearly visible.
* Figure 2: Top view of the Au (100) surface before and after heating at $`900K`$ (brighter atoms are higher, thus the top left terrace is the highest one). Before heating, the four subterraces have the same size (left); after heating, the horizontal step has retracted, signalling a lateral density increase (right). The numbers (0,-1,1) indicate the relative height of the four subterraces.
* Figure 3: First few layers of the simulated Au (111) slab under compressive strain along the $`1\overline{1}0`$ ($`y`$) direction. Both the two orthogonal steps and the externally imposed curvature can be observed($`T=800K`$).
* Figure 4: Change in periodicity (phase difference of the overlayer over the second layer) of the Au (111) surface reconstruction under compressive strain: the flat surface shows 6 solitons along the $`y`$ direction (lower curve), whereas under a compressive strain of about $`4\%`$ the period of reconstruction has doubled (upper curve).


text:
# Surface-state electron dynamics in noble metals

## 1 Introduction

Shockley surface states are known to exist near the Fermi level in the $`\mathrm{\Gamma }`$-L projected bulk band gap of the (111) surfaces of the noble metals Cu, Ag, an Au. The wave functions of these crystal-induced surface states are localized near the surface and decay into the solid, in contrast to the Bloch waves propagating into the bulk. Hence, these states form a quasi two-dimensional (2D) electron gas, which overlaps in energy and space with the three-dimensional (3D) substrate, and represent a playground for lifetime investigations.

Various techniques have been used to measure the lifetime broadening of surface-state electrons and holes near the Fermi level. High-resolution angle-resolved photoemission spectroscopy (ARPES) has been used to investigate the lifetime of surface-state holes. In contrast to this technique, which is restricted to lifetime measurements of occupied states, it has been demonstrated that scanning tunneling spectroscopy (STS) offers the possibility to measure the lifetime of long-living surface states above and below the Fermi level. The lifetime of excited holes at the band edge of the partially occupied Shockley surface states on Ag(111), Cu(111), and Au(111) was accurately determined by Li et al and by Kliewer et al with the use of STS. Recently, STS has been used to measure lifetimes of hot surface-state and surface-resonance electrons, as a function of excess energy.

On the (111) surfaces of Cu, Ag, and Au, the Shockley surface state at the center of the surface Brillouin zone ($`𝐤_{}=0`$) lies just below the Fermi level, with $`EE_F=445`$, $`67`$, and $`505\mathrm{meV}`$, respectively , $`E`$ and $`E_F`$ representing the surface-state energy and the Fermi level. The dispersion relation of these states is displayed in Fig. 1, which shows that they are free-electron-like with effective masses of $`0.42`$, $`0.44`$, and $`0.28m_e`$ ($`m_e`$ is the free-electron mass) that account for the potential variation parallel to the surface. This figure clearly shows that the decay of Shockley surface-state electrons and holes in these materials may proceed either through the coupling with bulk states (3D channel) or through the coupling, within the surface-state band itself, with surface states of different wave vector parallel to the surface (2D channel).

While 3D free-electron calculations of the decay rate of excited holes at the surface-state band edge ($`𝐤_{}=0`$) of the noble metals are known to predict significantly longer lifetimes than those observed with either ARPES or STS, recent measurements have shown that, in the case of surface-state electrons above the Fermi level, the experimental values are comparable to the calculated lifetimes of bulk free-electrons with the same energy. That the decay of surface-state holes is dominated by 2D electron-electron interactions, but screened by the underlying 3D electron system, has been demonstrated recently, showing an excellent agreement with experiment. In this paper, we extend these theoretical investigations to the case of hot surface-state and surface-resonance electrons in Cu(111), and show that major contributions to the electron-electron interaction of surface-state electrons above the Fermi level come from the underlying bulk electrons. We also present self-consistent calculations of the screened interaction, and investigate the role that both the partially occupied 2D surface-state band and the underlying 3D electron gas play in the relaxation mechanism.

The contribution to measured lifetimes of surface-state electrons in the noble metals arising from electron-phonon interactions has been discussed in . Here we focus on the investigation of the energy-dependent inelastic lifetimes that are due to electron-electron interactions.

## 2 Theory

We take an arbitrary Fermi system of interacting electrons and consider an excited electron or hole interacting with the Fermi sea. The electron-hole decay rate, i.e., the probability per unit time for the probe particle to scatter from an initial state $`\psi _i(𝐫)`$ of energy $`E_i`$ to some available final state $`\psi _f(𝐫)`$ of energy $`E_f`$ ($`|E_fE_F|<|E_iE_F|`$), by carrying the Fermi gas from the many-particle ground state to some excited many-particle state, may be obtained by using the golden rule of time-dependent perturbation theory. By keeping terms of first-order in the screened interaction, one finds (we use atomic units throughout, i.e., $`e^2=\mathrm{}=m_e=1`$):

$`\tau ^1=`$ $`2{\displaystyle \underset{f}{}}{\displaystyle d𝐫d𝐫^{}\psi _i^{}(𝐫)\psi _f^{}(𝐫^{})}`$ (3)
$`\times \mathrm{Im}\left[W(𝐫,𝐫^{};|E_iE_f|)\right]\psi _i(𝐫^{})\psi _f(𝐫),`$
where $`W(𝐫,𝐫^{};\omega )`$ is the frequency-dependent dynamically screened interaction

$`W(𝐫,𝐫^{};\omega )=`$ $`v(𝐫,𝐫^{})+{\displaystyle 𝑑𝐫_1𝑑𝐫_2v(𝐫,𝐫_1)}`$ (6)
$`\times \chi (𝐫_1,𝐫_2;\omega )v(𝐫_2,𝐫^{}),`$
$`v(𝐫,𝐫^{})`$ and $`\chi (𝐫_1,𝐫_2;\omega )`$ being the bare Coulomb interaction and the exact density-response function of the Fermi system, respectively.

In the case of a bounded 3D electron gas that is translationally invariant in the plane of the surface, which we assume to be normal to the $`z`$ axis, the probe-particle initial and final states are of the form

$$\psi _i(𝐫)=\varphi _i(z)\mathrm{e}^{i𝐤_{}𝐫}$$
(7)
and

$$\psi _f(𝐫)=\varphi _f(z)\mathrm{e}^{i(𝐤_{}𝐪_{})𝐫},$$
(8)
with energies

$$E_i=\epsilon _i+\frac{𝐤_{}^2}{2m_i}$$
(9)
and

$$E_f=\epsilon _f+\frac{(𝐤_{}𝐪_{})^2}{2m_f},$$
(10)
where the one-particle wave functions $`\varphi (z)`$ and energies $`\epsilon `$ describe motion normal to the surface. Using these wave functions and energies in Eq. (3), one easily finds

$`\tau ^1=`$ $`2{\displaystyle \underset{f}{}}{\displaystyle \frac{\mathrm{d}𝐪_{}}{(2\pi )^2}dzdz^{}\varphi _i^{}(z)\varphi _f^{}(z^{})}`$ (13)
$`\times \mathrm{Im}\left[W(z,z^{};𝐪_{},|E_iE_f|)\right]\varphi _i(z^{})\varphi _f(z),`$
$`W(z,z^{};𝐪_{},\omega )`$ being the 2D Fourier transform of the screened interaction.

## 3 Screened interaction

The main ingredient in the evaluation of electron-hole decay rates in a bounded 3D electron gas is the screened interaction $`W(z,z^{};𝐪_{},\omega )`$. When both $`z`$ and $`z^{}`$ are fixed far from the surface, a few atomic layers within the bulk, there is translational invariance in the direction normal to the surface, and $`W(z,z^{};𝐪_{},\omega )`$ can then be easily obtained from the knowledge of the dielectric function of a homogeneous electron gas $`ϵ(q,\omega )`$:

$$W^{bulk}(z,z^{};𝐪_{},\omega )=\frac{dq_z}{2\pi }\mathrm{e}^{iq_z(zz^{})}v(q)ϵ^1(q,\omega ),$$
(14)
where $`v(q)`$ represents the three-dimensional Fourier transform of the bare Coulomb interaction, and $`q^2=q_{}^2+q_z^2`$. In the random-phase approximation (RPA), $`ϵ(q,\omega )`$ is the Lindhard dielectric function.

Shockley surface-state electrons and holes are found to be very close to the Fermi level ($`|EE_F|<<E_F`$). Hence, as the energy transfer $`|E_iE_f|`$ entering Eq. (13) cannot exceed the value $`|EE_F|`$, one can use the low-frequency form of the Lindhard dielectric function, which yields the energy-loss function

$$\mathrm{Im}\left[ϵ^1(q,\omega )\right]=\frac{2}{q^3}ϵ^2(q,0)\omega \mathrm{\Theta }(2q_Fq),$$
(15)
where $`\mathrm{\Theta }(x)`$ represents the Heaviside function and $`q_F`$ is the Fermi momentum. If one further replaces the static dielectric function $`ϵ(q,0)`$ by the Thomas-Fermi (TF) approximation, then one finds

$$\mathrm{Im}\left[ϵ^1(q,\omega )\right]=\frac{2q}{(q^2+q_{TF}^2)^2}\omega \mathrm{\Theta }(2q_Fq),$$
(16)
$`q_{TF}=\sqrt{4q_F/\pi }`$ being the TF momentum.

In the bulk, the maximum magnitude of $`\mathrm{Im}\left[W(z,z^{};𝐪_{},\omega )\right]`$ occurs at $`z=z^{}`$ and is independent of the actual value of $`z`$. Fig. 2 shows this quantity for Cu(111), as obtained from Eq. (14) with either the Lindhard dielectric function or the approximated energy-loss function of Eq. (16), versus $`q_{}`$ for fixed values of $`\omega `$ \[$`\omega =0.2,0.3,0.4,0.5\mathrm{eV}`$\]. One clearly sees that for the momentum and energy transfers of interest (see Fig. 1), the approximated form of Eq. (16) yields a screened interaction which is close to that obtained with use of the Lindhard dielectric function, small differences being mainly due to the fact that the static dielectric function $`ϵ(q,0)`$ has been replaced by the TF approximation. Contributions to the surface-state decay rate coming from the penetration of the surface-state wave function into the solid are, therefore, expected to be well described by Eq. (16).

Nevertheless, coupling of Shockley surface states with the crystal may also occur, either through the evanescent tails of bulk states near the surface or within the surface-state band itself. These contributions to the surface-state decay rate are both dictated by the screened interaction at $`z`$ points near the surface, where the representation of Eq. (14) is not accurate.

For a realistic description of the screened interaction, we proceed along the lines reported in . First of all, we compute the single-particle eigenfunctions and eigenvalues of a 1D model potential. This model potential, which reproduces far outside the surface the classical image potential, is chosen so as to describe the width and position of the energy gap at the $`\mathrm{\Gamma }`$ point ($`𝐤_{}=0`$) and, also, the binding energies of both the Shockley surface state at $`\mathrm{\Gamma }`$ and the first ($`n=1`$) image-potential induced state. We then introduce these eigenfunctions and eigenvalues into Eq. (13), we also use them to compute the non-interacting density-response function, and finally solve an integral equation to derive the RPA density-response function $`\chi (z,z^{};𝐪_{},\omega )`$ and the screened interaction $`W(z,z^{};𝐪_{},\omega )`$.

As the maximum magnitude of $`\mathrm{Im}\left[W(z,z;𝐪_{},\omega )\right]/\omega `$ near the surface still occurs at $`zz^{}`$ (only far from the surface into the vacuum, where the Shockley surface-state amplitude is negligible, does the maximum of this quantity occur at $`zz^{}`$), we choose $`z=z^{}`$ and show $`\mathrm{Im}\left[W(z,z;𝐪_{},\omega )\right]/\omega `$ in Fig. 3, as a function of $`z`$ going from the bulk to the vacuum, for $`q_{}=0.2a_0^1`$ ($`a_0`$ is the Bohr radius) and for various values of $`\omega `$. We find that at the surface $`\mathrm{Im}[W]/\omega `$ is enhanced, for the smallest values of $`\omega `$, by a factor of 3 relative to the bulk. This is mainly due to the strong anisotropy at the surface, which enhances the electron-hole pair creation probability. We also find that for frequencies smaller than $`0.20.3\mathrm{eV}`$, the strong anisotropy of the surface does not modify the linear frequency scaling characteristic of the bulk screened interaction \[see Eq.(15)\]. However, for the larger but still small frequencies explored in Fig. 3, the screened interaction at the surface exhibits an important enhancement due to corrections to the linear scaling that are not present within the bulk. As a result, for the largest values of $`\omega `$ that we have considered, $`\mathrm{Im}[W]/\omega `$ is enhanced at the surface by a factor as large as 4.

Fig. 4 displays $`\mathrm{Im}\left[W(z,z;𝐪_{},\omega )\right]/\omega `$ as a function of $`z`$, for $`\omega =0.2\mathrm{eV}`$ and for various values of $`q_{}`$. While $`\mathrm{Im}[W]/\omega `$ is enhanced at the surface, for the smallest values of $`q_{}`$, by a factor of 3 or 4 in the frequency range $`\omega =0.20.5\mathrm{eV}`$, this enhancement becomes very small for the largest values of $`q_{}`$ that we have explored. These calculations indicate that only transitions with $`q<0.4a_0^1`$ are affected by the strong anisotropy at the surface.

## 4 Lifetime broadening

We compute the lifetime broadening of Shockley surface states from Eq. (13) with a realistic description of the RPA screened interaction, as described in the previous section. Although Eq. (13) has been derived by assuming that our electron system is translationally invariant in the plane of the surface, we account for the potential variation parallel to the surface through the introduction of a realistic effective mass into Eqs. (9) and (10), and also through the introduction into Eq. (13) of initial and final wave functions that change along the actual dispersion curve of each state. The effective mass of bulk states has been chosen to increase from our computed values of 0.31, 0.25, and $`0.21m_e`$ at the bottom of the gap in Cu(111), Ag(111), and Au(111), respectively, to the free-electron mass $`m_e`$ at the bottom of the valence band. The $`z`$-dependent initial and final wave functions $`\varphi _i(z)`$ and $`\varphi _f(z)`$ have been recalculated for each value of $`𝐤_{}`$ and $`𝐤_{}𝐪_{}`$, respectively, as was done at the $`\mathrm{\Gamma }`$ point, with use of the 1D hamiltonian described in the previous section.

As the relaxation of Shockley surface-state electrons and holes may proceed either through the interband coupling with bulk states (3D channel) or through the intraband coupling, within the surface state itself, with surface states of different wave vector parallel to the surface (2D channel), we also consider the decay of electrons and holes in 3D and 2D uniform systems. Introduction of Eq. (16) into Eq. (14) and then Eq. (14) into Eq. (13) with both $`\varphi _i(z)`$ and $`\varphi _f(z)`$ replaced by plane waves, yields

$$\tau _{3D}^1=\frac{\sqrt{\pi q_F}}{8q_F^2}\left[\mathrm{tan}^1\sqrt{\pi q_F}+\frac{\sqrt{\pi q_F}}{1+\pi q_F}\right]\frac{(EE_F)^2}{k},$$
(17)
$`k`$ being the momentum of the excited electron or hole. On the same level of approximation, the 2D decay rate is given by,

$$\tau _{2D}^1=\frac{E_F}{4\pi }\left[\mathrm{ln}\frac{|EE_F|}{E_F}+\frac{1}{2}+\mathrm{ln}\frac{2q_{2D}^{TF}}{q_F}\right]\left(\frac{EE_F}{E_F}\right)^2,$$
(18)
where $`q_{2D}^{TF}=2m`$ is the TF screening wave vector in 2D, and $`m`$ is the electron mass.

### A Surface-state holes at $`\mathrm{\Gamma }`$

The partially occupied Shockley surface-state band forms a uniform two-dimensional electron gas, with the 2D Fermi energy being given by the band-edge of the surface state, i.e., $`E_F`$=445, 67, and 505 eV for the (111) surfaces of Cu, Ag, and Au, respectively. Our full calculation for the intraband linewidth of the band-edge surface-state hole on the (111) surfaces of Cu, Ag, and Au is presented in Fig. 5 by full squares, together with the linewidth of Eq. (18) with $`E=0`$ and the electron mass $`m`$ chosen to be either the free-electron mass (solid line) or the surface-state effective mass (full circles). We find that calculations based on a pure 2D electron gas are larger than the actual intraband contribution to the linewidth by a factor of $`7`$. This large discrepancy is due to the fact that electron-electron interactions within the actual Shockley surface-state 2D band are strongly screened by the underlying 3D bulk electron system, thereby reducing the scattering probability.

Separate intraband and interband contributions to the linewidth of Shockley surface-state holes at the $`\mathrm{\Gamma }`$ point of the projected bulk band-gap of the (111) surfaces of Cu, Ag, and Au are displayed in Table I. The 3D linewidth of bulk free holes with the same energy, as obtained from Eq. (17), is also shown in this table. Differences between our full interband calculations and those obtained from Eq. (17) arise from (a) the enhancement of $`\mathrm{Im}[W]`$ at the surface, which increases the linewidth, (b) localization of the surface-state wave function in the direction perpendicular to the surface, and (c) the restriction that only bulk states with energy lying below the projected band-gap are allowed. Both localization of the surface-state wave function and the presence of the band-gap reduce the linewidth, and therefore tend to compensate the enhancement of $`\mathrm{Im}[W]`$ at the surface. In the case of Cu(111) this compensation is nearly complete, thereby yielding an interband linewidth that nearly coincides with the 3D linewidth of free holes. However, this is not necessarily the case for other materials, depending on the surface band structure.

The impact of the enhanced $`\mathrm{Im}[W]`$ at the surface on both interband and intraband contributions to the Shockley surface-state hole linewidth at the $`\mathrm{\Gamma }`$ point in Cu(111) is illustrated in Fig. 6. In this figure, our full calculations of the interband and intraband contributions to the total inelastic linewidth (full circles) are compared to the results we have also obtained from Eq. (13), but with the actual screened interaction replaced by that of Eq. (14) (full triangles). 3D and 2D free-electron gas calculations, as obtained from Eqs. (17) and (18), are represented by open squares. We find that the impact of the surface anisotropy on $`\mathrm{Im}[W]`$ (difference between full circles and triangles) is to largely increase both interband and intraband contributions to the linewidth. Fig. 6 clearly shows that the agreement, in the case of Cu(111), between our interband linewidth (full circle) and that obtained from Eq. (17) (open square) is due to a fortuitous cancellation between both localization of the surface-state wave function and the presence of the band-gap, on the one hand, and the enhancement of $`\mathrm{Im}[W]`$ at the surface, on the other. We also note from this figure that a large contribution from intraband transtions, strongly screened by the underlying 3D bulk electron system, is responsible for the large differences between the 3D free-electron gas prediction and the experimental results (represented by open triangles), as discussed in .

### B Surface-state electrons and holes with $`𝐤_{}0`$

Our full calculation of the total inelastic linewidth of Shockley surface-state electrons and holes on Cu(111), as obtained from Eq. (13), is shown in Fig. 7 (solid line with circles) versus the surface-state energy. Separate interband and intraband contributions to the linewidth are also represented in this figure, by solid lines with triangles (interband) and inverted triangles (intraband), and the linewidth of bulk free-electrons with the same energy, as obtained from Eq. (17), is represented by a dotted line. Measurements of the linewidth of excited holes at the surface-state band edge ($`𝐤_{}=0`$) are represented by open triangles, and the STS measurements for surface-state electrons reported by Burgi et al ($`𝐤_{}0`$) are represented by full squares.

We find that the fortuitous agreement between the actual interband linewidth and that derived from Eq. (17) is present not only for surface-state holes at the $`\mathrm{\Gamma }`$ point, but also for surface-state electrons and holes near the Fermi level. Hence, the interband linewidth shows the $`(EE_F)^2`$ energy dependence that is present in the case of a 3D free-electron gas.

Fig. 7 shows that at the surface-state band edge ($`𝐤_{}=0`$ the intraband linewidth represents an $`80\%`$ of the total linewidth (see also Table I). Nevertheless, as the surface-sate wave vector parallel to the surface increases, the surface-state wave function acquires a bulk-like character, with a larger penetration into the bulk, and the intraband contribution to the linewidth decreases very rapidly. This conclusion explains the experimental observation that while 3D free-electron gas predictions of surface-state hole lifetimes are too large, they are in the case of surface-state electrons (above the Fermi level) comparable to measured lifetimes.

Open triangles in Fig. 7 represent the inelastic contribution to measured linewidths of Shockley surface-state holes at the $`\mathrm{\Gamma }`$ point, which are to be be compared with our calculated linewidth at this point (full circle). Since our model, which correctly reproduces the behaviour of the $`s`$-$`p`$ valence states, does not include screening from $`d`$-electrons, differences between measured linewidths (open triangles) and our calculations (full circles) are expected to be due to the effect of virtual transitions, giving rise to additional screening by $`d`$-electrons, as discussed in .

The STS measurements represented in Fig. 7 by full squares include both inelastic and electron-phonon contributions to the surface-state linewidth. Hence, one must be cautious in the comparison of these measurements with our calculations. Electron-phonon linewidths of typically $`8.0`$ in Cu(111), which are essentially independent of the surface-state energy, yield theoretical predictions for the total linewidth that are above the experimental observation, especially for the lowest surface-state energies explored. This discrepancy between theory and experiment could be partially compensated by the screening of $`d`$-electrons, which as in the case of surface-state holes at the $`\mathrm{\Gamma }`$ point, reduce the inelastic linewidth.

## 5 Summary and conclusions

We have reported theoretical investigations of surface-state electron and hole lifetimes in the noble metals Cu, Ag, and Au. We have also presented self-consistent calculations of the screened interaction, and have investigated the role that the partially occupied 2D surface-state band and the underlying 3D electron gas play in the relaxation mechanism.

We have reached the conclusion that, while the linewidth of surface-state excited holes at the $`\mathrm{\Gamma }`$ point of the (111) surfaces of the noble metals is dominated by a 2D decay channel, major contributions to the electron-electron interaction of surface-state electrons above the Fermi level come from the underlying bulk electrons. This key dependence of the relative contribution of intraband transitions to the total decay rate explains the experimental observation that while 3D free-electron gas predictions of surface-state hole lifetimes are too large, they are in the case of surface-state electrons (above the Fermi level) comparable to measured lifetimes.

###### Acknowledgements.

We acknowledge partial support by the University of the Basque Country, the Basque Hezkuntza, Unibertsitate eta Ikerketa Saila, and the Spanish Ministerio de Educación y Cultura and the Max Planck Research Award funds.


text:
# Long-range effects in granular avalanching

## Abstract

We introduce a model for granular flow in a one-dimensional rice pile that incorporates rolling effects through a long-range rolling probability for the individual rice grains proportional to $`r^\rho `$, $`r`$ being the distance traveled by a grain in a single topling event. The exponent $`\rho `$ controls the average rolling distance. We have shown that the crossover from power law to stretched exponential behaviors observed experimentally in the granular dynamics of rice piles can be well described as a long-range effect resulting from a change in the transport properties of individual grains. We showed that stretched exponential avalanche distributions can be associated with a long-range regime for $`1<\rho <2`$ where the average rolling distance grows as a power law with the system size, while power law distributions are associated with a short range regime for $`\rho >2`$, where the average rolling distance is independent of the system size.

There is nowadays a massive evidence of scale invariant phenomena in nature. They appear in systems as diverse as geological (e.g., earthquakes), climatic (e.g., atmospheric turbulence), granular flow (e.g., rice piles, the topic of this work) and biological (e.g., biological evolution, cell growth, among many others). In many of these cases no particular tunning is perceived. In 1987 Per Bak and colaborators  advanced the hypothesis that this is so because the microscopic dynamics of the system makes it to spontaneously evolve towards a critical, scale-invariant state. This is known today as self-organized criticality (SOC). To illustrate the idea they introduced a sand pile model  that quickly became the paradigm ofSOC models. In this model a pile is gradually built by adding individual grains one by one into an open system. After a transient period the pile reaches a stationary critical state, where the slope of the pile fluctuates around a stationary value. At this stage avalanches of all possible sizes are observed, giving a power law avalanche size distribution. This is the SOC state. The experimental search for power laws in the avalanche size distribution in real granular systems presented a great challenge . Some years ago the Oslo group  found evidence of SOC in controlled experiments on the granular dynamics of rice piles. In 1996 Frette et al  performed an experiment where elongated rice grains were added slowly in a narrow gap between two plates. They found that the avalanche size distribution for grains with large aspect ratio presents a power law behavior, while a stretched exponential behavior is observed for rounder grains. These experimental results showed that SOC is not insensitive to the details of the system such as the shape of the grains and that inertia effects play an important role in the relaxation dynamics. The Oslo group introduced a “rice-pile” model to phenomenologically describe their experiments . This model describes the phenomenology of rice piles presenting SOC behavior. Moreover, further studies showed that it can be exactly mapped to a model for interface depinning . However, the Oslo model does not present a stretched exponential distribution in any region of its parameter space.

In a recent communication Head and Rodgers  presented a clustered based model for granular flow, which exhibits both stretched exponential and power law behavior over its parameter range. Through this model they showed that power law behavior appears to be related to the coherent motion of clusters of grains that slide along the surface, while the stretched exponential behavior results from a fluid like motion associated with the rolling of isotropic grains. However, the motion of individual grains is implicitly incorporated in this model, where the basic dynamical variables are related to cluster properties. Hence, it is interesting to investigate in a more microscopic level how the crossover from power law to stretched exponential behavior emerges as a consequence of a change in the transport properties of individual grains.

In this work we present a model where the ability of individual grains to roll a distance $`r`$ is described by a long-range rolling probability of the form

$$P(r)=\{\begin{array}{cc}A/r^\rho \hfill & if\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}rN\hfill \\ 0\hfill & otherwise\hfill \end{array}$$
(1)
where $`A`$ is a normalization constant, $`N`$ is the system size, $`r=1,2,\mathrm{},N`$ and $`0\rho \mathrm{}`$. The parameter $`\rho `$ controls the average rolling distance which is expected to depend on the aspect ratio of the grains. The $`\rho \mathrm{}`$ limit corresponds to a nearest neighbors movement ($`P(r)=\delta _{r,1}`$), thus describing the case where the grains do not roll. In the opposite limit $`\rho =0`$ the grains can move to any site of the pile (or even drop off in a single jump) with equal probability. Hence, this limit can be associated to the idealized situation of perfect spherical grains without friction. Moreover, studies on different dynamical systems with similar long-ranged interactions have shown that, except from some possible rescaling factor, the behavior of almost all the relevant properties of the systems in the full range $`0\rho <1`$ reproduce the corresponding ones for $`\rho =0`$. This appears to be related to the behavior of the first moment of the distribution (1), which in the present context represents the mean rolling distance of the grains. For large values of $`N`$ we have that

$$r\frac{1\rho }{2\rho }\frac{N^{2\rho }1}{N^{1\rho }1}\{\begin{array}{cc}\frac{(1\rho )}{(2\rho )}N& \text{if }\rho <1\hfill \\ & \\ \frac{N}{ln(N)}& \text{if }\rho =1\hfill \\ & \\ \frac{(\rho 1)}{(2\rho )}N^{2\rho }& \text{if }1<\rho <2\hfill \\ & \\ \frac{ln(N)}{N}& \text{if }\rho =2\hfill \\ & \\ \frac{(1\rho )}{(2\rho )}& \text{if }\rho >2\hfill \end{array}$$
(2)
We see that, for $`0\rho <1`$, $`r`$ is of the order of the system size, which means that a very large number of grains will drop off of the system in a single jump, for any system size. Since this behavior appears to be highly unrealistic we will concentrate our study on the $`\rho 1`$ regime .

Our model is based on the Oslo one  and it is defined as follows. We consider a one dimensional lattice of size $`N`$ ($`1iN`$), each site $`i`$ having associated an integer variable $`h(i)`$ representing the local height of the pile. The local slope is then given by $`\sigma (i)=h(i)h(i+1)`$. The grains enter into the system from the left ($`i=1`$) and may drop off at the rightmost site $`i=N+1`$, imposing $`h(N+1)=0`$ for all times. Every time the local slope $`\sigma (i)`$ of a site $`i`$ exceeds a local critical value $`\sigma (i)>\sigma _c(i)`$ the topmost grain at site $`i`$ rolls $`r`$ sites to the right with probability $`P(r)`$ given by Eq.(1). Then, the heights of the sites $`i`$ and $`i+r`$ are recalculated as $`h(i)h(i)1`$ and $`h(i+r)h(i+r)+1`$ and the corresponding local slopes are modified accordingly. Each time a grain leaves a column $`i`$ we assign ita new critical slope $`\sigma _c(i)=1`$ or $`\sigma _c(i)=2`$ with equalprobability. This process is repeated until all the local slopes satisfy $`\sigma (i)\sigma _c(i)`$. An avalanche starts when $`\sigma (1)>\sigma _c(1)`$ and when it stops (i.e., $`\sigma (i)\sigma _c(i)`$ $`i`$) new grains are added until a new avalanche is initiated. Notice that in the limit $`\rho \mathrm{}`$ the present model reduces to Oslo model.

Once the system attains the stationary state, fluctuations in the slope of the pile appear in the form of avalanches that redistribute the grains, thereby changing the profile. By considering the profile just before an avalanche is initiated and the final profile when it stops, the avalanche size is defined as the total energy dissipated between both profiles . The local energy change in a site $`i`$ is calculated as the difference between the potential energies in both profiles, where the units were chosen such that $`mg=1`$. We study for different values of $`\rho `$ the probability density $`P(E,N)`$, $`P(E,N)dE`$ being the probability that an avalanche with energy dissipation between $`E`$ and $`E+dE`$ will occur in a system of size $`N`$. Numerical data were smoothed using a local average procedure in order to diminish large statistical fluctuations present for high values of the ratio $`E/N`$. We verified for a wide range of values of $`\rho `$ (even in the non physical region $`\rho <1`$) that $`P(E,N)`$ displays the finite-size scaling behavior $`P(E,N)=N^1f(E/N)`$, in agreement with the experimental data , as shown by the data collapse of Fig. 1 for $`\rho =1.4`$ and $`\rho =4`$. The scaling function $`f(x)`$ depends on the value of $`\rho `$. We found two clearly distinct regimes. For $`\rho >2`$ (hereafter referred as the “sliding regime”) $`f(x)`$ is almost constant for small values of $`x`$ and presents a power law dependency $`f(x)x^\alpha `$ for large values of $`x`$, as shown in the example of Fig.2 for $`\rho =4`$. The exponent $`\alpha `$ increases smoothly from $`\alpha 1`$ in the limit $`\rho 2^+`$ to $`\alpha =1.33\pm 0.05`$ for values of $`\rho 4`$, where it becomes independent of $`\rho `$. This value is consistent with that obtained by Head and Rodgers  for the cluster based rice pile model.

For $`1<\rho <2`$ (hereafter referred as the “rolling regime”) the best fitting of the numerical data is obtained by a function of the form:

$$f(x)\frac{1}{x^\alpha }\mathrm{exp}\left(\left(\frac{x}{x_0}\right)^\gamma \right)$$
(3)
as shown in the example of Fig.3 for $`\rho =1.2`$. The parameters $`\gamma `$, $`\alpha `$ and $`x_0`$ depend on $`\rho `$, as shown in table I, although the small variation of $`\gamma `$ in this regime suggests that it may not depend on $`\rho `$. Notice that this function seems to describe the behavior of the whole distribution, not just for large values of $`E/N`$ (where it is dominated by the stretched exponential) but also for small values of it.

We have shown that the crossover from power law to stretched exponential behavior observed in the granular dynamics of rice piles can be well described as a long-range effect resulting from a change in the transport properties of individual grains. Within this scenario two distinct regimes appear, regarding the qualitative behavior of the avalanche distribution: a short-range “sliding” regime and a long-range “rolling” one. A direct inspection of snapshots of the profiles obtained during the simulations shows that, indeed grains group into clusters that move coherently in the sliding regime, while in the rolling regime they move almost independently. These behaviors are consistent with what is observed experimentally  and with Head and Rodgers results . Details of these simulations will be published elsewhere. This change in the transport properties of individual grains is related to different scaling properties of the average rolling distance with the system size: while in the rolling regime $`d`$ scales as a power law $`dN^{2\rho }`$, in the sliding regime $`d`$ does not depend on $`N`$.

In the whole sliding regime $`\rho >2`$ the system displays the qualitative behavior expected for a short-range model, that is, a power law for large values of the energy density $`E/N`$ and constant for small values of it. However, it is worth noting that for intermediate values of $`\rho `$ ($`2<\rho <4`$) the characteristic exponent $`\alpha `$ of the power law distributions depends on $`\rho `$, while for large values of $`\rho `$ it becomes independent of $`\rho `$. This results suggests that non-universal behavior may be expected for intermediate values of the grains aspect ratio. The value $`\rho =4`$ is only indicative and it seems not to be a critical value of $`\rho `$ separating the universal from the non-universal behaviors.

In the rolling regime the whole avalanche distribution is well described by a stretched exponential modulated by a power law. This correction to the usual stretched exponential behavior is only important in the small energy density region of the spectrum, since for $`xx_0`$ the function is always dominated by the stretched exponential. Notice that for $`1.4<\rho <2`$ the exponent $`\alpha `$ changes of sign taking small absolute values, suggesting logarithmic rather than power law corrections to the stretched exponential in this region. Indeed, a slightly better fitting is obtained with a function of the form

$$f(x)\mathrm{ln}(x)\mathrm{exp}\left(\left(\frac{x}{x_0}\right)^\gamma \right)$$
(4)
without any significant change in the fitting parameter values of $`x_0`$ and $`\gamma `$ shown in table I.

From table I we see that the best agreement  with the experimental values $`x_0=0.45\pm 0.09`$ and $`\gamma =0.43\pm 0.03`$ is obtained for $`\rho =1.01`$. Both functions Eq.(3) and (4) cease to give good fittings of the simulation data as we approach the border values $`\rho 1^+`$ and $`\rho 2^{}`$ respectively. This fact, together with an increase in the statistical fluctuations in the neighborhood of those values, makes very difficult to determine how the crossover from stretched exponential and power law occurs around $`\rho =2`$. Possibly the calculation of the distributions of other quantities (like transit times) may be useful. Some work along this line is in progress.

This work was partially supported by grants from Consejo Nacional de Investigaciones Científicas y Técnicas CONICET (Argentina), Consejo Provincial de Investigaciones Científicas y Tecnológicas (Córdoba, Argentina) and Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (Argentina).