text:
# Nonextensivity of the cyclic Lattice Lotka Volterra model

## Abstract

We numerically show that the Lattice Lotka-Volterra model, when realized on a square lattice support, gives rise to a finite production, per unit time, of the nonextensive entropy $`S_q=\frac{1_ip_i^q}{q1}`$ $`(S_1=_ip_i\mathrm{ln}p_i)`$. This finiteness only occurs for $`q=0.5`$ for the $`d=2`$ growth mode (growing droplet), and for $`q=0`$ for the $`d=1`$ one (growing stripe). This strong evidence of nonextensivity is consistent with the spontaneous emergence of local domains of identical particles with fractal boundaries and competing interactions. Such direct evidence is for the first time exhibited for a many-body system which, at the mean field level, is conservative.

Keywords: Lattice-Lotka-Volterra model, Fractals, Reaction limited processes, Nonextensivity.

Pacs Numbers: 5.40.-a, 05.10.Ln, 05.65.+b, 82.40.Np

Many natural and artificial systems are known today to be hardly, or not at all, tractable within Boltzmann-Gibbs statistical mechanics, hence the usual thermodynamics. Such is the case of systems which include long-range interactions or long-range microscopic (or mesoscopic) memory, or other sources of (multi)fractality. Phenomena where many spatial and/or temporal scales are involved typically exhibit power laws. The celebrated Boltzmann-Gibbs (BG) entropy $`S_{BG}=_ip_i\mathrm{ln}p_i`$ appears to be inadequate for handling the thermostatistics associated with such situations. This is due to the fact that the corresponding stationary states do not emerge through ergodic dynamics. An ubiquitous class of the above anomalous systems have nonlinear dynamics which generate weak chaos, in the sense that the sensitivity to the initial conditions is less-than-exponential in time. Such situations quite naturally accomodate with an entropy which generalizes the BG one, namely

$$S_q=\frac{1_ip_i^q}{q1}(q;S_1=S_{BG}).$$
(1)
For independent systems $`A`$ and $`B`$ (i.e., such that $`p_{ij}^{A+B}=p_i^Ap_j^B`$), this entropy satisfies $`S_q(A+B)=S_q(A)+S_q(B)+(1q)S_q(A)S_q(B)`$. It is due to this property of nonextensivity that the thermostatistical formalism based on Eq. (1) is usually referred to as nonextensive statistical mechanics  (see  for recent reviews). This theory has received many applications in areas such as self-gravitating polytropes , electron-positron annihilation , turbulence , motion of Hydra viridissima , anomalous diffusion , classical  and quantum  chaos, long-range-interacting many-body Hamiltonians , option pricing , particular biological processes involving a large number of degrees of freedom , among others.

The lattice Lotka-Volterra (LLV) model  is known to well mimic simple chemical reactions, predator-prey systems, and other biological and echological phenomena. It has recently been studied with cyclic interactions amongst three species, as a modification to the original Lotka-Volterra model. The LLV has exhibited spatial clustering and fractality, even when implemented on low dimensional lattices . Since fractality is a distinct sign of nonextensivity, it is natural to pose the question whether the LLV model is indeed consistent with the nonextensive premises. The aim of the present paper is to verify that this model exhibits a strong and direct evidence of having $`q1`$. In particular we study the LLV model with various initial conditions, namely at a) the domain formation mode, b) the nucleus (or droplet) growth mode, c) the stripe growth mode and d) the roughening mode. Modes (b) and (c) enable, as we shall see, the direct calculation of $`q`$. Several analytical or numerical calculations of $`q`$ exist already in the literature, but this is the first time such evidence is directly provided, through the time evolution of $`S_q`$ itself, on a many-body system which is conservative at the Mean Field (MF) level.

The LLV model is a minimal complexity model, with MF conservative dynamics which can be directly implemented on lattice and involves only two reactive species $`X_1`$ and $`X_2`$ (adsorbed on a lattice support) and the empty sites of the support $`S`$. All reactive steps are bimolecular and the reaction occurs via hard core interactions. Schematically, the LLV model has the following form:

$`X_1+X_2`$ $`\stackrel{k_s}{}2X_2`$ (3)
$`X_2+S`$ $`\stackrel{k_1}{}2S`$ (4)
$`S+X_1`$ $`\stackrel{k_2}{}2X_1`$ (5)
In particular, a particle $`X_1`$ adsorbed on a lattice site changes its state into $`X_2`$ when it is found in the neighborhood of another $`X_2`$ particle. This step (2a) is an autocatalytic reactive step. A particle $`X_2`$ desorbs leaving an empty site $`S`$, if in the neighborhood another empty site $`S`$ is found. This step (2b) is a cooperative desorption step. Finally, a particle $`X_1`$ can be adsorbed on an empty lattice site $`S`$ if in the neighborhood another $`X_1`$ particle is found. This step (2c) is a cooperative adsorption step.

