format: text

text:
# Simulations of the kinetic friction due to adsorbed surface layers

## 1 Introduction

One of the major challenges for researchers in the field of nanotribology is to uncover the molecular origins of macroscopic behavior. For example, why do many systems obey Amontons’ laws that friction is proportional to normal load $`L`$ and independent of the apparent geometric area of the surfaces $`A_{\mathrm{app}}`$? Unfortunately many simple molecular scale models give results that are qualitatively inconsistent with macroscopic measurements . Indeed they suggest that static friction should almost never be observed ! A possible origin for such discrepancies  is that the models consider contacts between two bare surfaces. Real surfaces almost always have an intervening layer of material called ”third bodies” . This layer may be dust, grit, wear debris, or the few angstroms of hydrocarbons and water that adsorb rapidly on surfaces that are exposed to air.

In a recent paper He et al.  explored the effect of adsorbed layers on the static friction $`F_s`$. They showed that such layers naturally lead to a non-vanishing $`F_s`$ that exhibits much of the behavior seen in macroscopic experiments. For example the calculated value of $`F_s`$ is insensitive to many parameters that are not controlled in experiments, including crystallographic alignment, the thickness of the adsorbed layer, and the precise size of the molecules. More importantly, the simulation results provided a molecular basis for Bowden and Tabor’s phenomenological explanation of Amontons’ laws .

Bowden and Tabor had noted that $`A_{\mathrm{app}}`$ is usually much larger than the area of intimate molecular contact $`A_{\mathrm{real}}`$. They pointed out that both Amontons’ laws and many exceptions to them could be explained if the local shear stress $`\tau `$ in the contacts increased linearly with the local pressure $`P`$:

$$\tau =\tau _0+\alpha P.$$
(1)
Summing over $`A_{\mathrm{real}}`$ and dividing by load gives a friction coefficient

$$\mu F/L=\alpha +\tau _0/P.$$
(2)
This is independent of load if $`P`$ is constant or $`\tau _0/P\alpha `$. Later work showed that both elastic and plastic models of surfaces with many contacts yield a constant $`P`$ at high loads. At low loads, $`P`$ decreases, explaining why many systems show an increase in $`\mu `$ in this limit .

He et al.’s work  showed that adsorbed layers lead to a static friction that obeys Eq. 1. However, most experiments measure the kinetic friction. Although experimental values for the two quantities are often correlated, it is not obvious that they should be. Static friction is the threshold force that must be overcome before one solid slides over another. It arises because the two bodies have locked together into a local energy minimum, and $`F_\mathrm{s}`$ is the force needed to break them free from this minimum. Kinetic friction $`F_k(v)`$ is the force needed to maintain sliding at a given velocity $`v`$. As such it is proportional to the work that must be done due to energy dissipation during sliding. Thus kinetic and static friction involve different phenomena and they have very different values in many simple models .

In this paper we consider the effect of adsorbed layers on kinetic friction. We first show that the linear relation in Eq. 1 is obeyed for each velocity. The intercept $`\tau _0`$ is relatively independent of velocity below about 1 m/s. The slope $`\alpha `$ shows a logarithmic dependence on velocity over more than two decades in velocity. This same functional dependence is observed in experiments, and used in rate-state models of friction . We next examine the effect of interaction potentials, temperature, and crystallographic alignment on kinetic friction. The variations with these parameters are consistent with a simple hard sphere model for lateral forces. We conclude by considering the strong correlation between calculated values of $`F_k`$ and $`F_s`$. A qualitative explanation for this correlation is suggested based on results for the Tomlinson model , and this analogy is confirmed by detailed analysis of the dynamics of individual molecules.

The following section describes the model we use for the surfaces and the intervening molecular layer. The next section presents results, and the final section provides a summary and conclusions.

## 2 Computational Methodology

Since our goal is to establish general mechanisms, we use relatively simple interaction potentials in our simulations. These allow us to quickly span a wide range of sliding velocities, system sizes, geometries and interactions. They also allow treatment of longer time and length scales than more detailed potentials.

We find similar behavior for both crystalline and disordered walls . For the results described below, each solid contains discrete atoms forming the (111) surface of an fcc crystal. The coordinate system is chosen so that the surfaces lie in the $`xy`$ plane, and periodic boundary conditions are applied in these directions to minimize edge effects. Note that the periodic boundary conditions also prevent us from considering perfectly incommensurate systems. However, we find that surfaces exhibit incommensurate behavior as soon as the common period becomes longer than a few lattice constants .

The two opposing surfaces have approximately the same nearest-neighbor spacing $`d_{nn}`$, and the effect of commensurability is explored by rotating the top surface by an angle $`\theta `$ about the $`z`$-axis. A discrete set of angles is studied that allow both surfaces to retain their hexagonal symmetry. Except where noted, the bottom wall has $`d_{nn}=1.2\sigma `$. The value of $`d_{nn}`$ for the top wall is adjusted slightly ($`2.1`$%) in order to satisfy the periodic boundary conditions. No correlation between the size of this adjustment and the calculated friction was seen at the level of the statistical fluctuations.

Atoms in the bottom solid are fixed in place, and atoms in the top solid translate together as a rigid body of mass $`M`$. Previous work shows that including elastic deformations of the solids produces a slightly lower static friction , and we find a similar reduction in the kinetic friction. A constant downward pressure is applied to the top solid, and its height varies in response to this external pressure and molecular interactions. A constant lateral force or velocity is applied along the nominal sliding direction, usually $`\widehat{x}`$. The wall is free to move in the other lateral direction, i.e. the external force is zero in this direction. The kinetic friction is determined by imposing a constant lateral velocity $`v`$, and calculating the average opposing force. This average is performed over a displacement of at least 25$`d_{nn}`$ in order to reduce statistical fluctuations. The static friction is obtained by applying a constant lateral force, and determining the value needed to initiate sliding .

The layer of molecules between the solids is described with a bead-spring model. Previous studies have shown that this model yields realistic dynamics for polymer melts , and shown how to map between it and detailed chemical models of polymers . These detailed models take at least an order of magnitude more computer time.

Each molecule contains $`n`$ spherical monomers of mass $`m`$. We considered simple spherical molecules $`n=1`$, and short chains that model airborn hydrocarbons, $`n=3`$ and 6. All monomers interact through a truncated Lennard-Jones potential :

$$V_{LJ}(r)=4ϵ\left[(\sigma /r)^{12}(\sigma /r)^6\right]\mathrm{for}r<r^c,$$
(3)
where $`r`$ is the distance between molecules and the potential is zero for $`r>r^c`$. Adjacent monomers along a chain are coupled by an additional potential that prevents chains from crossing or breaking:

$$V_{CH}(r)=(1/2)kR_o^2\mathrm{ln}[1(r/R_o)^2],$$
(4)
where $`R_o=1.5\sigma `$ and $`k=30ϵ/\sigma ^2`$ . The parameters $`ϵ`$, $`\sigma `$, and $`t_{LJ}(\sigma ^2/mϵ)^{1/2}`$ are characteristic energy, length, and time scales, respectively. These characteristic scales are used to normalize other quantities. Values that would be representative of hydrocarbons are $`ϵ30`$meV, $`\sigma 0.5`$nm and $`t_{LJ}3`$ps . The unit of pressure, $`ϵ\sigma ^3`$, corresponds to about 40MPa.

Wall atoms interact with monomers through a Lennard-Jones potential with a different cutoff, $`r_{wf}^c`$, and different energy and length scales, $`ϵ_{wf}`$ and $`\sigma _{wf}`$, respectively. This allows us to increase or decrease the amount of adhesion, and the effective surface corrugation of the walls. In most cases we considered purely repulsive interactions: $`r^c=2^{1/6}\sigma `$ and $`r_{wf}^c=2^{1/6}\sigma _{wf}`$. In this limit, atoms on opposing walls are far enough apart that they do not interact. As shown below, the main effect of increasing the interaction range is to add an extra adhesive pressure between the walls.

All of our runs are done in a constant temperature ensemble. We have performed extensive tests on the effects of thermostats . In the appropriate limits, thermostats prevent a gradual rise in temperature due to frictional dissipation without having a direct influence on the friction force. To ensure that one is in these limits, the thermostat must be applied normal to the sliding direction and be sufficiently weak, and the velocity must be low ($`0.1\sigma /t_{LJ}`$). At high velocities, the structure and other properties of the film may be affected . In the following simulations, a Langevin thermostat with damping rate $`0.4t_{LJ}^1`$ is coupled to the equations of motion in both directions normal to the sliding direction .