We now recall briefly some of the mean field (MF) and lattice properties of the LLV, which have been studied in detail in previous works.

In the MF approximation the LLV model, Eqs. (1), can be described via the kinetic, rate equations:

$`{\displaystyle \frac{dx_1}{dt}}`$ $`=`$ $`x_1(k_sx_2+k_2s)`$ (4)
$`{\displaystyle \frac{dx_2}{dt}}`$ $`=`$ $`x_2(k_sx_1k_1s)`$ (5)
$`{\displaystyle \frac{ds}{dt}}`$ $`=`$ $`s(k_2x_1+k_1x_2)`$ (6)
where $`x_1`$, $`x_2`$ and $`s`$ correspond to the mean coverage of the lattice with particles $`X_1`$, $`X_2`$ and empty sites $`S`$ respectively. In Eqs. (3), the mean coverages satisfy identically the conservation condition $`x_1+x_2+s=C`$, where $`C`$ is a constant which can be chosen equal to unity, corresponding to interpreting $`x_1`$, $`x_2`$ and $`s`$ as fractions of the overall lattice respectively occupied by $`X_1`$ particles, $`X_2`$ particles or being empty. Using $`C=1`$ it is possible to eliminate one of the three variables, say $`s=1x_1x_2`$, and, to reduce system (3) to two equations. This reduced system admits four steady state solutions, three of which are trivial, and one non-trivial:

$`x_{1s}`$ $`=`$ $`0,x_{2s}=0(\mathrm{empty}\mathrm{lattice})`$ (5)
$`x_{1s}`$ $`=`$ $`1,x_{2s}=0(\mathrm{lattice}\mathrm{poisoned}\mathrm{by}\mathrm{X}_1)`$ (6)
$`x_{1s}`$ $`=`$ $`0,x_{2s}=1(\mathrm{lattice}\mathrm{poisoned}\mathrm{by}\mathrm{X}_2)`$ (7)
$`x_{1s}`$ $`=`$ $`{\displaystyle \frac{k_1}{k_1+k_2+k_s}},x_{2s}={\displaystyle \frac{k_2}{k_1+k_2+k_s}}`$ (8)
A linear stability analysis shows that the trivial states are saddle points while the nontrivial one is a center compatible with an additional constant of motion $`C^{}=x_1^{k_1}x_2^{k_2}(1x_1x_2)^{k_s}`$ at the MF level. Fig. 1a depicts the temporal evolution of the system for typical values for $`(k_1,k_2,k_s)`$ and initial conditions. The black solid line represents the concentration of $`X_1`$ and the dashed line the concentration of $`X_2`$. The motion is periodic but non-harmonic. The amplitude of the periodic motion, for given parameter values, depends solely on the initial conditions. At this level of description the system size does not enter into the calculations since the MF approximation involves only average concentrations.

To mesoscopically describe the system on a lattice, many details enter: lattice size and geometry, number of nearest neighbors (coordination number), interaction range, etc. To realize the square lattice LLV we adopt a typical Monte Carlo (MC) algorithm (details in ), namely (1) At every microscopic step one lattice site is randomly chosen; (2) One of the nearest neighbors is also selected randomly; (3) If the original chosen site is $`X_1`$ ($`X_2`$) and the selected neighbor is $`X_2`$ ($`S`$) then the chosen site changes to $`X_2`$ ($`S`$) with probability $`k_s`$ ($`k_1`$); If the original chosen site is $`S`$ and the selected neighbor is $`X_1`$ then the chosen site changes to $`X_1`$ with probability $`k_2`$; Otherwise the system remains as it is; (4) Return to step 1.