The equations of motion are integrated using a fifth-order Gear predictor-corrector algorithm with a time step $`dt=0.005\sigma /t_{LJ}`$. Unless otherwise noted, the area of a periodic unit cell is $`𝒜=720\sigma ^2`$ and the mass of the wall is $`M/m=0.8𝒜/\sigma ^2`$. Previous studies show that finite-size effects are very small at this system size . Most simulations are done at temperature $`T=0.7k_B/ϵ`$, where $`k_B`$ is Boltzmann’s constant.

## 3 Results

### 3.1 Variation with pressure and velocity

Figure 1 shows the pressure dependence of the kinetic friction at the indicated velocities for an incommensurate system, $`\theta =30^{}`$. At each velocity the shear stress $`\tau _k`$ rises linearly with pressure over the entire range of pressures. Using realistic values of $`ϵ`$ and $`\sigma `$, the highest pressures correspond to about 1.5GPa. The linear relation continues up to $`P=100ϵ\sigma ^3`$ in the cases we have tested. However, the larger forces at these pressures require a reduction in the integration time step, making the calculations much slower.

Values of $`\alpha `$ and $`\tau _0`$ were obtained from linear fits to data like that shown in Fig. 1. In each case the correlation coefficient was greater than .999. The results are plotted against the logarithm of velocity in Figure 2. The value of $`\alpha `$ increases as the logarithm of velocity over the entire range of Fig. 2 – a factor of 200 in $`v`$. In contrast, the value of $`\tau _0`$ becomes independent of velocity for $`v<0.01\sigma /t_{LJ}1`$m/s. It is negative because these simulations used a purely repulsive potential $`r^c/\sigma =2^{1/6}`$. We show later that including adhesive interactions leads to positive values of $`\tau _0`$.

A logarithmic variation of friction with velocity is also seen in many experimental systems . Measured friction coefficients fit a “rate-state” equation:

$$\mu =\alpha +\tau _0/P=\mu _0+A\mathrm{ln}(v/v^0)+B\mathrm{ln}(\mathrm{\Theta }/\mathrm{\Theta }^0),$$
(5)
where $`\mathrm{\Theta }`$ is a variable that describes the state of the system, and $`\mu _0`$ is the friction at the reference velocity $`v^0`$ and state $`\mathrm{\Theta }^0`$. In most cases $`\mathrm{\Theta }`$ is treated as a phenomenological parameter whose evolution is described by:

$$d\mathrm{\Theta }/dt=1\mathrm{\Theta }v/D_c$$
(6)
where $`D_c`$ has units of length. Dieterich and Kilgore have shown that changes in $`\mathrm{\Theta }`$ are proportional changes in the true area of contact between rough surfaces and that $`D_c`$ corresponds to the typical diameter of contacts .

Our simulations give the dependence of friction on $`v`$ at fixed contact area  and thus the slope of Fig. 2(a) gives $`A=0.00110\pm 0.00005`$. This is smaller than typical experimental values for glass or rocks where $`A`$ is .005 to .015. However, one might expect $`A`$ to increase with $`\alpha `$. While our values of $`\alpha `$ are comparable to values for flat crystalline surfaces like mica or MoS2, they are much smaller than values for rocks. The ratio $`A/\alpha `$ is about $`.05`$ for our system, and varies from .008 to .025 for rocks. We do not know of any measurements of $`A`$ for mica or MoS2, but it would be interesting to see how results for these systems compare to our simulations.

### 3.2 Variation of $`\mu `$ with interaction potential, $`T`$ and geometry

A variety of interaction potentials and geometries have been considered to determine which factors influence $`\alpha `$ and $`\tau _0`$ and which leave them unchanged. As shown in Fig. 3, doubling the interaction between wall atoms and fluid monomers has almost no effect on the friction. In contrast, the friction increases rapidly with the ratio of $`d_{nn}`$ to $`\sigma _{wf}`$. For $`d_{nn}=1.2\sigma `$ we find $`\alpha `$ rises from 0.012 to 0.034 as $`\sigma _{wf}`$ drops from 1.2 to 0.9. The value of $`\tau _0`$ increases in magnitude from -0.07 to -0.16. Although not shown, increasing $`d_{nn}`$ at fixed $`\sigma _{wf}`$ also increases $`\alpha `$.

These trends in $`\alpha `$ can be understood from a simple geometrical picture. At the pressures of interest, the interactions are dominated by the repulsive $`r^{12}`$ term in the Lennard-Jones potential. Monomers and wall atoms can be thought of as hard spheres with an effective diameter $`D`$ that is determined by the separation $`r`$ at which the Lennard-Jones force balances the force from the applied pressure. This picture gives a diameter $`D\sigma _{wf}(ϵ_{wf}/\sigma _{wf}P)^{1/13}`$, that varies rapidly with $`\sigma _{wf}`$ but very little with $`ϵ_{wf}`$ or $`P`$. As illustrated in Fig. 4, monomers can penetrate deep between surface atoms when $`D/d_{nn}`$ is small. This increases the friction because monomers must be lifted up the ramp formed by the surface of closest approach in order to translate. The required lateral force is the normal load times the slope of the surface which increases as $`D/d_{nn}`$ decreases. Thus the hard-sphere model explains the linear dependence of $`\tau `$ on pressure as well as the insensitivity to $`ϵ_{wf}`$ and the variation with $`\sigma _{wf}`$.

The simulations just described were for strictly repulsive potentials $`(r_{wf}^c=2^{1/6}\sigma _{wf})`$, and one may wonder whether the hard-sphere picture is more generally applicable. We have done simulations with $`r_{wf}^c=2.2\sigma _{wf}`$ to test this. When direct wall/wall interactions are neglected, we find that the kinetic friction is fit by $`\alpha =0.0211\pm 0.0003`$ and $`\tau _0=+0.007\pm .008`$ over a pressure range from 0 to 36$`ϵ\sigma ^3`$ for the parameters of Fig. 1. This value of $`\alpha `$ is only 8% smaller than for a purely repulsive potential. In contrast, $`\tau _0`$ increases substantially because the attraction between monomers and wall atoms increases the effective pressure holding them together. Adding direct wall/wall interactions leads to a slightly larger adhesive pressure, and a slightly larger value of $`\tau _0`$.

The effect of temperature on $`\alpha `$ is also fairly weak. So far we have shown results for $`T=0.7ϵ/k_B`$, which is near the triple point for spherical Lennard-Jones molecules, and roughly 30% above the glass transition temperature for chain molecules. This temperature was chosen because many airborn hydrocarbons are above their melting point at room temperature. Figure 5 shows the variation of $`\alpha `$ and $`\tau _0`$ with temperature at $`v=0.005\sigma /t_{LJ}`$ for the system of Figs. 1 and 2. The value of $`\alpha `$ changes only weakly with temperature, falling about 20% from $`k_BT/ϵ=0.3`$ to 1.2. The high end of this temperature range is more than twice the glass transition temperature for chain molecules and close to the critical temperature for spherical molecules.

The weak variation of $`\alpha `$ with $`T`$ is consistent with our geometrical picture, since the hard-sphere diameter is insensitive to temperature. In contrast, the value of $`\tau _0`$ is roughly proportional to temperature over the studied range. Equation 1 implies that the friction vanishes for $`P<\tau _0/\alpha `$. Of course this can not happen, but there is a qualitative change in the type of friction as the pressure drops below $`\tau _0/\alpha `$. At high pressures we observe static friction and a slightly smaller kinetic friction that rises logarithmically with velocity. For $`P<\tau _0/\alpha `$ the static friction vanishes, and the kinetic friction rises from zero as a power of the velocity.

This change in frictional behavior reflects a change in the state of the film that has been studied previously in multi-layer systems . At high pressures the molecules lock into a glassy state, while at low pressures they diffuse freely and act like a viscous lubricant. Fig. 5 implies that the pressure required to lock molecules into a glassy state rises linearly with $`T`$. This result is specific to our choice of purely repulsive potentials ($`r^c=2^{1/6}\sigma `$). Adding an attractive tail to the potential shifts the entire temperature dependent curve upwards, and gives glassy behavior at lower (even negative) pressures. Note that earlier results on multi-layer films indicated that $`\tau `$ was independent of velocity in the glassy state . However, these results covered a relatively small range of $`v`$ and are not inconsistent with the logarithmic variation seen in Fig. 2.