In the MC procedure the unit of time is chosen as $`1/N`$, where $`N`$ is the total number of lattice sites (occupied and empty). For example, for square lattice, $`N=L^2`$, where $`L`$ is the linear size of the lattice. With this choice of micro-time, in one MC step all lattice sites are, on the average, scanned once. In Fig. 1b typical behavior of the temporal evolution of the MC concentrations is shown. In particular, the concentrations of $`X_1`$ is depicted on the full lattice of size $`L\times L=2^8\times 2^8`$ (solid line) and on a sublattice of size $`l\times l=2^5\times 2^5`$ (dotted line). Periodic boundary conditions are used in all simulations. It is clear that while on the sublattice the concentrations show oscillatory behavior with added noise, on the entire lattice the oscillations shrink. Fig. 2a gives the typical evolution of a system starting from random initial conditions (Fig. 2a (t=0 MC)). As time increases the system develops local domains and each domain behaves as a local oscillator with specific characteristic frequency. Because the various domains have different phases, globally, no oscillation are observed, in contrast with the MF predictions. Moreover, it has been shown that the different species organize in local domains which present competing interactions and they have fractal boundaries. In this figure and hereafter the $`X_1`$ particles are depicted in grey color, the $`X_2`$ in black and the empty sites in white. The initial condition was a homogeneous infinite lattice with equal concentrations of $`X_1,X_2`$ particles and empty sites $`S`$. Periodic boundary conditions are used. The fractal properties of the spatial structures can be used to measure the size of the local oscillators and point out to a nonextensive formalism for the calculation of its entropy.

To describe the temporal evolution of the entropy with respect to one of the species, e.g. $`X_1`$, we start from a given configuration, with specific initial conditions on lattice and let the system evolve according to the MC algorithm. The choice of the particular species does not play any role in the entropy calculations, since the model is cyclic and all the species are equivalent. As time increases the system passes through various configurations which we record at regular temporal intervals. Let us call $`C(t)=\{C_{ij}(t)\},i=1,\mathrm{}L;j=1,\mathrm{}L`$ the specific configuration of the lattice at time $`t`$, while $`C_{ij}(t)`$ denotes the state of site $`(i,j)`$ at time $`t`$ and

$`C_{ij}=\{\begin{array}{cc}1\hfill & \text{if site }(i,j)\text{ is occupied by }X_1\hfill \\ 1\hfill & \text{if site }(i,j)\text{ is occupied by }X_2\hfill \\ 0\hfill & \text{if site }(i,j)\text{ is occupied by }S.\hfill \end{array}`$ (12)
Within each configuration we introduce a set of M non-overlapping windows $`\{W_i\},(i=1,\mathrm{}M)`$, of size $`l\times l`$ which cover completely the lattice. The number of windows is $`M=n^2=(L/l)^2`$. Consequently $`C(t)=_{i=1}^MW_i`$. Let us denote with $`p_i`$ the probability that window $`i`$ is occupied by particles $`X_1`$. If $`n_1(i,t)`$ is the number of particles $`X_1`$ inside window $`i`$ at time $`t`$ and $`n_1(t)`$ is the total number of $`X_1`$ particles on the lattice, then $`p_i(t)=n_1(i,t)/n_1(t)`$. This probability set into Eq. (1) provides $`S_q(t)`$. For short times $`S_q(t)`$ scales as a nonlinear function of the time, while for large times depends non-linearly on the system size. This non-linear dependence on the system size is the basic indication of non-extensivity. The various values of $`q`$ highlight characteristics on different length scales in the system. As an example, rare events are characterized by low values of $`p_i`$. For $`q<1`$, the term $`p_i^q`$ takes relatively large values and gives important contribution to the function $`S_q`$. In contrast, if $`q>1`$, then $`p_i^q<<p_i`$ and the contribution of rare events is negligible.

It is well known that scaling behavior is proper to systems which present fractality. Especially in monofractals only one level of scaling is detected while in multifractal structures the different scales grow with different power laws. The $`S_q`$ entropy is then the appropriate measure of complexity because it addresses the complexity in different length scales by appropriate tuning of the $`q`$ value. We study $`S_q`$ with different initial conditions. Depending on the degree of organization of the initial state, the entropy may increase going to a more disordered state or decrease going to a more ordered state.

Consider first the case of the ”domain growth mode”, where initially the system contains particles $`X_1`$, $`X_2`$ and empty sites $`S`$ randomly distributed and with equal probability, as in Fig 2a (t=0MC). The initial state is statistically uncorrelated. This is the state of maximum entropy. It is easy to calculate $`S_q`$ for this state as a function of the window size. Assume that the $`X_1`$ particles are the information carrying sites, while the $`X_2`$ and $`S`$ are the medium. This assumption does not reduce the system complexity since all three species are equivalent for cyclic reactions. If the three species have the same concentrations then the number of $`X_1`$ particles will be on the average equal to $`x_1(t=0)=L^2/3`$. If we divide the lattice in non-overlapping windows of size $`l\times l`$ then the average number of $`X_1`$ particles in a box is $`l^2/3`$ while the number of windows is $`n^2=L^2/l^2`$. The initial values of $`p_i`$ are:

$`p_i={\displaystyle \frac{l^2/3}{L^2/3}}=\left({\displaystyle \frac{l}{L}}\right)^2(i),`$ (13)
hence $`S_q(0)=[\left(L/l\right)^{2(1q)}1]/(1q)`$. As time increases, the various species organize in domains with fractal boundaries and we call this mode “domain formation mode”. The entropy gradually decreases since the system evolves from a random to a more organized state: See Fig. 2a. In Fig. 2b the corresponding temporal evolution of $`S_q`$ is shown for different values of $`q`$. $`S_q`$ initially undergoes a few oscillations, while the system organizes in domains of homologous (identical) species, and then stabilizes into a lower entropy state.

To calculate the entropy production rate at the nucleus growth mode, we start with a fully organized state consisting only of particles $`S`$ and we include a nucleation droplet of infinitesimal radius $`r`$ placed on the lattice. The droplet contains particles $`X_1`$, $`X_2`$ and $`S`$ homogeneously and randomly distributed within the droplet area. As time increases the droplet grows forming spontaneously several rings of particles $`X_1`$, $`X_2`$ and $`S`$ sequentially. The widths of the rings shrink with the distance from the pure $`S`$ region and when their width becomes zero the typical LLV fractal pattern appears in the middle as can be seen in Fig. 3a. This type of spreading is called the “nucleus growth mode” because an initially small droplet grows in size and finally covers the entire system. This 2-dimensiona ($`d=2`$) growth leads to a reorganization of the species which at the beginning were randomly distributed within the infinitesimal droplet, while they eventually present fractal patterns as in Fig. 2a. As time increases this typical pattern will cover the entire lattice and $`S_q`$ attain the values calculated in the previous case (see Fig. 2b).

As seen in Fig. 3b, only the case $`q=0.5`$ shows a linear increase with time during the entropy production duration; behavior is sub-linear for $`q>0.5`$ and superlinear for $`q<0.5`$. In Fig. 3c we see $`S_{0.5}`$ for various system sizes. They all start linearly, coincide during the entropy production period, and saturate at different values, in accordance with Eq. (5). The $`q`$-value does not change with variations of the lattice size $`L`$, of the window size $`l`$, and of the initial concentrations of reactants within the original droplet. Note the difference in ordinate scale between Figs. 2b and 3b. If we inspect closer the $`L=500`$ steady state of Fig. 3b, the entropy lines for $`t>150`$ present fluctuations similar to the ones in Fig. 2b. An interesting data collapse is shown in Fig. 3d. In Fig. 4a,b the entropy of the LLV model, at the “stripe growth mode” is shown. This is a 1-dimensional growth ($`d=1`$) mode. The initial state of the system consists of a stripe of randomly distributed $`X_1`$, $`X_2`$ and $`S`$ embedded in a lattice containing only $`S`$ particles otherwise (Fig. 4a). The entropy features at this mode are similar to the “nucleus growth mode” but now only $`q=0`$ produces the linear behavior, hence $`q`$ depends on $`d`$.

To explore the entropy production due to surface roughening (or “roughening mode”), we investigate the case of an interface separating two stripes of identical particles. The setup of the system at the initial state is as follows: On a $`500\times 500`$ square lattice we create one stripe of size $`S\times L=50\times 500`$ consisting only by $`X_1`$ particles followed by a stripe of the same size but consisting only of $`X_2`$ particles, while the rest of the lattice is covered by $`S`$, see Fig. 5a (t=0MC). In the Figs. 5a, corresponding to times 0 MC, 20 MC, 60 MC and 100 MC respectively, the originally linear interfaces roughen and the stripes are deformed. The process is dynamical and all interfaces move to the right with the same average velocity. While the size of the stripes is on the average kept constant, fluctuations make them vary significantly. In fact, after sufficiently long times (depending on the width $`S`$ and the size $`L`$ of the stripes) the stripes will mix and the typical fractal patterns of Fig. 2 will reappear. The entropy increase due to the roughening of the interfaces is shown in Fig. 5b. After an initial increase due to roughening the entropy remains constant with statistical fluctuations around the steady state. This mode is clearly not adequate for extracting the physical value of $`q`$.

In the current study the nonextensive entropic properties of the LLV model are examined. The special value of $`q`$ which produces a linear increase of $`S_q(t)`$ depends on the dimensionality $`d`$ of the growth: $`q=11/d(d=1,2)`$. These non-trivial values of $`q`$ are in accordance with the appearance of fractal spatial structures observed in earlier studies for the same system; also, one expects $`q=1`$ in the $`d\mathrm{}`$ limit. Further studies (e.g., sensitivity to the initial conditions, multifractal function $`f(\alpha )`$, entropy relaxation, aging) are welcome.

One of us (CT) acknowledges warm hospitality at, as well as financial support by, the Demokritos Center and the Physics Department of the University of Athens.