### 3.3 Relation to static friction

The changes in kinetic friction with temperature, interaction potential, geometry and other factors are very similar to what we have observed for static friction . Indeed the ratio between the value of $`\alpha _k(v)`$ at a given velocity and the value of $`\alpha _s`$ for static friction is relatively constant. For example, over the range of temperatures and potential parameters shown in Figures 3 and 5, $`\alpha _k(v=0.005\sigma /t_{LJ})`$ is between 75 and 85% of $`\alpha _s`$. This variation is comparable to our statistical uncertainties.

Figure 6 illustrates the correlation between static and kinetic friction as $`\theta `$ is varied. Only the range $`\theta =0`$ to 30 is shown, since the values at other $`\theta `$ are related by symmetry. The simulation data in previous figures was for $`\theta =30^{}`$ and the commensurate case corresponds to $`\theta =0^{}`$. We confirmed that the friction is linear in $`P`$ for several $`\theta `$, and show results for a single, relatively high pressure ($`P=24ϵ\sigma ^3`$) to minimize the number of calculations, and emphasize the contribution from $`\alpha `$. Data for a monolayer film is shown because it gives the largest variation in friction with $`\theta `$ and thus the best test of correlations between $`F_k`$ and $`F_s`$ . As $`\theta `$ changes from 2 to 30, $`F_k`$ varies by almost a factor of two, yet the ratio $`F_k/F_s`$ stays constant within our statistical uncertainties. For $`v=0.005\sigma /t_{LJ}`$ we find $`F_k/F_s=0.88\pm 0.10`$, while $`F_k/F_s=1.05\pm 0.10`$ for $`v=0.05\sigma /t_{LJ}`$.

The correlation between $`F_k`$ and $`F_s`$ is lost when the walls become commensurate. As $`\theta `$ decreases from 2 to $`0^{}`$ there is a sharp increase in $`F_s`$ by about a factor of seven, and a drop in $`F_k`$ by about a factor of two. In addition, there is a qualitative change in the velocity dependence of the kinetic friction: $`F_k`$ decreases with velocity in the commensurate case, but increases with $`v`$ for incommensurate walls.

The observed correlation between static and kinetic friction is somewhat surprising given that they reflect different physical effects. As noted in the introduction, the static friction is related to the force needed to pull the system out of a local potential energy minimum, while the kinetic friction is related to the energy dissipated during sliding. Nevertheless, the above comparisons of $`F_k`$ and $`F_s`$ indicate that the two arise from similar processes when the walls are incommensurate.

### 3.4 Comparison to the Tomlinson model

Some insight into the relation between static and kinetic friction is provided by a simple one-dimensional ball and spring model called the Tomlinson model . This model considers two bare surfaces in contact. The bottom surface is modeled by a sinusoidal potential with amplitude $`V_0`$ and period $`d_{nn}`$. Atoms in the top surface are treated as independent, overdamped oscillators that are coupled to the center of mass $`X_{CM}`$ by a spring of stiffness $`k`$.

The behavior of the Tomlinson model depends on the dimensionless parameter $`\lambda 4\pi ^2V_0/kd_{nn}^2`$ that measures the relative stiffness of the potential and spring. When $`\lambda `$ is small, the spring is so stiff that each atom has only one metastable state at any $`X_{CM}`$. In this limit atoms move smoothly over the periodic potential from the bottom wall. This nearly reversible motion leads to little dissipation, and $`F_k`$ goes to zero linearly with the center of mass velocity . When $`\lambda `$ is large, each atom can get locked in many metastable states at any center of mass position. Each atom starts in one such state. As $`X_{CM}`$ increases, that state eventually becomes unstable and the atom pops rapidly to the next state. During these pops the system is far from equilibrium and a substantial amount of energy is dissipated. The sequence of pops and the amount of energy dissipated during a given wall displacement $`\mathrm{\Delta }X_{CM}`$ is independent of wall velocity in the limit $`V_{CM}0`$. Since the dissipated energy equals $`F_k(v)\mathrm{\Delta }X_{CM}`$, the kinetic friction approaches a constant value, $`F_k(0)`$, in this limit. One can show that $`F_k(0)`$ approaches the maximum force exerted by the periodic potential $`F_{max}2\pi V_0/d_{nn}`$ times the number of atoms $`N`$ as $`\lambda \mathrm{}`$ .

The static friction in the Tomlinson model depends on whether the walls are commensurate or incommensurate . In the commensurate case all atoms feel the same periodic potential at each $`X_{CM}`$. The maximum force that must be overcome as they advance by a period is $`F_s=NF_{max}`$. The ratio $`F_k/F_s`$ vanishes for $`\lambda <1`$ and approaches unity as $`\lambda \mathrm{}`$. In the incommensurate case $`F_s=0`$ for $`\lambda <1`$. For $`\lambda >1`$, the static friction is equal to the low velocity kinetic friction. This link between $`F_s`$ and $`F_k`$ can be seen by calculating $`F_k(0)`$ from instantaneous forces rather than the dissipated energy. At very low velocities almost all of the atoms are stuck in a metastable minimum at any instant in time. The force resisting motion, $`F_k(0)`$, comes from the force exerted by the periodic potential on these atoms. One can show that summing this force or calculating the energy dissipated during atomic pops gives the same answer for $`F_k`$ . For incommensurate walls the distribution of metastable states and thus the lateral force is independent of position. Hence the static friction equals the force $`F_k(0)`$ needed to advance the wall at a very low velocity.

In most cases two bare surfaces will be incommensurate and the interactions between them will be weaker than their internal interactions. This corresponds to the case $`\lambda <1`$ and one expects the static friction and low velocity kinetic friction to vanish . In contrast, the layer of adsorbed molecules included in our calculations acts like a very compliant solid with $`\lambda >1`$. As expected from the Tomlinson model, this produces a finite static friction and a comparable kinetic friction. The lateral force is predominantly determined by the hard sphere interactions which can be expressed as the normal load times the local slope of the surface of closest approach (Figure 4). Since the surface is nearly independent of load and other parameters, the lateral force is proportional to load and Amontons’ laws are obeyed.

We have confirmed the connection between our results and the Tomlinson model by analyzing the motion of individual monomers at low velocities. The number of particles that move by more than a distance $`c`$ over a time interval of $`2t_{LJ}`$ is calculated as a function of $`c`$. This time interval is chosen because it is longer than the velocity-velocity autocorrelation time and the typical period of oscillations of individual monomers about their metastable minima. We focus on the case where the film is in a glassy state and static friction is observed ($`P>\tau _0/\alpha `$).

In equilibrium ($`v=0`$), the distribution of monomer displacements is consistent with thermal motion about metastable minima plus extremely slow diffusion of monomers between metastable minima. When the top wall is advanced at a low velocity, almost all of the monomers are locked in local metastable minima during any short time interval. As in the Tomlinson model, they provide almost all of the force on the top wall and thus determine $`F_k`$. Monomers that advance by significantly more than thermal displacements are far fewer in number and contribute a smaller force per atom to $`F_k`$. For example, with $`v=0.0005\sigma /t_{LJ}`$ at $`k_BT/ϵ=0`$ only 0.3% of the monomers move by more than 0.12$`\sigma `$ in 2$`t_{JL}`$. These monomers contribute less than 0.1% to the lateral force on the walls, and the percentage decreases as the velocity is lowered.

When the walls are commensurate, the distribution of metastable states depends on the relative positions of the walls. The number of atoms that pop between metastable states is a periodic function of the wall displacement. No pops occur when the walls are aligned directly above each other, and the total energy of the system is minimized in this configuration. The number of pops rises rapidly as the walls are moved out of alignment, and there are pronounced clusters of pops at unfavorable wall positions. The static friction is the maximum force exerted on the top wall at any position, while the kinetic friction represents an average over all wall positions. As a result, the two can be very different and the ratio depends on many factors including the length of the chains, number of chains, etc..

When the walls are incommensurate, the distribution of metastable states is independent of wall position. Thus the kinetic and static friction are correlated. The analogy to the Tomlinson model would suggest that the kinetic and static friction should be equal in this limit. As illustrated above, this ratio depends on the velocity at which $`F_k`$ is calculated and it may also depend on how $`F_s`$ is calculated. Quoted results for $`F_s`$ were obtained by determining the minimum force needed to produce a displacement of $`\sigma `$ in a time $`100t_{LJ}`$. One might thus expect the results to correspond to an effective velocity of order 0.01$`\sigma /t_{LJ}`$, and the static friction is indeed between the values of $`F_k`$ for $`v=0.005`$ and 0.05$`\sigma /t_{LJ}`$. Rate-state models (Eqs. 5 and 6) also imply that the static friction depends on the time spent at zero velocity before a lateral force is applied. We are currently exploring the effect of waiting time on static friction and its relation to the velocity dependent kinetic friction.

## 4 Summary and Conclusions

We have examined the kinetic friction of surfaces separated by a layer of adsorbed molecules as a function of velocity, temperature, interaction potentials and surface alignment. The results are qualitatively consistent with many macroscopic measurements, and can be understood by analogy to simple model systems.

One of our main results is that the kinetic friction rises logarithmically with velocity. This dependence is commonly observed in experimental systems and is often associated with thermal activation . Our analysis of the dynamics of individual monomers is consistent with this idea. We found that the kinetic friction was given by the force exerted by monomers that were trapped in metastable states. A simple argument shows that this force will decrease as $`v`$ decreases.

The force from the top wall on a monomer in a given local potential well is zero at the bottom of the well. As it is displaced away from the minimum the force rises. When the force reaches its maximum value the monomer becomes unstable and pops to a new well. Since atoms that exert the largest force are closest to becoming unstable, they are the ones that are most likely to be thermally activated. When they hop to a new well, they contribute a smaller lateral force. As $`v`$ decreases the time for thermal activation increases as $`\sigma /v`$. Thus more and more of the large contributions to $`F_k`$ are lost due to thermal activation. As a result, $`F_k`$ drops as $`v`$ decreases.

Another important result is that adsorbed layers naturally lead to a kinetic friction that satisfies Eq. 1. He et al. had previously established this relation for the static friction . Thus adsorbed layers provide a molecular basis for Bowden and Tabor’s explanation of Amontons’ laws in both contexts.

Changes in the kinetic friction with interaction potential and temperature closely parallel changes in the static friction. Both $`F_k`$ and $`F_s`$ are insensitive to parameters that are not controlled in typical experiments. Examples include wall orientation ($`\theta `$), the number and precise size of adsorbed molecules, and even the strength of adsorption ($`ϵ_{wf}`$). This insensitivity is consistent with the fact that variations in the measured friction between nominally dry surfaces are relatively small ($`25`$%). In contrast, simulations and analytic studies of bare surfaces show dramatic variations with surface orientation .

The linear relation between shear stress and pressure in Eq. 1 can be understood from a simple hard sphere picture (Fig. 4). The wall atoms create a surface of closest approach that adsorbed molecules must climb up in order to move laterally. The lateral force and normal force are related by the local slope of this surface, and $`\alpha `$ represents some average of the slope. Since $`ϵ_{wf}`$ has little affect on the hard sphere diameter, changing it has little affect on $`\alpha `$. In contrast, changing the ratio $`\sigma _{wf}/d_{nn}`$ has a direct effect on both the slope and $`\alpha `$ (Figs. 3 and 4).

The simple hard sphere model also explains the insensitivity of $`\alpha `$ to the number of adsorbed molecules per unit surface area $`\rho `$ . Increasing the density of molecules spreads the normal pressure $`P`$ and lateral stress $`\tau `$ over more molecules. However, the ratio between them, $`\alpha `$, is still determined by the slope of the surface of closest approach. This surface is determined solely by the solids when the film is less than a monolayer thick. One only expects variations in $`\alpha `$ when the film becomes thicker than a monolayer and can shear internally. Preliminary studies of the thickness dependence do show such effects as the thickness increases to two or three layers.

The close connection between our results for kinetic and static friction is typical of experiments on nominally dry surfaces and surfaces that are purposely coated with thin hydrocarbon films . We discussed how this connection can be understood from results for the simple Tomlinson model. Although the ratio $`F_k(v)/F_s`$ is relatively constant for a given velocity, there are noticeable changes with $`v`$ due to the logarithmic variation in $`F_k(v)`$. We expect that the value of $`F_s`$ will also vary logarithmically with the measurement time, and this remains an avenue for future research.

Another important focus of future research will be to improve the quantitative comparison between theory and experiment. The simulations described here use simple interaction potentials and surface geometries that capture qualitative features of the friction due to adsorbed layers. More detailed models of molecular structure and bonding, as well as surface roughness and elasticity will be needed to allow direct comparison to experiment.

###### Acknowledgements.

The authors wish to thank Drs. M. H. Müser and P. M. McGuiggan for useful discussions. Support from National Science Foundation Grant No. DMR-0083286 and from Intel Corporation through the donation of workstations is gratefully acknowledged.


text:
# Bond-disordered Anderson model on a two dimensional square lattice — chiral symmetry and restoration of one-parameter scaling

## I Introduction

A quantum particle moving in a random potential undergoes Anderson localization quantum phase transition in three dimensions with increasing of the strength of disorder . The order parameter characterizing the localized phase is the inverse localization length $`\xi ^1`$ , describing the exponential decay of the envelope of eigenstates. When the critical point is being approached from the localized phase, localization length, which depends for a given energy only on the strength of the disorder, increases with decreasing of the disorder strength and finally diverges as a power-law at a particular disorder strength. Further decrease of disorder strength then makes the eigenstate extended throughout the whole system. Simultaneously, on length scales smaller than the localization length, eigenstates exhibit multifractal scaling behavior characterized by anomalous scaling of the inverse participation numbers (for definitions and references, see Sec. V).

This basic phenomena, together with the work of Thouless on the scaling of conductance in finite-size systems , led to the scaling theory of localization , one of the main consequences of which is the absence of extended states in two dimensional disordered systems, with two dimensions being the lower critical dimension of the transition. If the spin-orbit interaction is present, however, picture changes and two dimensional systems from symplectic ensemble exhibit localization transition even in two dimensions as opposed to systems from orthogonal ensemble which have all states localized . Presence of the strong magnetic field in two dimensional disordered systems, on the other hand, leads to a completely different behavior – the integer quantum Hall effect (IQHE) – where critical states are present at the middle of each of disorder-broadened Landau levels .

Another class of models exhibiting localization properties different from the systems mentioned above are systems with chiral (particle-hole) symmetry. Such systems are defined on a bipartite lattice with only hopping (off-diagonal or bond-) disorder. Wegner was first to realize the importance of this symmetry in disordered systems , and even one-dimensional systems with this symmetry are known to have peculiar properties, such as diverging DoS at the band center , where the eigenstate decays as $`\mathrm{exp}(\gamma \sqrt{r})`$ , in contrast to one-dimensional site-disordered systems which have DoS bounded  and all states localized.

There are several models with chiral symmetry that have been extensively studied. The simplest two, in the sense that only one orbital per site and nearest neighbor hopping are included, time-reversal symmetry present and the spin not relevant, are the Anderson bond-disordered model (ABD)  and the random Dirac fermion model (RDF) . The main difference between these two models is that, in the non-disordered case, ABD model has a line of points as the Fermi surface at half-filling while RDF has a point Fermi surface and linear dispersion of energies.

This work is concerned with the ABD model on a square lattice of size $`L`$ and periodic boundary conditions, defined by the Hamiltonian:

$$H=ϵ_0\underset{i,j}{}\left(t_{i,j}c_i^{}c_j+\text{H.c}\right),$$
(1)
where brackets denote neighboring sites on the lattice, $`c_i`$ is annihilation operator of the electron at site $`i`$, and $`t`$’s are uniformly distributed random variables $`t_{i,j}(12w,1)`$, with $`0<w1`$. They represent random hopping energies between nearest neighbors, expressed in units of energy $`ϵ_0`$, which is set to $`1`$ hereafter.