format: text

text:
# Dynamics of the 7×7 DAS-reconstructed Silicon (111) Surface

## Abstract

We studied the dynamics of the reconstructed Si (111) surface using a total-energy-vibrational-spectrum approach based on a non-orthogonal tight-binding Hamiltonian. We first established the size of the supercell sufficient to yield a reliable determination of surface parameters by the structural optimization. The site density of vibrational states (SDOS) for the semi-infinite system (the optimized slab on top of the bulk) was then calculated using the method of real space Green s function. A decomposition scheme for the SDOS that identifies directly coupled vibrations for a given mode was also proposed. Our study has uncovered and elucidated, for the first time, all the important surface dynamical features associated with the dimers, adatoms, stacking faults, and rest atoms.

The $`7\times 7`$ Dimer-Adatom-Stacking-Fault (DAS)-reconstructed (111) surface of silicon is one of the most complicated surface systems (Fig. 1). The intricacies of this reconstructed surface structure render the study of its surface dynamics challenging as well as rewarding. Experimentally, an EELS measurement Da87  of the clean $`7\times 7`$ surface revealed dynamical features at 25-33 meV and 71 meV. A high resolution HAS study of the low energy surface modes showed Rayleigh waves (RW) peaked at 8, 10, and 15 meV La98 . Theoretical studies of surface vibrations of the reconstructed clean $`7\times 7`$ silicon (111) surface include an ab initio cluster calculation Da87 , a total energy and vibrational spectrum-study based on empirical potentials Li88 , a molecular dynamics (MD) simulation study of the phonon spectrum based on Car-Parrinnello fictitious Lagrangian and a tight-binding Hamiltonian Ki95 , and a first-principles finite-temperature MD study of the surface dynamics St95 ; St96 . The two MD studies had yielded similar surface dynamical features. Since both MD studies employed a supercell of 400 atoms and the convergence criterion for the force of 0.1 eV/Å, the comparison suggests that tight-binding approaches will give comparable result as that obtained by first-principles approaches, provided that both approaches used the same computational framework. However, with the imposed inversion symmetry and atoms in the central two layers in the supercell held at bulk positions in the two MD studies, only 151 atoms (adatoms plus atoms in the first three layers) were allowed to relax in the structural determination or to move in the dynamical simulations. Thus the extent of the relaxation in the surface structural determination and the coupling of the surface to the bulk vibrational modes are not adequately taken into account in the two MD studies. Furthermore, the MD trajectory was followed for only 1.2 ps in the first-principles MD study, suggesting that the first-principles MD study may not have provided sufficiently accurate description of the low energy modes. Because of the complexity of the surface structure, experimentalists need to have guiding input to efficiently measure surface dynamical features. Therefore, a reliable and comprehensive theoretical investigation of the surface dynamics of the clean $`7\times 7`$-reconstruted (111) surface is urgently needed.

We carried out a total-energy-vibrational-spectrum study of the clean silicon reconstructed (111) surface. Our convergence test showed that the energy optimization using a supercell composed of 10 layers and adatoms (494 atoms), with atoms in the bottom 4 layers held at bulk positions to mimic the effect of the bulk and periodic boundary conditions imposed horizontally, is sufficient to ensure a reliable surface structural determination. The energy optimization was implemented with a non-orthogonal tight-binding (NOTB) HamiltonianMe97  because of the system size (494 atoms). This Hamiltonian has been found to yield the energetics and the reconstruction of the clean Si (001) surface in excellent agreement with those obtained by DFT calculations as well as experimental resultsLi00 . Our simulation has yielded structural parameters for the Si (111) surface comparable to previous studiesKi95 ; St92  with only one notable exception. We found the distance between the adatom and the atom in layer 2 directly underneath it to be about 2.97Å, compared with a result of 2.60Å of Ref.  and a result of 2.40Å of Ref. . This substantial difference can be traced to the fact that previous studies  relaxed only adatoms and atoms in the first three layers (151 atoms) at a convergence criterion for the force of 0.1 eV/Å, with only one layer of atoms held at bulk positions while our optimization relaxes 298 atoms with the force criterion at 10<sup>-3</sup> eV/Å and with four layers held at bulk positions. Thus our optimization is expected to provide a more accurate determination of the surface structure. We merged the optimized structure of the 10-layer slab with the bulk structure to form the semi-infinite $`7\times 7`$ reconstructed DAS (111) surface system. The vibrational spectrum of this semi-infinite system was calculated layer by layer using the method of real space Green’s functionWu95 . In the calculation, the force constant matrix for the adatoms and the atoms in the first 6 layers (298 atoms) in the semi-infinite system and the coupling force constant matrix between atoms in the first 6 layers and those in layers 7-10 (interfacial layers) were computed numerically as second partial derivatives of the total energy of the original 10-layer slab about its equilibrium configuration. The force constant matrix for the atoms in layers 7-10 and those between these atoms and atoms in subsequent bulk layers were again computed numerically as second partial derivatives of the total energy of a slab composed of layers 3-10 in the original slab plus 4 bulk layers. The bulk force constant matrix was used as that for the remaining atoms in this semi-infinite system. In this way, the dynamical matrix of the semi-infinite system was constructed continuously and smoothly from the surface region, to the interfacial region, and to the bulk region. The vibrations of all the atoms in the system were therefore treated on equal footing. Our study has identified and elucidated, for the first time, all the important dynamical features associated with the dimers, adatoms, stacking faults, and rest atoms.