Interest in this model mainly comes from its unusual scaling properties at the band center, where Soukoulis et. al  have found critical state at the band center using Green’s function and transfer matrix method (TMM) . More recent TMM calculation by Eilmes et. al confirmed this result with a higher accuracy and showed the validity of one-parameter scaling not too close to the band center. Nevertheless, Miler and Wang  have found in their study of two models with chiral symmetry an apparent band of extended states near the band center, and it remained unclear what is the fate of these states in the infinite 2D system. Another study  yet showed that scaling exponent of the average participation number changed discontinuously near the band center, and the explicit dependence of this energy on the system size proposed by authors implied the existence of another diverging length scale in the problem. The last effect is rather subtle to calculate and led to a different participation number scaling exponent of the ABD model at the band center in Ref.  compared to the one calculated here, as discussed in detail in Sec. V below. Furthermore, Brouwer et al.  have calculated conductance distribution of quantum wires described by (1), and showed its non-universality and necessity to introduce an additional microscopic parameter.

It is thus goal of this paper to present a detailed study of the scaling of localization length on the approach to the band center for an infinite 2D square lattice, and test the validity of one-parameter scaling, as well as to calculate the multifractal properties of the electron probability density on length scales smaller than the localization length. Also a new analytical expression for the DoS of the infinite two-dimensional system near the band center is proposed. This paper is organized as follows: Some general properties and exact results are presented in Sec. II; Calculation of DoS is presented and analyzed in Sec. III; Sec. IV contains analysis of level spacing distributions between the nearest neighbors; Multifractal properties of eigenstates are studied in Sec. V; Scaling of the Lyapunov exponents of transfer matrices of long strips and the scaling of localization length are studied in Sec. VI; and, finally, Sec. VII summarizes results of this work.

## II Some general properties of lattice Hamiltonians with chiral symmetry

Suppose that the lattice is composed of two sublattices $`A`$ and $`B`$ with, respectively, $`N_A`$ and $`N_B`$ sites. The corresponding bond-disordered Hamiltonian with chiral symmetry then has the form:

$$H=\underset{iA,jB}{}\left(t_{i,j}c_i^{}c_j+\text{H.c}\right),$$
(2)
It is easy to show that for every eigenstate $`|\psi `$ with energy $`E`$ there is an eigenstate with energy $`E`$ with a wavefunction that has the opposite sign at each site of one of the two sublattices.

If the total number of cites $`N=N_A+N_B`$ is odd and open boundaries condition is applied (in order to keep the symmetry), then, since all eigenstates come in the opposite energy pairs, there will be exactly one state with eigenenergy $`0`$. This can be furthermore generalized, and if $`m=N_AN_B>0`$, there exist exactly $`m`$ zero-energy eigenstates that have vanishing amplitude on the sublattice $`B`$ .

On the other hand, if $`m=0`$, the electron has equal probability of occupying each of the two sublattices. To show this, (2) is represented in the basis where the first and second half of basis vectors are eigenstates of the position operator on sites of sublattice $`A`$ and $`B`$, respectively. The Hamiltonian is then represented as

$$H=\left(\begin{array}{cc}0& M\\ M^{}& 0\end{array}\right),$$
(3)
where $`M`$ is a square matrix of hopping elements from one sublattice to the other. Eigenstate $`|\psi =\left(\begin{array}{c}|\psi _A\\ |\psi _B\end{array}\right)`$ satisfies:

$$E=\psi |H|\psi =\mathrm{\hspace{0.17em}2}\text{Re}\psi _A|M|\psi _B.$$
(4)
On the other hand,

$$H|\psi =\left(\begin{array}{c}M|\psi _B\\ M^{}|\psi _A\end{array}\right)=\left(\begin{array}{c}E|\psi _A\\ E|\psi _B\end{array}\right).$$
(5)
From (4) and (5) now follows that, for ABD model, $`\psi _A|\psi _A=1/2`$.

In this work only even $`L`$ finite size systems on a square lattice with periodic boundary conditions are studied, because one of the main goals of this work is to understand the vicinity of the critical point of ABD model on the infinite square lattice, which in turn has $`m=0`$, while the limit $`L\mathrm{}`$ for odd $`L`$ and open boundary conditions has $`m=1`$.

## III Density of states near the band center

Density of states (DoS) is calculated by exact numerical diagonalization of finite size Hamiltonians for various $`L`$ for many configurations of disorder, and binning of eigenenergies. Obtained DoS for each system size, $`\rho _L(E)`$, are normalized to $`1`$. The $`L`$ dependent parts of such obtained $`\rho _L(E)`$ are then removed leaving $`L`$ independent DoS $`\rho (E)`$, which is therefore expected to be correct in the $`L\mathrm{}`$ limit. The removal of finite size dependency is based on an observation that DoS converges quickly away from the band center with increasing of $`L`$. Thus, only a small number of eigenenergies (up to 20) closest to the band center and corresponding DoS histograms has been calculated for each $`L`$. The calculated $`\rho _L(E)`$ plotted on a single graph revealed that three bins closest to the band center are where the system size dependence sets in. Their removal thus led to the DoS $`\rho (E)`$ of the infinite system.

Results for $`\rho _L(|E|)`$ are given in Fig. 1, for system sizes $`L=10,20,\mathrm{},60`$ and number of disordered configurations ranging, respectively, from $`160000`$ to $`4100`$. They are fitted to a powerlaw divergence $`\rho _L(E)=C_L|E|^{\alpha _L}`$, as well as to the Gade’s result :

$$\rho _L(E)=C_L\frac{1}{|E|}\mathrm{exp}\left(\kappa _L\sqrt{\mathrm{ln}|E|}\right).$$
(6)
All the calculations were done for several different number of bins, and the obtained values of fitting parameters were the same within error bars.

The figure 1 shows that (6) describes $`\rho _L`$ very accurately for $`L40`$, including the size-dependent part. The powerlaw, on the other hand, also describes the data accurately for same system sizes, but fails to describe the $`L`$-dependent part of $`\rho _L`$. Despite this, neither of the two forms describe the whole $`\rho (E)`$ accurately. Instead, the expression found to best fit the obtained $`L`$-independent DoS, given in Fig. 2, is

$$\rho (E)=C\frac{1}{\sqrt{|}E|}\mathrm{exp}\left(\kappa \sqrt{\mathrm{ln}|E|}\right),$$
(7)
with $`\kappa =1.345\pm 0.005`$ and $`C=1.30\pm 0.03`$, represented by the full line in the same figure. The observed range in which (7) is accurate is for all the energies studied smaller than $`6\times 10^2`$.

## IV Distribution of the nearest-neighbor level spacings

In the localized regime, an eigenstate is determined mainly by a local configuration of disorder where the wavefunction is localized, and two eigenstates close in energy are spatially far apart. Level repulsion is therefore absent and the distribution of the nearest neighbor level spacings $`sE_{i+1}E_i`$ is Poissonian ,

$$D_P(s)=\frac{1}{\delta }\mathrm{exp}\left(\frac{s}{\delta }\right).$$
(8)
where $`\delta s`$ is the mean level spacing.

In the delocalized phase, on the other hand, eigenstates are extended throughout the system and level repulsion becomes significant for eigenstates with close energies. In the infinite 3D Anderson site-disordered model (ASD),

$$H=\underset{i}{}ϵ_in_i\underset{i,j}{}\left(c_i^{}c_j+c_j^{}c_i\right),$$
(9)
with uniformly distributed $`ϵ_i(W/2,W/2)`$, distribution of level spacings becomes that of Gaussian orthogonal ensemble (GOE) , very accurately described by the Wigner surmise

$$D_W(s)=\frac{\pi }{2}\frac{s}{\delta ^2}\mathrm{exp}\left(\frac{\pi }{4}\left(\frac{s}{\delta }\right)^2\right).$$
(10)
In finite-size systems, localized states are on average at a distance $`L`$ rather than infinitely far apart. This leads to a repulsion between adjacent energy levels and non-universal distribution $`D_L(s)`$ . Shklovskii et. al.  have shown that $`D_L(s)`$ of the 3D site-disordered Anderson model exhibits linear dependence on $`s`$ characteristic for $`D_W(s)`$ for small $`s`$ and exponential tail characteristic for $`D_P(s)`$ for large $`s`$. They were able, from the finite size scaling analysis of the tail, to accurately determine the critical point and exponent. In the infinite size limit, they have recovered not only $`D_P(s)`$ in the insulating phase and $`D_W(s)`$ in the conducting phase, but also a system size independent non-universal distribution at the critical point which was furthermore shown by Braun et. al. to be dependent on boundary conditions. This method was also used for an accurate determination of localization length in two dimensional ASD model , confirming the absence of delocalized states following the scenario of the insulating phase from Ref.  described above.

To see the effect of symmetry of the Hamiltonian (2) on the distribution of level spacings, let us for a moment consider $`i`$-th eigenenergy $`E_i`$ of the ASD model. Upon averaging over disorder, the $`E_i`$ will be distributed between $`E_i^{min}`$ and $`E_i^{max}`$ according to some distribution. Some of the eigenenergies, for $`i`$ close to $`N/2`$, will have $`E_i^{min}<0<E_i^{max}`$. This is, however, forbidden for eigenstates of the ABD model since every $`E_i`$ of (2) is negative for $`i<N/2`$ and positive for $`i>N/2`$. This means that eigenenergies of (2) close to the band center are effectively pushed away from it due to the symmetry. If $`L`$ is much smaller than the localization length, states will be repelled among themselves due to their large spatial overlap. But the two states closest to the band center, being simply related to each another by the symmetry, will not repel at all, i.e. the state closest to the band center is at the (high energy) end of the spectrum. Thus, these two states are distributed around zero, where distributions of all other individual levels go to zero.

To explore consequences of this simple analysis, the level spacing distribution is calculated between each pair of adjacent levels separately. Let us denote by $`D_i(s)`$ level spacing distribution between $`i`$ and $`i+1`$ energy level after unfolding of spectrum , i.e. expressing level energies in units of mean level spacing, where $`0<E_1<E_2<\mathrm{}<E_{N/2}`$. Figure 3 shows $`D_1(s),\mathrm{},D_5(s)`$, for $`L=20,40`$ and, respectively $`150000,120000`$ configurations, and it can be seen that $`D_1(s)`$ is distinctly different than $`D_2(s),\mathrm{},D_5(s)`$. The same effect was also present for $`L=10,30`$, while for $`L=50`$ and $`60`$ the number of disorder configurations was insufficient for accurate enough determination of individual $`D_i(s)`$ . This illustrates how the presence of chiral symmetry can profoundly influence spectral characteristics near the band center, despite the fact that DoS of ABD and ASD models seem to have same shape for adequately chosen pairs of disorder parameters $`w`$ and $`W`$ away from the band center (and after rescaling of $`ϵ_0`$) .

## V Multifractality of eigenstates

Eigenstate of an electron in the random potential fluctuates from site to site and it was proposed that the eigenstate at the mobility edge in disordered systems in general should have fractal structure, and shown that even localized states in one and two dimensions exhibit fractal character on length scales smaller than the localization length .

Inverse participation numbers $`Z_q`$ (IPN) are particularly convenient quantities to describe scaling properties of probability distribution of the electron. IPN of an eigenstate $`\mathrm{\Psi }`$ are defined as:

$$Z_q(\mathrm{\Psi })\underset{i=1}{\overset{L\times L}{}}|\mathrm{\Psi }(𝐫_i)|^{2q}.$$
(11)
Intuitively, their meaning can be seen by looking at the participation number $`Z_2(\mathrm{\Psi })^1`$ : it is equal to 1 for a state localized at one site and to $`N`$ for plane-waves. Participation number thus gives generally the number of sites at which the wavefunction is significantly different than zero. Participation numbers $`Z_q(\mathrm{\Psi })^1`$ then generalize this by giving the number of sites where probability distribution of electron is very high (for large positive $`q`$’s), very low (for large negative $`q`$’s), and in between these extrema, continuously parameterized by $`q`$.

More convenient, with the advantage of being defined as averages over disorder at a given energy $`E`$, are IPN defined as functions of $`E`$ and system size $`L`$,

$$Z_q(L,E)Z_q(\mathrm{\Psi })\delta (E(\mathrm{\Psi })E),$$
(12)
where the brackets denote averaging over disorder. $`Z_q(E,L)`$ can be numerically calculated by averaging (11) over all eigenstates from $`M`$ configurations of disorder belonging to an energy interval of width $`\mathrm{\Delta }E`$ around $`E`$, and studying the limit $`\mathrm{\Delta }E0`$ for large $`M`$ .

Wegner  pioneered this kind of investigations, and Castellani and Peliti  proposed that eigenstates near the critical point are multifractal on length scales smaller than $`\xi `$. The most important feature of IPN of eigenstates is their scaling with system size and energy :

$`Z_q(L,E)`$ $``$ $`L^{\tau _q},`$ (13)
$`Z_q(L,E)`$ $``$ $`|EE_c|^{\pi _q},`$ (14)
where $`E_c`$ is the critical energy. The former scaling is present at any $`E`$ for $`L\xi (E)`$, while the latter holds in the critical region of the transition .

Within the framework of multifractality , electron probability density is characterized by several quantities that can be derived from $`\tau _q`$ — the generalized dimension $`D_q`$ and the singularity strength $`\alpha _q`$ of the $`q`$-th singularity with the fractal dimension $`f_q`$ :

$$(q1)D_q\tau _q,\alpha _q\frac{d\tau _q}{dq},f_q(\alpha _q)\alpha _qq\tau _q.$$
(15)
$`D_q`$ represents generalization of the fractal dimension, and it is constant and equal to the fractal dimension for ordinary fractals, while $`f_q(\alpha _q)`$ is the singularity strength spectrum describing multifractal as an interlaced set of fractals with fractal dimensions $`f_q`$, where the measure on the $`q`$-th fractal scales as a powerlaw with exponent $`\alpha _q`$. These quantities have several general properties: $`D_0`$ is the fractal dimension of the support (two in this work); $`D_1`$ is called information dimension since it describes scaling of the entropy of the measure, and there exist finite $`D_{min}=D_q\mathrm{}`$ and $`D_{max}=D_q\mathrm{}`$.

This work is concerned mainly with the spectrum of generalized dimensions $`D_q`$ characterizing the spatial structure of eigenstates on the length scales smaller than the localization length, while the properties of $`\pi _q`$ will be discussed elsewhere. Before detailed discussion of the results, an overview of the main results of this part of the paper is given. Scaling of IPN at the band center is calculated first, and shown that scaling properties (that is, whole spectrum of generalized dimensions $`D_q`$) changes discontinuously near the band center, at an energy $`E^{}(L)`$ for a range of $`L`$ studied. It is then shown that $`E^{}`$ can be quite accurately identified with the half of the width of the energy range around band center within which two states occur on average in the ensemble of disordered systems. Existence of this energy reveals the existence of a length scale $`\xi ^{}(E)`$, diverging when $`E0`$, that is the system size at which IPN change their scaling properties from one powerlaw dependence on $`L`$ to another. This change is then explained as a finite size manifestation of the critical point, and the critical exponent calculated by a finite-size analysis.

Calculation of $`Z_q(E,L)`$ starts with calculation of $`Z_q(E,L,\mathrm{\Delta }E)`$, which is just IPN averaged over all eigenstates from an energy interval $`(E\mathrm{\Delta }E/2,E+\mathrm{\Delta }E/2)`$, taken from $`N_\mathrm{\Omega }`$ realizations of disorder, followed by studying the limit $`\mathrm{\Delta }E0`$ . Results at the band center for system $`L=80`$ and several different $`q`$’s, are presented in Fig. 4. The error bars in the figure are taken to be the standard deviation of average value.

The figure suggests the existence of an energy $`E^{}`$ independent of $`q`$ (and therefore defined by the whole multifractal measure) such that decreasing $`\mathrm{\Delta }E`$ below $`E^{}`$ does not change $`Z_q(E,L,\mathrm{\Delta }E)`$ significantly. Decreasing of $`Z_q`$ to a smaller extend, however, is still present for $`\mathrm{\Delta }E<E^{}`$, and the main source of this is the mismatch between average and typical value of IPN at a given energy. Thus, the effect should become smaller as the number of disorder configurations that are averaged over is increased, and $`Z_q(E,L,\mathrm{\Delta }E<E^{})Z_q(E,L)`$ up to the corresponding statistical error. This can be seen in Fig. 4 and 5, where values of $`Z_q`$ for $`\mathrm{\Delta }E<E^{}`$ are approximately constant within the error bars (as indicated by the horizontal dashed lines), while for $`\mathrm{\Delta }E>E^{}`$ there is approximately linear dependence of $`Z_q`$ on the bin size $`\mathrm{\Delta }E`$. In this sense Figure 5 suggests that $`E^{}`$ exists for all the systems studied (and can be shown to be independent of $`q`$ for each of them analogously as shown in Fig. 5).