Fig. 2 presents SDOS for certain adatoms and certain atoms in the first 4 layers to highlight important surface features and illustrate the underlying physics. They are all presented in the same scale, i.e., $`\rho _{i\alpha }(\omega )=(2\omega /\pi )Im_{\epsilon 0}G_{i\alpha ,i\alpha }(\omega ^2+i\epsilon )`$ with $`G(\omega ^2)=(\omega ^2\varphi )^1`$, $`\varphi `$ being the mass-normalized force constant matrix, $`i`$ denoting the site, and $`\alpha `$ the direction, so as to allow meaningful comparisons between them. We have calculated SDOSs for atoms in the 21st layer and found them to be very close to each other and in good agreement with the bulk DOS fitted to experiment We77 . Hence, in Fig. 2(a) we use the SDOS for an atom in the 21st layer as the SDOS for a typical atom in the bulk. The SDOSs in Fig. 2 are all decomposed along $`x(\overline{1}\overline{1}2)`$, $`y(1\overline{1}0)`$, and $`z(111)`$ directions. Fig. 2(a) shows that the SDOSs along $`x`$, $`y`$, and $`z`$ directions for a bulk atom are practically all the same. Hence one can use the extent of deviation of SDOSs along $`x`$, $`y`$, and $`z`$ directions from each other as a gauge of the surface effect. To provide useful guidelines for the experimental investigation, we present in Figs. 2(b) through 2(i) the SDOSs for atoms (see Fig. 1) that play significant roles to surface modes. Figs. 2(b) and 2(c) give the SDOSs for the corner adatom (COA or atom L0A) and central adatom (CEA or atom L0B) in the unfaulted subcell respectively. The SDOSs for atom L0A are similar to those of atom L0B, with only a slight shift in peak frequencies. They both have prominent surface features at about 14 (in-plane), 17 (in-plane), 62 (out-of-plane), 63 (out-of-plane) and 75 (out-of-plane) meV. We did not show SDOSs for adatoms in the faulted subcell because they are very similar to the corresponding SDOSs in the unfaulted subcell. In fact, the features of the SDOSs for the corresponding sites in both halves are in general similar, with a few notable exceptions to be discussed later. A comparison of Figs. 2(c), 2(f) (SDOSs for atom L2B in layer 2 directly underneath atom L0B), and 2(h) (SDOSs for atom L3B directly underneath atom L2B) immediately identifies the split-off mode at 75 meV as originating from the $`z`$-vibration of the compressed pair of atoms L2B and L3B and attributable to the 71 meV-mode observed by EELS, as also noted by previous theoretical studies Ki95 ; St95 ; St96 . An examination of Figs. 2(b) and 2(d) (SDOSs for the rest atom L1C with a dangling bond in the unfaulted subcell of layer 1) reveals that the mode at 63 meV actually originates from the rest atom L1C while Figs. 2(b) and 2(g) (SDOSs for atom L2E in a dimer in layer 2 adjacent to the COA L0A) show that the mode at 62 meV originates from the dimer. The rest atom mode at 63 meV can be traced from its source to the COA L0A by examining Figs. 2(d), 3(b) (SDOSs of the atom L2G in layer 2 that is a nearest neighbor of the rest atom L1C), 3(a) (SDOSs of the atom L1F in layer 1 that is a nearest neighbor of both the COA L0A and the atom L2G), and 1(b). The pronounced mode at 63 meV with almost equal $`x`$\- and $`y`$-polarization originating from the rest atom L1C is seen to induce a vibration at its neighboring atom L2G polarized almost in the $`xz`$ plane, which in turn excites its neighboring atom L1F to vibrate almost in $`y`$-direction, yielding eventually an almost $`z`$-polarized vibration for the COA L0A. From Figs. 2(g), 3(a), and 2(b), the pronounced mode at 62 meV polarized along the dimer axis ($`y`$-direction) of the dimer atom L2E is seen to have excited its neighboring atom L1F to vibrate with its polarization dominated by the $`x`$-vibration, which in turn induces the $`z`$-polarized vibration at the COA L0A. Our calculation has also shown that these three modes are localized surface modes, with the 62 (dimer) and 75 meV (compressed atom pair) modes barely discernable at layer 4. However, the rest atom mode at 63 meV extends beyond layer 4 and exhibits different polarization in layer 4 in the faulted subcell (Fig. 2(i)) as compared to the unfaulted subcell. This can be understood as a feature associated with the stacking fault. For the rest atom L1D in the faulted subcell of layer 1, because of the stacking fault, there is an atom directly underneath it in the faulted subcell of layer 4 whose SDOSs are shown in Fig. 2(i). The 63 meV mode still shows up as an in-plane mode polarized almost equally along the $`x`$ and $`y`$-axis, the same as the mode at its rest-atom origin (see Fig. 2(e)). However, there is no atom in the unfaulted subcell of layer 4 directly underneath the rest atom in the unfaulted half of layer 1. This mode shows up in the unfaulted half of layer 4 as a mode polarized close to the $`x`$-axis. The effect of the stacking fault, however, shows up most prominently in the low frequency strongly $`z`$-polarized mode at about 16 meV in the SDOSs of the rest atom L1D in the faulted half of layer 1. A comparison of Figs. 2(d) and 2(e) reveals that such a mode is not present in the SDOSs of the rest atom L1C in the unfaulted half. To shed light on this unusual contrast, we decompose the SDOS in terms of off-diagonal density matrix elements ($`\rho _{i\alpha ,j\beta }`$) as follows Li01 :