Analogous analysis is carried out for energies away from the band center and $`Z_q(E,L)`$ determined accordingly, where it turns out that the convergence for these energies is slightly easier to establish and occurs at larger $`\mathrm{\Delta }E`$ than at the band center. Such obtained energy intervals used in calculations of IPN for different $`E`$ were small enough so that there was practically no overlap among them.

The vertical dashed lines in Fig. 4 and Fig. 5 represent half of the energy $`E_2(L)`$, defined as the width of the energy interval around zero in which every system from the ensemble of disordered systems has two states on average,

$$1=L^2_{E_2/2}^0\rho _L(ϵ)𝑑ϵ.$$
(16)
From the figures one can see that $`E^{}E_2/2`$ for all system sizes studied except $`L=40`$ (the smallest system studied, not shown in Fig. 5), where the convergence seems to be somewhat slower. Equation (16) defines a new length scale $`\xi ^{}(E)`$ that can be described as the system size $`L`$ for a given energy $`E`$ such that the number of states within the energy interval between $`E`$ and $`E`$ is two on average. This can be defined as

$$\xi ^{}(E)=E_2^1(E),$$
(17)
where $`E_2^1(E)`$ is the inverse function of the $`E_2(L)`$, which is defined by (16). Since $`E_2(L)`$ goes to zero when $`L\mathrm{}`$, the new length scale diverges when $`E0`$.

It is tempting to integrate results for $`\rho _L`$ from Sec. III to obtain $`\xi ^{}(E)`$ explicitly. This, however, does not give the correct result since the whole analysis of DoS from Sec. III is done for energies larger than the width of the distribution of two states closest to the band center, which in turn defines $`\xi ^{}`$ in (17). In other words, there is an energy cutoff, vanishing when $`L\mathrm{}`$, below which the fits are not accurate, most obviously seen by noticing that all of the assumed analytical forms of $`\rho _L`$ are diverging at the band center, while the actual $`\rho _L`$ is not.

Nontrivial feature connected with the existence of the new length scale $`\xi ^{}`$ is that scaling exponents $`\tau _q(E)`$ are different for $`L\xi ^{}(E)`$ and $`L\xi ^{}(E)`$. This is a generalization of findings from Ref. , where a deviation from the powerlaw scaling of the average participation number for different $`\mathrm{\Delta }E`$ at the band center appeared whenever $`\mathrm{\Delta }E`$ exceeded $`E^{}(L)`$ in several models with chiral symmetry. It can be straightforwardly shown that same happens with scaling of $`Z_q(E=0,L,\mathrm{\Delta }E)`$. This suggests that a new scaling characteristic should be attributed to the two states closest to the band center, and that $`E^{}`$ scales as $`E_2`$, with a coefficient of proportionality close to one. This also implies that the corrections to the constant $`Z_q(E,L,\mathrm{\Delta }E)`$ for $`\mathrm{\Delta }E<E^{}`$ are small. Therefore, the length at which change of scaling occurs should be close to and depend on the energy proportionally to $`\xi ^{}(E)`$, while the change of scaling should be a narrow crossover, as opposed to much broader crossover from powerlaw to constant IPN that occurs at the termination of multifractal scaling (13) for $`L\xi `$.

This is compatible with the results presented in Fig. 6, which gives calculated $`Z_q(L,E)`$, for $`q=0.9,2`$ as well as for all $`L`$ and $`E`$ studied. Smaller $`q`$’s allow for more accurate determination of the scaling properties, and it can be seen from the data that for both $`q`$’s there are two characteristic scaling behaviors — one occurs at the band center and nearby energies for smaller $`L`$, while the other scaling holds at energies further away from the band center, and for larger $`L`$ at energies close to the band center.

Scaling exponents $`\tau _q`$ are determined form the linear regression of the data for system sizes $`L=50,60,\mathrm{},100`$. Goodness of the powerlaw fit is quantitatively characterized by a coefficient $`\gamma _q(E)`$ next to the quadratic term from an additional quadratic fit of the data. Such obtained $`\tau _q`$ and $`\gamma _q`$ from the data in Fig. 6 are presented in Fig. 7. Results show the existence of three different cases: (i) away from the band center, $`\gamma _q0`$ and $`\tau _q`$ is independent of $`E`$, indicating powerlaw dependence of IPN on $`L`$; (ii) approaching the band center, $`\gamma _q`$ becomes different than zero, indicating that powerlaw is not obeyed. This is due to the emergence of the new scaling for system sizes $`L\xi ^{}(E)`$ discussed above and present in Fig. 6; and (iii) for $`E=0`$ powerlaw is obeyed again ($`\gamma _q0`$), but with a different $`\tau _q`$ than for energies away from the band center.

Discontinuity of $`\tau _q(E)`$ for all the other $`q`$’s studied naturally leads to two different spectra of generalized dimensions, and Fig. 8 shows calculated $`D_q`$ for all energies except the three nonzero energies closest to the band center (which cannot give $`\tau _q`$ from the fitting procedure used here due to the change of scaling properties discussed above). In particular, the participation number grows with the number of sites $`L^2`$ as a powerlaw with exponents $`\beta (E=0)=0.25\pm 0.02`$, and $`\beta (E0)=0.55\pm 0.05`$. This should be compared with the result for ABD model of Ref. , $`\beta (E=0)=0.50\pm 0.06`$. It is easy to explain this discrepancy since the bin size $`\mathrm{\Delta }E=4\times 10^4`$ used in Ref.  was, depending on the system size, roughly over an order of magnitude too large to detect the correct scaling behavior, and therefore $`\beta (E0)`$ was obtained instead.

$`D_q`$ is calculated for only three negative values of $`q`$, mainly because IPN for negative $`q`$’s are determined mostly by the parts of eigenstates with the smallest probability to find the electron, which in turn acquire the highest relative error during numerical diagonalization of the Hamiltonian. Difficulties in calculating $`D_q`$ in this regime even arose suspicion that multifractality might break down for negative $`q`$ . It is thus important to show that $`D_q`$ is defined for negative $`q`$’s as well as for positive ones. Accuracy of all the calculations of IPN done in this section can be straightforwardly improved by increasing the number of configurations of disorder that was averaged over. This would lead to smaller error bars of all the quantities calculated, as well as to the wider range of $`q`$’s for which $`D_q`$ can be calculated.

Results of this section give the following picture of the scaling of IPN with system size: for any energy $`E`$ close enough to the band center, there exist two powerlaw scalings: one for $`L_0<L\xi ^{}(E)`$ described by the set of exponents $`\tau _q(E=0)`$, and another one for $`\xi ^{}(E)L\xi (E)`$, described by a different set of exponents $`\tau _q(E0)`$, which leads to the two different spectra of generalized dimensions $`D_q`$.

Dependence of the new length scale on energy, $`\xi ^{}(E)`$, can be easily determined from (17) by integrating the actual numerical data for $`\rho _L(E)`$, and the result, obtained from Fig. 9, gives $`\xi ^{}(E)|E|^\nu ^{}`$, with

$$\nu ^{}=0.35\pm 0.01.$$
(18)
The meaning of this new length scale and corresponding exponent $`\nu ^{}`$ can be understood by assuming that the additional scaling of the two states is due to the finite size effect of “smearing” of the $`E=0`$ critical point of the infinite system, because exactly the two states closest to the band center become critical when $`L\mathrm{}`$. The critical energy then changes by $`\mathrm{\Delta }ϵ`$ which is, in a finite system of size $`L`$ (in units of the lattice spacing), equal to $`\mathrm{\Delta }E_2(L)L^\delta =L^{1/\nu ^{}}`$. If furthermore $`\xi |E|^\nu `$ in the infinite system, the shift $`\mathrm{\Delta }ϵ_c`$ of the critical energy in the system of size $`L`$ is $`\mathrm{\Delta }ϵ_cL^{1/\nu }`$, from the general theory of the finite size scaling . Therefore, $`\nu ^{}`$ equals $`\nu `$, the critical exponent of the correlation length of the infinite two-dimensional system. It should be noticed that the finite size scaling applied here is somewhat different than usual, where $`\mathrm{\Delta }ϵ_c|EE_c|/E_c`$ . Here, $`E_c=0`$ and $`\mathrm{\Delta }ϵ_c|EE_c|`$, where both energies are expressed in units of $`ϵ_0`$, as discussed in the Introduction. Thus, energy $`ϵ_0`$ appears in the denominator of $`\mathrm{\Delta }ϵ_c`$ rather than the critical energy itself.

## VI Localization length near the band center

In order to calculate the localization length of the infinite two dimensional system, finite-size scaling analysis of MacKinnon and Kramer  (FSS) is applied to TMM of Pichard and Sarma . In this analysis, inverse Lyapunov exponents (ILE) of the transfer matrix of the ABD model  of a long quasi one dimensional strip of width $`M`$ are calculated for several energies near the band center and one parameter scaling analysis applied to the largest ILE, from which the correlation length of 2D system is calculated. The scaling analysis consists of assuming that the change of the largest ILE, $`\mathrm{\Lambda }(E,M)`$, due to rescaling $`MbM`$ can be compensated by an appropriate change of energy, after which $`\mathrm{\Lambda }`$ will remain the same, which implies that

$$\mathrm{\Lambda }(E,M)=\mathrm{\Lambda }(M/\xi (E)).$$
(19)
Figure 10 shows the calculated $`\mathrm{\Lambda }(E,M)`$ for several energies $`E`$ close to the band center and for various strips up to $`128`$ sites wide. The figure also gives the second largest (dashed line) renormalized ILE, $`\mathrm{\Lambda }_2`$, for the two lowest non-zero energies studied ($`E=10^5`$ and $`10^6`$) and for $`M20`$. All values are obtained with relative error of $`1\%`$ or better.

At $`E=0`$ and for $`M`$ even, all ILE become doubly degenerate due to the presence of chiral symmetry , and they scale linearly with $`M`$ for $`M16`$, reflecting the scale invariance of $`\mathrm{\Lambda }`$ characteristic of a critical state. To see the effect this degeneracy has on scaling properties of ILE, we should recall that ILE of transfer matrices of disordered systems repel each other in general  and become self-averaging quantities for sufficiently long stripes . The symmetry now enforces “dimerization” of pairs of ILE, acting as an effective attractive force between each pair of ILE that become degenerate at zero energy, as in Fig. 10, where, for a fixed $`M`$, $`\mathrm{\Lambda }`$ and $`\mathrm{\Lambda }_2`$ are closer together for the smaller energy. The largest ILE thus decreases in a strip of width $`M`$ on approaching the band center, as in models with chiral symmetry studied in Ref. . On the other hand, at any given energy, every pair of ILE becomes more repelling with increasing of $`M`$. Such increased number of ILE thus diminishes the attractive effect of the symmetry, and, depending on the relative strengths of the two effects, there are two regimes: (i) for smaller $`M`$, the largest ILE increases approximately linearly with $`M`$, and, since the slope of the rise changes with energy, FSS cannot be done; (ii) for sufficiently large $`M`$, on the other hand, the repulsion among ILE due to disorder dominates and FSS is possible.

The one-parameter universal function $`\mathrm{\Lambda }(M/\xi (E))`$, obtained for system sizes $`M50`$ for energies $`|E|>10^4`$, and for system sizes $`M64`$ for energies $`10^5|E|10^4`$, is presented in Fig. 11. The obtained localization length $`\xi (E)`$, in units of the localization length of the smallest energy studied ($`E=0.1`$), is shown in the inset of the same figure. Calculated $`\xi (E)/\xi (0.1)`$ is fitted to the powerlaw for energies $`10^5E10^3`$ , and the result suggests powerlaw diverging localization length at the band center, with the exponent

$$\nu =\mathrm{\hspace{0.17em}0.335}\pm 0.034,$$
(20)
in agreement with the result (18) of Sec. V.

In Ref. , the authors also calculated one-parameter scaling curves for energies further away from the band center than in this work and found that there is a universal curve for a large set of ABD systems with different $`w`$ that they studied. The obtained one-parameter function suggested that $`\mathrm{\Lambda }(x)`$ either grows indefinitely when $`E0`$, which is similar to the two dimensional ASD model  or becomes constant. If this result were correct, however, for the smaller energies than studied there, it would imply that localization length at the band center either diverges logarithmically, as in the 2D ASD model, or that it has a finite limiting value. We see therefore that, although correct in the energy range studied in Ref. , the picture of a single scaling curve for different $`w`$ from the same work should break down for smaller energies, since the powerlaw diverging $`\xi `$ requires that $`d\mathrm{ln}\mathrm{\Lambda }/d\mathrm{ln}M`$ goes to zero when $`E0`$, and separate one-parameter scaling curves exist for each $`w`$.

Some additional analysis of the results can be done by introducing $`\mathrm{\Lambda }_{max}(M)`$, defined as the maximal $`\mathrm{\Lambda }(E,M)`$ for a given strip width $`M`$. Importance of this quantity comes from the fact that points where $`\mathrm{\Lambda }`$ reaches its maximum for various energies cannot be described by FSS, and, at the same time, $`\mathrm{\Lambda }_{max}`$ limits by its definition possible values that $`\mathrm{\Lambda }(x)`$ (obtained from the FSS) can have. The main observation is that $`\mathrm{\Lambda }_{max}`$ seems to grow slower than linearly in Fig. 10. Linear (or faster) growth would imply existence of an additional energy scale below which FSS would break down for all large $`M`$ and $`\mathrm{\Lambda }(E,M)`$ would grow linearly and indefinitely with $`M`$ in the same figure, implying the existence of a whole band of extended states where one-parameter scaling would not hold. Slower than linear growth of $`\mathrm{\Lambda }_{max}`$ seen, therefore, furthermore suggests the existence of a single critical point at the band center of the system on an infinite square lattice.

## VII Conclusions

In conclusions, main results of this work are summarized. Density of states of finite size systems are calculated and shown the validity of Gade’s expression (6). The calculated part of the DoS for the system on an infinite square lattice, however, suggests a different dependence on energy near the band center,

$$\rho (E)=C\frac{1}{\sqrt{|}E|}\mathrm{exp}\left(\kappa \sqrt{\mathrm{ln}|E|}\right),$$
(21)
with $`\kappa =1.345\pm 0.005`$ and $`C=1.30\pm 0.03`$.

Other calculated quantities share in common a qualitative feature of a discontinuous change near or at the band center. The nearest neighbor level spacing distributions, for instance, between the state closest to the band center and the next one seem to be distinctly different than the level spacing distributions between other neighboring states. This was argued to be connected to the chiral symmetry of the model, which places the two states closest to $`E=0`$ at the (high energy) end of the spectrum. These two states are furthermore found to play a crucial role in explaining the discontinuous change of the scaling properties of IPN near the band center. Extrapolation to the limit of infinite square lattice then led to the two different spectra of generalized dimensions $`D_q`$ — one for $`E0`$ (present at the length scales smaller than the localization length), and another one for $`E=0`$ (present at all length scales). Finite size scaling associated with the effects that finite $`L`$ has on the critical band center states of the infinite system led to the value $`\nu =0.35\pm 0.01`$ of the critical exponent of the localization length.

Additional scaling analysis of ILE of transfer matrices of long quasi-1D systems gave the value $`\nu =\mathrm{\hspace{0.17em}0.335}\pm 0.034`$. This suggests that there is only one critical state at the band center with all other states localized in the system on the infinite square lattice, in agreement with findings of Ref. . One parameter universal function $`\mathrm{\Lambda }(x)`$ is calculated and another discontinuity found, since $`\mathrm{\Lambda }(x(E0))\mathrm{\Lambda }(E=0,M)`$.

A puzzling feature of the critical exponent $`\nu `$ of the ABD as well as of the RDF model  is their apparent disagreement with the rigorous theorem of Chayes and Chayes et. al , which states that $`\nu 1`$ in two-dimensional quantum disordered systems in general. This question, however, requires further study and will be addressed elsewhere.

###### Acknowledgements.

This work is partially supported by the Department of Physics and Astronomy at Michigan State University. Author is thankful to S. D. Mahanti for suggestions on improving the manuscript and to S. A. Trugman, R. Bhatt, T. A. Kaplan, V. Zelevinski, M. Mostovoy, I. Herbut, D. Mulhall, B. Nikolić and M. Milenković for inspiring discussions.