$$\rho _{i\alpha }(\omega )=\frac{_{\beta (\alpha )}\varphi _{i\alpha ,i\beta }\rho _{i\beta ,i\alpha }+_{\beta ,j(i)}\varphi _{i\alpha ,j\beta }\rho _{j\beta ,i\alpha }}{\omega ^2\varphi _{i\alpha ,i\alpha }}.$$
(1)
Because of the presence of $`\varphi _{i\alpha ,j\beta }`$, only directly coupled vibrations will contribute to the SDOS. Note that $`\rho _{i\alpha ,j\beta }`$ also determines the phase relation between the vibration at $`i\alpha `$ and that at $`j\beta `$. Fig. 4 gives the SDOSs along the $`z`$-axis in the vicinity of 16 meV (solid line), decomposed according to $`R_{ij}3.15\AA `$ (long line segment), $`3.15\AA <R_{ij}4.15\AA `$ (short line segment), and $`4.15\AA <R_{ij}`$ (dot-dashed line), for rest atoms in the unfaulted (Fig. 4(a)) and faulted (Fig. 4(b)) halves respectively. Figs. 4(b) and 2(i) show that the most significant contribution to the 16 meV-mode originated from the rest atom in the faulted half comes from its coupled vibration with the atom in layer 4 directly underneath it (distance = 2.29Å), with some participation of its three nearest neighbors and the three nearest neighbors in layer 3 of the atom in layer 4. Our results also reveals that this mode is localized in the vicinity of the configuration composed of the eight atoms mentioned above and anchored by the out-of-phase $`z`$ vibrations of the rest atom L1D and the layer 4-atom L4D underneath it. Since there is no such configuration for the rest atom in the unfaulted half, no such mode can exist for the rest atom in the unfaulted half (see Fig. 4(a)). Hence the presence and the absence of the 16 meV-mode at the rest atom L1D in the faulted and the rest atom L1C in the unfaulted half respectively is a direct consequence of the stacking fault. Figs. 2(f) and 2(c) show another interesting surface mode at about 57 meV. This mode originates as a $`x`$-polarized vibration at the atom L2B in layer 2 directly underneath the CEA L0B in the unfaulted half and emerges as a $`z`$-polarized vibration of CEA L0B, with no appreciable penetration to atoms in layer 3 and beyond. The decomposition of the SDOS at atom L2B where the mode originates using Eq. (1) indicates that atoms at a distance $`R>3.15\AA `$ from atom L2B give insignificant contributions to the SDOS at atom L2B. We therefore present in Fig. 5(a) the decomposition of the SDOS along the $`x`$-axis of atom L2B in the vicinity of 57 meV in terms of off-diagonal density matrix elements only for those atoms at a distance less than 3.15Å from atom L2B. It can be seen that almost the entire contribution to $`\rho _{L2Bx}`$ comes from its three neighboring atoms in layer 1, with no meaningful contribution from either atom L0B (directly above) or atom L3B (directly underneath). The absence of the contribution from atom L0B to the SDOS of atom L2B and the presence of the 57 meV- mode as a $`z`$-polarized vibration in the SDOS of L0B (Fig. 2(c)) is a clear indication that the $`z`$-vibration of atom L0B is not directly coupled to the $`x`$-vibration of atom L2B. Thus the 57 meV-mode is a mode localized in the first two layers. It originates from the coupled vibrations of the $`x`$-polarized vibration of atom L2B with its three neighbors in layer 1 and, in turn, induces the $`z`$-polarized vibration of CEA L0B. One interesting feature of the method of analysis based on Eq. (1) may be exemplified by the examination of the $`z`$-polarized mode at 45 meV of Fig. 2(h). This is apparently a mode restricted to the immediate neighborhood of atom L3B as there is no evidence of such a mode in the SDOSs of its neighboring atoms, including atom L2B, its nearest neighbor. We present in Fig. 5(b) the decomposition of $`\rho _{L3Bz}`$ in terms of $`_{\beta ,j(L3B)}\varphi _{L3Bz,j\beta }\rho _{j\beta ,L3Bz}/(\omega ^2\varphi _{L3Bz,L3Bz})`$ for atoms at distances $`R<3.15\AA `$ from atom L3B. Figs. 4 and 5 show that there is a common feature among all the decomposition curves, namely, the simultaneous appearance of pronounced peaks and valleys at the frequency corresponding to $`\omega =\sqrt{\varphi _{i\alpha ,i\alpha }}`$, the natural frequency of vibration for the atom $`i`$ along the $`\alpha `$ direction (see Eq. (1)). The combined peaks and valleys usually lead to a complete cancellation at the natural vibrational frequency and hence no contribution to the SDOS as the atom is not isolated. A comparison of Figs. 5(a) and 5(b) however reveals that the natural vibrational frequency of atom L3B along the $`z`$-direction is about 47 meV while that of atom L2B along the x-direction is about 45 meV. Thus the $`z`$-polarized mode at 45 meV of atom L3B is the consequence of the perturbation of the natural vibration of atom L3B along $`z`$-direction by the natural vibration of atom L2B along the $`x`$-direction (see Fig. 5(b)). Finally, by comparing SDOSs of corresponding atoms layer-by-layer, we have identified RWs at 15, 9, and 7 meV, in good agreement with the experimental result of Ref. .

Our comprehensive study of the dynamics of the $`7\times 7`$ DAS-reconstructed (111) has succinctly explained all the important dynamical features associated with the structural features of the reconstructed surface. These findings should provide useful guidelines for experimental studies of the dynamics of the $`7\times 7`$ reconstructed surface.

###### Acknowledgements.

This work was supported by the NSF (DMR-011284) and the DOE (DE-FG02-00ER45832) grants.
FIGURE CAPTIONS Fig.1 A schematic presentation of the $`7\times 7`$ DAS-reconstructed Si (111) surface. The sites used in Figs 2 through 5 are labeled by L0A, L1C, $`\mathrm{}`$, a combination of the layer marker (L0, adatom; L1, layer 1, $`\mathrm{}`$) and the position marker (A, B, $`\mathrm{}`$). Fig. 2: SDOSs of a typical bulk atom (obtained from an atom in the 21st layer) (a) and of atoms at sites L0A (b), L0B (c), L1C (d), L1D (e), L2B (f), L2E (g), L3B (h), and L4D (i) respectively. Fig. 3: SDOSs of atoms at sites L1F (a) and L2G (b) respectively. Fig. 4: The decomposition of SDOSs $`\rho _{L1Cz}`$ (a) and $`\rho _{L1Dz}`$ (b) (solid line) using Eq. (1) for $`R_{ij}<3.15\AA `$ (long line segment), $`3.15\AA <R_{ij}<4.15\AA `$ (short line segment), and $`R_{ij}>4.15\AA `$ (dot-dashed line) respectively. Fig. 5: The decomposition of SDOSs $`\rho _{L2Bx}`$ (a) and $`\rho _{L3Bz}`$ (b) (solid line) using Eq. (1) for $`R_{ij}<3.15\AA `$ respectively. In (a), long line segment: coupling between L2B and the atom in layer 1 to the left of L2B; short line segment: the average of the coupling to the two atoms in later 1 to the right of L2B; dot-dashed line: coupling to atom L3B. In (b), long line segment: coupling between atom L3B to atom L2B; short line segment: the average of the coupling to the three atoms in layer 4.