text:
# Kinetic Friction and Atomistic Instabilities in Boundary-Lubricated Systems

## I Introduction

The every-day phenomenon friction is of great practical and economical importance, which is one of the motivations to improve our understanding of tribological processes persson98 ; bhushan99 . Friction between two solids differs from that between a solid and a fluid in that both static and kinetic friction appear finite, while the force between a solid and a fluid vanishes linearly with sliding velocity $`v_0`$ at small $`v_0`$. Static friction $`F_\mathrm{s}`$ is the externally applied force necessary to initiate relative sliding motion between two solids, whereas kinetic friction $`F_\mathrm{k}`$ is the force needed to maintain the sliding motion. Phenomenological friction laws, which date back to da Vinci, Amontons, and Coulomb dowson79 , often provide a good description on the macroscopic scale.

The microscopic origins of kinetic friction are still a matter of debate, even though it has long been recognized that kinetic friction must be due to dynamical instabilities brillouin09 ; caroli98epjb . While there can be many different processes leading to instabilities, they all have in common that potential energy is converted abruptly into kinetic energy and ultimately lost as heat muser03acp . Although instabilities can occur on many different time and length scales, there has been an enhanced interest in identifying those that occur on atomic scales. This quest is not only motivated by the miniaturization of technical devices down to the nanometer scale, but also by the desire to better understand macroscopic friction. The understanding of single-asperity contacts is needed as basis for the full description of macroscopic friction, where the bulk-mediated coupling between contacts gives rise to additional effects.

Load-bearing, simple-asperity contacts are often in the order of microns. According to Hertzian contact mechanics and generalizations thereof, the pressure is rather constant in the contact with the exception of the areas close to the circumference, where pressure gradients are large. In the centre of the contact, most of the lubricant is squeezed out. One may assume that these boundary-lubricated areas often account for most of the energy dissipation when two solids are slid against each other, unless the solids are very compliant, in which case elastic instabilities may also contribute a significant amount to the net dissipation. If wear was the main source of friction, material would have to rub off from the surfaces much faster than observed experimentally gnecco02prl . Hydrodynamic lubrication would likewise result in values for friction orders of magnitude too small, if it were assumed to be the dominant dissipation process.

While the crucial role of surfactants for friction has certainly been recognized, relatively little attention has been paid to characterize dynamical instabilities in boundary lubricants. Most of the work on instabilities leading to friction is devoted to elastic processes, which are most simply described in the Prandtl-Tomlinson (PT) model prandtl28 ; tomlinson29 . In the PT model, an atom is pulled over a substrate by a spring that moves at constant velocity $`v_0`$. If the spring stiffness is below a certain critical value, the atom’s instantaneous velocity can exceed $`v_0`$ by many orders of magnitude, see i.e. the discussion in Refs. caroli98epjb, ; muser03acp, . This process results in non-vanishing $`F_\mathrm{k}`$ in the limit of zero $`v_0`$ as long as thermal fluctuations are absent. There is, however, a crucial difference between instabilities in boundary lubricants and instabilities occurring in elastic manifolds that are modelled in terms of the PT model and related approaches such as the Frenkel-Kontorova model frenkel38 ; braun98 . In boundary lubricants, atoms are only weakly connected to each other and to the confining walls. As a consequence, bond breaking can occur, whereas in elastic models, bonds are treated as unbreakable. This seemingly subtle difference leads to different tribological behavior.

Two different avenues have been pursued in the recent past to study dynamics in boundary lubricants and its consequences for tribological properties. One is a minimalist approach, in which one single lubricant atom embedded between two shearing plates is considered rozman96prl ; rozman96pre . In the following, we will refer to this approach as the impurity limit. The other avenue incorporates a large ensemble of lubricant atoms thompson90sci ; thompson90pra . This approach can eventually include surface curvature and elastic deformation of the surfaces making it possible to study what effect the interplay of surface curvature and elastic deformations have on dry or boundary-lubricated friction wenning01 ; persson02jcp .

In this paper, we intend to analyze what features of simplistic models appear robust as the level of complexity in the description of the boundary lubricant is increased. Since kinetic friction is intimately connected with instabilities, we focus on the analysis of instabilities. In a precedent paper by one of us muser02prl , it was found that the existence of instabilities in the impurity limit and as a consequence the friction-velocity relationship $`F_\mathrm{k}(v_0)`$ depends on the ‘details’ of the model. For instance, it was found that for 1D, commensurate interfaces, the sign of the first higher harmonic in the lubricant-wall potential determines: (a) whether or not the athermal kinetic friction remains finite in the zero-velocity limit, and (b) the exponent $`\beta `$ that describes the finite-velocity corrections by

$$F_\mathrm{k}(v)F_\mathrm{k}(0)v_0^\beta .$$
(1)
Note that Eq. (1) changes its form when thermal noise is included into the treatment, i.e., it becomes linear at small velocities muser02prl . Depending on the ratio of the relevant energies and temperature, thermal effects may be negligible down to very small values of $`v_0`$.

While Ref. muser02prl,  is mainly focused on the impurity limit, we intend to extend the analysis in a systematic fashion to less idealized situations. For example, instead of simple spherical impurities, dimers and hexamers (6-mers) will be studied. Moreover direct interactions between lubricant particles will be included and the effects of increasing coverage will be discussed. The central assumption of our analysis is the existence of instabilities or ‘pops’ of certain degrees of freedom. A pop is a sudden, seemingly erratic motion of a particle (or a collective degree of freedom) characterized by a velocity much larger than the associated thermal velocity or the drift velocity of the atom. Pops heat the lubricant or alternatively they couple directly to the confining solid walls, i.e., by inducing phonons in the walls. They will eventually induce more dramatic effects such as generation of dislocations or abrade the surfaces. However, as argued above, these extreme processes are rare and hence presumably they are not responsible for the main part of the energy dissipation. This is the motivation to concentrate on the energy transfer to the phonon bath that is due to elementary process in the lubricant. The underlying idea of the presented approach can be described as follows. Sliding-induced instabilities make the velocity probability distribution (PD) of the lubricant atoms deviate from the thermal equilibrium PD. This alters the balance of energy flow from and to the lubricant. The energy missing in this balance is provided by the external driving device.

In this paper, we will develop a simple kinetic theory that connects the energy dissipation with the velocity PD (Section II). After discussing the numerical techniques in Sec. III, we will apply the theory to models of boundary lubrication of various complexity. This will include both the impurity limit, which is discussed in Section IV, and more complex situations that include interaction between lubricant atoms (Section V). Section VI contains the conclusions.

## II Theory

### II.1 General Comments

The most fundamental assumption in this paper is that the interaction between the lubricant atom $`i`$ and the confining wall can be decomposed into one conservative part $`V_\mathrm{w}(𝐫_i)`$ and one non-conservative term consisting of a damping force plus thermal noise. $`V_\mathrm{w}(𝐫_i)`$ depends only on the difference between the position $`𝐫_i`$ and the positions of top wall $`𝐫_\mathrm{t}`$ and bottom wall $`𝐫_\mathrm{b}`$. It can be written as

$$V_\mathrm{w}(𝐫_i)=V_\mathrm{b}(𝐫_i𝐫_\mathrm{b})+V_\mathrm{t}(𝐫_i𝐫_\mathrm{t}),$$
(2)
where depending on the model under consideration, the vectors $`𝐫`$ can be one-, two-, or three-dimensional. Unless otherwise noted, the relative motion of the walls is imposed externally, i.e., by constant separation (or constant load) and constant relative velocity $`v_0𝐞_x=(\dot{𝐫}_\mathrm{t}\dot{𝐫}_\mathrm{b})`$ of the walls parallel to the sliding direction indicated by the unit vector $`𝐞_x`$. We assume the normal pressure variations to be small, which means that the coupling to each individual confining (crystalline) wall is periodic parallel to the interface, i.e., it is periodic in the $`xy`$-plane.

In the theoretical part of our treatment, we assume that the non-conservative force $`𝐅_{\mathrm{t},i}^{(\mathrm{nc})}`$ that a wall (here the top wall) exerts on the lubricant atom $`i`$ consists of a simple, viscous damping term $`\gamma _\mathrm{t}(\dot{𝐫}_i\dot{𝐫}_\mathrm{t})`$ plus a random force $`𝚪_\mathrm{t}(t)`$, thus

$$𝐅_{\mathrm{t},i}^{(\mathrm{nc})}=m_i\gamma _\mathrm{t}\left(\dot{𝐫}_i\dot{𝐫}_\mathrm{t}\right)+𝚪_\mathrm{t}(t),$$
(3)
where the random force $`𝚪(t)`$ is chosen such that detailed balance is obeyed when the external stress is absent. Thus the usual $`\delta `$ correlation of the random forces is assumed, namely risken84

$$𝚪_\mathrm{t}(t)𝚪_\mathrm{t}(t^{})=2D\gamma _\mathrm{t}m_ik_\mathrm{B}T\delta (tt^{}),$$
(4)
where $`D`$ is the physical dimension, $`k_\mathrm{B}T`$ is the thermal energy, and $`m_i`$ the mass of lubricant atom $`i`$. Random forces plus damping term $`m_i\gamma _\mathrm{t}𝐯`$ (we dub the sum thermostat) mimic the interactions with phonons and/or other excitations, which are not treated explicitly. Typically, the time scales associated with these excitations are short compared to the motion of a particle from one minimum to another, which justifies the assumption of $`\delta `$ correlated random forces for our purposes. Of course, damping can and will be different normal and parallel to the interface. However, this detail does not have any significant consequences for the conclusions presented in this paper. Similarly, the explicit treatment of internal elastic deformations does not alter the major conclusions either.

We will now be concerned with the derivation of a formal equation for the friction force. In any steady-state of the system, the average force on the upper wall (or the lower wall) must be zero. If the time average was different from zero, the upper wall would be accelerated in contradiction to the steady-state assumption, as pointed out for instance, by Thompson and Robbins thompson90pra . The net force on the upper wall consists of three contributions: The externally applied force $`𝐅_{\mathrm{ext}}`$, the conservative force between lubricant and wall, and the non-conservative force done by the thermostat. The external force $`𝐅_{\mathrm{ext}}`$ does the work $`W_{\mathrm{ext}}`$ on the upper wall given by

$$W_{\mathrm{ext}}=𝑑𝐫_\mathrm{t}𝐅_{\mathrm{ext}}.$$
(5)
The kinetic friction force $`F_\mathrm{k}`$ follows from that expression as it equals the work done on the system by external forces divided by the distance moved.

Since the conservative potential is assumed to be periodic, it cannot do any net work on a steady-state system and we may not consider it in our energy balance. This implies that the work must be done on the thermostat. The power dissipated into the damping term is proportional to $`\gamma m(\dot{𝐫}\dot{𝐫}_\mathrm{t})^2`$, however, parts of that contribution can be provided by the stochastic random force $`𝚪(t)`$. Hence, if we want to account only for the power $`P_{\mathrm{ext}}`$ that is dissipated into the damping term due to the externally applied force, we have to integrate over the velocity distribution $`P(v)`$ but we have to subtract the contribution that is due to the random force. The latter contribution is very difficult to calculate. We assume that this heat flow from the random force into the impurity system is identical to that in thermal equilibrium, in which the equilibrium (Maxwell Boltzmann) distribution $`P_{\mathrm{eq}}(v)`$ applies. This yields

$$P_{\mathrm{ext}}=N_{\mathrm{fl}}\gamma m_0^{\mathrm{}}v^2\left\{P(v)P_{\mathrm{eq}}(v)\right\}𝑑v,$$
(6)
where $`N_{\mathrm{fl}}`$ is the number of lubricant atoms, $`v`$ is the velocity of a particle relative to the center-of-mass motion of the upper wall. The net external driving force (or in other words the kinetic friction force $`F_\mathrm{k}`$) can now be associated with

$$F_\mathrm{k}=P_{\mathrm{ext}}/v_0.$$
(7)
We want to emphasize that Eqs. (5) through (7) allow one to calculate friction forces under more general conditions than those of our particular model, for instance, if the thermostat only acts on the atoms in the outermost layers of the walls as e.g. employed in Ref. braun01 . The approach can also be extended in a straightforward manner if generalized forms of the thermostat are employed such as in dissipative particle dynamics espanol95  or if the thermostat is based on a Mori Zwanzig formalism mori65 ; zwanzig65 . The main limitation of Eq. (6) in the present context is that effects due to heating of the walls are not included. Again, a minor modification would allow one to include heating of the walls into the presented framework as well. However, as we will mainly focus on small velocities, the mentioned effects will be small and shall be neglected in the following.

Note that an alternative way of determining the friction force in the steady state is to time average the conservative plus the non-conservative force that the upper wall exerts on the lubricant. The observation that the work done by the conservative force is essentially zero does not imply that its time-average must be zero. A formal derivation of the conclusions from this section is given in the appendix.

### II.2 Effect of Instabilities

As discussed in the introduction, the externally imposed relative motion of the confining walls may induce sudden, dynamic instabilities or ‘pops’ during which the particles’ velocities greatly exceed both their thermal velocities and the relative sliding velocity $`v_0`$ of the walls. This means that at a time $`t+\delta t`$ the atom does not find a stable position in the $`𝒪(v_0\delta t)`$ vicinity of the old stable position at time $`t`$. The continuous trajectory ends at $`t`$ and the particle has to move to the next mechanically stable position to resume its path. The particle will then pop into the next local potential minimum and for low sliding velocities, its peak velocity $`v_{\mathrm{peak}}`$ will be solely determined by the energy landscape and consequently $`lim_{v_00^+}v_{\mathrm{peak}}/v_0`$ diverges. Its kinetic energy will be dissipated (e.g. by phonons) and lead to friction. This process will lead to a deviation of the velocity distribution $`P(v)`$ from the thermal equilibrium distribution $`P_{\mathrm{eq}}(v)`$ valid for $`v_0=0`$. Fig. 1 shows such instabilities for a model system that is described in detail in Sect. III.

The velocity distribution $`P(v)`$ and hence the friction force $`F_\mathrm{k}`$ can be calculated in principle, once the precise form of the lubricant’s interaction is known. Risken’s book on the Fokker-Planck equation risken84  gives an excellent overview of methods that allow one to treat models like ours, namely externally-driven systems that are mainly deterministic but also contain a certain degree of thermal noise. An analytical approach remains difficult in our case, due to the potential’s complex time dependence. Therefore a different, phenomenological approach will be pursued.

An instability will invoke a trajectory during which potential energy is abruptly converted into kinetic energy. The kinetic energy will then be dissipated into the thermostat, i.e., the phonon bath of the confining walls. After some time, which depends on the coupling strength to the thermostat, the Maxwell-Boltzmann probability distribution (PD) will be resumed, provided no new instability has been invoked in the meantime. An instability will thus create a typical velocity PD that will show up as a tail in the Maxwell-Boltzmann PD. Unless the two confining walls are identical and perfectly aligned (thus commensurate), there is a class of instabilities in which the energy lost during the ‘pop’ shows a broad distribution, see also Ref. daly02cm, . Every pop, characterized for instance by the energy dissipated, will contribute to $`P(v)`$ in its own way. We assume that the net sum $`P_{\mathrm{tail}}`$ of all these individual tails shows exponential dependence on velocity, thus

$$P_{\mathrm{tail}}\mathrm{exp}\left(B𝐯\overline{𝐯}\right),$$
(8)
where $`\overline{𝐯}`$ is the average velocity of the impurity under consideration, typically $`\overline{𝐯}=v_0𝐞_x/2`$, and $`B`$ is a constant. The motivation for this particular choice of $`P_{\mathrm{tail}}`$ partly stems from Jaynes’ principle of information theory jaynes82 . From it follows that the most likely normalized PD on $`[0,\mathrm{})`$ with given mean value about which we do not have more knowledge is the exponential distribution HonerkampJaynes , thus an exponential ansatz for $`P_{\mathrm{tail}}`$ is plausible. Further restraints lead to deviations from an exponential form. While these considerations are heuristic, our choice of $`P_{\mathrm{tail}}`$ happens to be a quite accurate description for the velocity PDs of impurities between 2D, incommensurate surfaces. This will be demonstrated later in the result section.

At small sliding velocity $`v_0`$, the statistical weight of the tails must increase linearly with velocity. Hence, the normalized PD function for the $`x`$ component is given by

$`P(v_x)`$ $`=`$ $`\sqrt{{\displaystyle \frac{m}{2\pi k_\mathrm{B}^{}T}}}\left(1{\displaystyle \frac{2A^{}v_0}{B^{}}}\right)\mathrm{e}^{m(v_x\overline{v}_x)^2/2k_\mathrm{B}T}`$ (9)
$`+A^{}v_0\mathrm{e}^{\mathrm{B}^{}|\mathrm{v}_\mathrm{x}\overline{\mathrm{v}}_\mathrm{x}|}.`$
Here $`A^{}`$, and $`B^{}`$ are phenomenological parameters that can (and will) depend on the externally applied load $`L`$ that an impurity has to counterbalance, damping $`\gamma `$, and other parameters. However, they should depend only weakly on temperature $`T`$ and sliding velocity $`v_0`$ at small $`T`$ and small values of $`v_0`$. This is because $`A^{}v_0`$ is a measure for the rate of the fast processes (which should be proportional to $`v_0`$ at small $`v_0`$), while $`B^{}`$ characterizes the instability related velocity PD. For 1D systems $`B^{}=B`$, the projection of a 2D exponential PD on one axis however leads to different $`A^{}`$ and $`B^{}`$ note\_on\_PD . Inserting Eq. (9) into Eqs. (6) and (7) and integrating over $`v_x`$ yields the following friction force per impurity atom $`F_\mathrm{k}/N_{\mathrm{fl}}`$:

$$\frac{1}{N_{\mathrm{fl}}}F_\mathrm{k}=4\gamma m\frac{A^{}}{B^3}2\gamma \frac{A^{}}{B^{}}k_\mathrm{B}T,$$
(10)
from which the friction coefficient $`\mu _\mathrm{k}=F_\mathrm{k}/(N_{\mathrm{wall}}L)`$ follows. (Here, $`L`$ is the average load carried per atom belonging to an outermost wall atom, thus $`F_\mathrm{k}`$ and $`N_{\mathrm{wall}}L`$ represent respectively the net friction force and the net load.) Of course, Eq. (10) can only be valid as long as Eqs. (8) and (9) give an accurate description of the non-equilibrium velocity PD and provided that the heat flow from the thermostat into the impurities is close to the thermal equilibrium heat flow. At extremely small $`v_0`$, two arguments show that the assumption of exponential tails cannot persist. First, the energy $`\mathrm{\Delta }E_{\mathrm{diss}}`$ that is dissipated during a pop has an upper bound, which in turn implies an upper bound for the peak velocity. Second, close to equilibrium, thermal noise is sufficient to invoke (multiple) barrier crossing and recrossing. The ratio of sliding and noise-induced instabilities becomes small, which in turn makes the non-equilibrium corrections be less significant.

Eq. (10) is based on the one-dimensional distribution functions $`P(v_x)`$. For two-dimensional interfaces and boundary lubricants, the relevant PD is $`P(v_{})`$ with $`v_{}=\sqrt{v_x^2+v_y^2}`$. We want to note in passing that accumulating the histogram $`P(v_{})`$ contains the same required information as $`P(v_x,v_y)`$, however, it requires less data storage. This is why for 2-dimensional boundary lubricants, we monitor $`P(v_{})`$ instead of $`P(v_x,v_y)`$. Assuming rotational symmetry, Eq. (9) can be replaced with

$$P(v_{})=2\pi v_{}Av_0e^{Bv_{}}+\left(1\frac{2\pi Av_0}{B^2}\right)P_{\mathrm{eq}}(v_{}),$$
(11)
where $`A`$ and $`B`$ are phenomenological coefficients with similar meanings as their counterparts $`A^{}`$ and $`B^{}`$ in the one-dimensional description.

Inserting Eq.(11) into Eq. (6) yields:

$$\frac{1}{N_{\mathrm{fl}}}F_\mathrm{k}=12\pi \gamma m\frac{A}{B^4}4\pi \gamma \frac{A}{B^2}k_\mathrm{B}T$$
(12)
The parameters $`A`$ and $`B`$ will be obtained by fitting the PDs accumulated during MD simulations. If such fits turn out to be good approximations (and they do for incommensurate surfaces), then one has to ask into question descriptions of frictional interfaces that are based on local effective temperature (and local pressure). Local effective temperatures would imply Gaussian rather than exponential velocity tails. As we will argue later, this difference might matter with regard to chemical reactivity in a frictional interface.

Of course, the fits will never be perfect, and one has to address the question whether one can obtain information experimentally on $`A`$ and $`B`$, i.e., by measuring $`F_k`$ and the average kinetic energy of the lubricant atoms. We will therefore compare the values of $`F_k`$ that are calculated directly (by averaging the force on the top wall) with those that are obtained indirectly with Eq. 12 after $`A`$ and $`B`$ are obtained through fits.

## III Main Model and Methods

In this paper, we analyze the trajectories of atoms and molecules embedded between two walls in relative sliding motion by means of molecular dynamics. Different models with varying degree of complexity are investigated ranging from rather simple, 1-dimensional impurity models to 3-dimensional systems, in which the interaction between the lubricant particles are taken into account. In the latter case, lubricant particles are not only simple atoms but may also represent short polymers. Here, we will only describe the methods relevant to the full 3-dimensional simulations, as all other cases only require ’dumbed-down’ versions of that method or small alterations thereof, such as suppressing the interaction between lubricant atoms.

In our model, lubricant atoms interact with each other and with wall atoms via a truncated Lennard Jones potential

$$V_{\mathrm{LJ}}(r)=\{\begin{array}{cc}\hfill 4ϵ\left[\left(\sigma /r\right)^{12}\left(\sigma /r\right)^6\right]+C,& r<r_\mathrm{c}\hfill \\ \hfill 0,& rr_\mathrm{c}\hfill \end{array},$$
(13)
where $`r`$ is the distance between two atoms. $`ϵ`$ defines the energy scale and $`\sigma `$ the length scale of the system. Both quantities are set to unity. A constant value $`C`$ is added to the potential for inter-atomic distances smaller than the cut-off radius $`r_c`$, which assures the continuity of the potential. $`r_\mathrm{c}`$ was chosen as the minimum of the LJ potential, $`r_\mathrm{c}=r_{\mathrm{min}}=2^{1/6}`$, unless otherwise noted. This choice corresponds to a purely repulsive interaction and can be justified by the observation that at large pressures the essential behavior is caused by the repulsion of the particles. The main effect of including the attractive LJ contribution in the present context would be to add an adhesive pressure. Throughout the paper quantities are measured using LJ units, such as time in units of $`t_{\mathrm{LJ}}=(m\sigma ^2/ϵ)^{1/2}`$ and forces in units of $`ϵ/\sigma `$. Atomic masses $`m`$ are also set to unity.

Both top and bottom wall lie in the $`xy`$ plane and consist of $`N_{\mathrm{wall}}`$ discrete atoms arranged in the hexagonal (1,1,1) plane geometry of an fcc crystal. The nearest neighbor spacing $`d_{\mathrm{nn}}`$ in the walls is $`1.20914\sigma `$ unless noted otherwise. This choice of $`d_{\mathrm{nn}}`$ does not match with other length scales in the system. Moreover, the relatively large value for $`d_{\mathrm{nn}}`$ enhances the effect of surface corrugation. In most simulations presented here, the normal load is kept constant. The load $`L`$ will be stated in terms of normal load per atom in the upper wall, hence a unity normal load corresponds to a pressure of about 0.79 in reduced quantities.

Commensurate wall geometries were realized by orienting the two walls in parallel, whereas incommensurability was achieved by rotating the upper solid surface by 90. The use of periodic boundary conditions in the wall plane required a slight distortion of the perfect hexagonal geometry in order to obtain two quadratic walls. Therefore, walls were not perfectly incommensurate anymore, but quasi-incommensurate (as every setup realized with finite number precision, strictly speaking). A wall unit cell consists of two atoms, at positions $`(0,0)`$ and $`(d_{\mathrm{nn}}/2,\sqrt{3}d_{\mathrm{nn}}/2)`$. By choosing the ratio of the wall unit cells in $`\widehat{x}`$ and $`\widehat{y}`$ close to the ideal value $`\sqrt{3}`$ this distortion was minimized. We do not use other relative wall rotations in the full 3-dimensional simulations, as it was found in the study of similar models that the influence of the rotation angle is weak if it exceeds $`5^{}`$ he01tl . We note that when two solids in an experiment come in contact they will most likely be incommensurate, as it would take utmost care to have two identical defect-free crystals and orient them perfectly. As detailed calculations show, elastic deformations do not generally alter this argument provided the solids are treated as 3-dimensional objects muser03acp ; shinjo93 .

While our analysis is focused on simple fluids, we include some work on small chain molecules in order to study aspects of molecular lubricants. To this end, we used a well established bead-spring model proposed by Kremer and Grest KremerGrest1990 , which models individual monomers as LJ particles while chain connectivity is ensured by a FENE (finitely extensible nonlinear elastic) potential given by

$$V_{\mathrm{FENE}}=\frac{1}{2}k_{\mathrm{ch}}R_{\mathrm{ch}}^2\mathrm{ln}\left[1(r/R_{\mathrm{ch}})^2\right],$$
(14)
with $`R_{\mathrm{ch}}=1.5\sigma `$ and $`k_{\mathrm{ch}}=30ϵ/\sigma ^2`$ KremerGrest1990 . Typical values for hydrocarbons are $`ϵ30\text{meV}`$, $`\sigma 0.5\text{nm}`$, resulting in a typical time scale of $`t_{\mathrm{LJ}}3\text{ps}`$ KremerGrest1990 .

Simulations were done using a fifth order Gear-predictor corrector algorithm (see e.g. Haile1997 ) with an integration step of $`dt=0.005`$. To maintain constant temperature, a stochastic Langevin thermostat was employed AllenTildesley . It consists of ideal white noise random forces and damping forces acting on all thermostatted particles, which obey the fluctuation-dissipation theorem AllenTildesley . A damping constant $`\gamma =0.5`$ was used in all simulations. In the presence of instabilities, the precise choice of $`\gamma `$ is usually quite irrelevant for friction forces at small velocities fisher85 , i.e., for the choice $`\gamma _\mathrm{b}=0.5`$ and $`\gamma _\mathrm{t}=0`$ we may expect similar friction forces as for the perhaps more natural choice of $`\gamma _\mathrm{b}=\gamma _\mathrm{t}=0.25`$, if $`v_0`$ is sufficiently small. However, in either case, one must ensure that the random forces satisfy the fluctuation dissipation theorem.

## IV Impurity Limit

### IV.1 1D model systems

Here, we want to discuss and extend those results from Ref. muser02prl,  that are relevant to this study. In Ref. muser02prl, , the following potential $`V_\mathrm{t}`$ for the interaction between the impurity and the (one-dimensional) top wall was employed

$`V_\mathrm{t}`$ $`=`$ $`V_{\mathrm{t},0}\mathrm{cos}(2\pi (xx_\mathrm{t})/b_\mathrm{t})`$ (15)
$`+`$ $`V_{\mathrm{t},1}\mathrm{cos}(4\pi (xx_\mathrm{t})/b_\mathrm{t}),`$
$`b_\mathrm{t}`$ being the lattice constant of the upper wall. A similar impurity-wall coupling is used to describe the interaction for the bottom wall, however, the indices $`t`$ (for top) have to be replaced with $`b`$ (for bottom). This is a generalization of the interactions suggested in Refs. rozman96prl,  and rozman96pre,  in that a non-zero first higher harmonic $`V_{\mathrm{t},1}`$ is considered in Ref. muser02prl . Moreover, the lattice constant of the bottom wall $`b_\mathrm{b}`$ is allowed to differ from that of the top wall $`b_\mathrm{t}`$. We note in passing that Refs. rozman96prl,  and rozman96pre,  were concerned with the interplay between external driving and embedded system, while Ref. muser02prl,  focused on the constant-velocity friction in such systems.

In Ref. muser02prl, , it was found that the behavior of the steady-state, low-velocity kinetic friction is surprisingly rich. For instance, athermal, zero-velocity friction turned out to vanish for $`V_{\mathrm{t},1}=V_{\mathrm{b},1}0`$. In that case, athermal, small-velocity friction can be described as a powerlaw $`F_\mathrm{k}v_0^\beta `$ with a non-universal exponent $`0<\beta <1`$. For $`V_{\mathrm{t},1}=V_{\mathrm{b},1}>0`$, $`F_\mathrm{k}`$ remains finite as $`v_0`$ approached zero. The case of $`V_{\mathrm{t},1}=V_{\mathrm{b},1}=0`$ is particularly intriguing, as the minima move at constant velocity $`v_0/2`$, then at some points in time (when there is perfect deconstructive interference of $`V_\mathrm{t}`$ and $`V_\mathrm{b}`$), the location of the potential energy minima make a ‘phase jump’ by a distance $`b_\mathrm{t}/2`$. This phase jump, however, does not result in significant energy dissipation, as the location from where the impurity stems and the location where the impurity ends up are symmetrically equivalent.

For incommensurate walls ($`b_\mathrm{t}b_\mathrm{b}`$), the behavior is even richer. If the first higher harmonic is not included, one wall exerts a maximum force on the impurity and drags the impurity along. As a consequence, $`F_\mathrm{k}`$ is linear in $`v_0`$, which we call Stokes friction. For one certain value of the first higher harmonic $`V_{\mathrm{t},1}^{}`$ (at a fixed ratio $`b_\mathrm{s}=b_\mathrm{b}/b_\mathrm{t}`$), $`F_\mathrm{k}`$ can best be described as a power law in the limit of small $`v_0`$. For $`V_{\mathrm{t},1}>V_{\mathrm{t},1}^{}`$, $`F_\mathrm{k}`$ remains finite in the limit $`v0`$, again provided thermal fluctuations are absent. Fig. 2 shows a friction diagram for the 1D, impurity-lubricant model. In order to yield a more complete picture than that given in Ref. muser02prl , additional calculations have been carried out, in order to determine $`V_{\mathrm{t},1}^{}`$ for different values of $`b_\mathrm{s}=b_\mathrm{b}/b_\mathrm{t}`$.

For a more detailed discussion of the effect of thermal fluctuations and the occurence of a Stokes friction regime in one dimension and its potential applications, we refer the reader to Ref. muser02prl, .

### IV.2 2D model systems

#### IV.2.1 Model details

We now allow the lubricant atoms to move within the $`xy`$-plane, but motion normal to the interface in $`z`$ direction is still neglected. Eq. (15) must then be replaced with a new model potential. As in other studies, we consider the symmetry of the confining walls to be triangular, i.e., (111) surfaces of an fcc crystal, for which the potential $`V_\mathrm{t}`$ between top wall and an impurity can be written as:

$$V_\mathrm{t}(𝐱,z)=\underset{𝐠}{}\stackrel{~}{V}_\mathrm{t}(𝐠,z)\mathrm{exp}\left\{i𝐠(𝐱𝐱_\mathrm{t})\right\}.$$
(16)
In Eq. (16), the sum goes over all two-dimensional reciprocal lattice vectors of the trigonal lattice $`𝐠`$, $`𝐱`$ is the position of the lubricant particle in the $`xy`$ plane and $`𝐱_\mathrm{t}`$ the in-plane position of the (top) wall. $`z`$ denotes the (fixed) distance between (top) wall and impurity.

The Fourier coefficients $`\stackrel{~}{V}_\mathrm{t}(𝐠,z)`$ between chemically non-bonding species often decay exponentially with $`|𝐠|`$ and increasing distance from the surface $`z`$, thus $`\stackrel{~}{V}_\mathrm{t}(𝐠,z)`$ can be written as

$$\stackrel{~}{V}_\mathrm{t}(𝐠,z)=\stackrel{~}{V}_\mathrm{t}(𝐠,0)\mathrm{exp}\{\alpha |𝐠|z\},$$
(17)
where both parameters $`V_\mathrm{t}(𝐠,0)`$ and $`\alpha `$ depend on the chemical nature of impurity atom and confining wall. Potentials of this form are known as Steele potentials steele73 . They have proven to describe the potential energy landscape of atoms on crystalline surfaces reasonably well Patrykiejew2000 . The fundamental harmonic in this potential is related to the smallest non-zero lattice vector $`𝐠=(2\pi /\sqrt{3},0)`$ and its five symmetrically equivalent counterparts, which are obtained by rotating $`𝐠`$ successively by 60. (Note that the distance between ’atoms’ in the walls is set to unity, which differs from the choice for the full 3D simulation model, see Sect. III.) First higher harmonics are related to reciprocal lattice vectors that are the sum of a suitable pair of two different fundamental $`𝐠`$’s and so on. The fundamental harmonic will be dominant at small loads, however, as the external pressure increases (which makes $`z`$ decrease), the relative importance of higher harmonics will increase due to Eq. (17). The coupling between impurities and bottom wall is similar to that between impurities and top wall, however, the reciprocal vectors $`𝐠`$ are rotated an angle $`\theta `$ with respect to the top wall’s $`𝐠`$’s. Two walls are called commensurate if $`\theta `$ is an integer multiple of 60. An equivalent 2D model without higher harmonics, was used recently by Daly et al. daly02cm  for a study similar to that presented here.

In the following, we will be concerned with an analysis of mechanically stable position for the impurity atoms and their motion as the walls slide against each other. The goal is to identify situations, where the trajectory of a mechanically stable position suddenly disappears, which would lead to a dynamical instability. Such an analysis was given for 1D lubricants in in Ref. muser02prl , see Fig. 1 in that paper, and also in Ref. daly02cm  for 2D systems. For our analysis, the bottom wall’s lateral position $`𝐱_\mathrm{b}`$ is kept fixed, while $`𝐱_\mathrm{t}`$ is moved in small constant increments $`d𝐱`$ with $`|d𝐱|dx=10^5`$ to $`10^2`$. After identifying an initial relative minimum in the impurity wall potential $`V_{\mathrm{i},\mathrm{w}}=V_\mathrm{b}+V_\mathrm{t}`$, a steepest descent algorithm (we used the Mathematica function FindMinimum\[ \] MathematicaFindMinimum ) searches for the closest minimum in $`V_{\mathrm{i},\mathrm{w}}`$ in the vicinity of the previous one. If the distance $`\mathrm{\Delta }d`$ between the location of the new and the old minimum is greater than $`\mathrm{\Delta }d=0.1`$, we say that we identify a pop. While this choice of $`\mathrm{\Delta }d`$ is somewhat arbitrary, we ensured that our conclusions remained unaltered when $`\mathrm{\Delta }d`$ was varied in reasonable bounds and $`dx`$ was further decreased. Since both methods employed in this study (simulations and steepest descent) are exact and identical within controllable errors, their results mutually agree within these margins.

#### IV.2.2 Commensurate Walls

Impurity atoms between commensurate walls only have a finite number of non-equivalent minima in their potential energy landscape. Once a minimum is identified, symmetrically equivalent minima will exist at periodically repeated positions that follow from the lattice of the confining walls. Various mechanically stable ’stacking’ geometries can be envisioned for our walls of trigonal symmetry, for instance hexagonal close packed (hcp) and face cubic centered (fcc) type configurations best characterized as respectively $`ABA`$ and $`ABC`$ layering structures. The boundary lubricant reflects the middle layer. While it does not correspond to an ideally crystalline layer, the probability for a lubricant atom to sit at a certain position would be indeed periodic, i.e. it would have a maximum in every single $`B`$ position.

As the two walls are slid with respect to each other, the situation is akin of the relative sliding of two commensurate, one-dimensional surfaces muser02prl . The ‘trajectories’ of mechanically stable positions bifurcate and recombine at certain relative, lateral displacements of the two solids, see Fig. 1a in Ref. muser02prl . This scenario invokes so-called continuous instabilities. The peak velocities during continuous instabilities tends to zero as $`v_0`$ tends to zero, however, this does not happen linearly. This will ultimately lead to the following behavior if the walls are slid parallel to a symmetry axis and the top wall is allowed to move freely in transverse direction: After every half lattice constant moved, the system will convert from an fcc type structure to an hcp type structure or vice versa. This behavior is also found in our 3D default system (see Fig. 3), in which impurities interact with wall atoms through Lennard Jones potentials rather than through Steele potentials. The situation changes, when the top wall is not allowed to move in the transverse direction, which was the choice in Ref. daly02cm, . In that case, particles will occupy positions with high potential energy, which leads to instabilities. In this study, however, we focus on the case of zero transverse force.

As the mechanically stable positions of the embedded impurities show no discontinuities, the kinetic friction force will tend to zero at small $`v_0`$ even if thermal fluctuations are absent. From the comparison to the 1D model systems, one would expect a powerlaw behavior as in Eq. (1) with $`F_\mathrm{k}(0)=0`$. This behavior does not depend on the sliding direction. It is also observed for ratios of $`V_\mathrm{t}(𝐠,z)/V_\mathrm{b}(𝐠,z)`$ different from unity.

A central issue in the present paper is the question how robust the property of the simple impurity model is as more complexity is added to the model. In the present case, one may argue that lubricant atoms would be able to move in a correlated fashion up to a coverage of one monolayer. Above this coverage, the impurity model breaks down for obvious reasons.

#### IV.2.3 Incommensurate Walls

Impurity atoms between incommensurate walls have an infinite number of inequivalent minima in their potential energy landscape for a given relative wall displacement. This means that at a given moment in time, it is impossible to find two different positions where the value of the potential and all its derivatives are identical. Yet, the number of inequivalent trajectories of (meta)stable positions can be small, because in most cases, they will all be identical up to temporal shifts when the walls are in relative sliding motion. See also the discussion of the dynamics of the incommensurate Prandtl-Tomlinson model by D. S. Fisher fisher85 .

We analyze the instabilities by varying randomly the relative orientation $`\theta `$ between the two walls as well as the sliding direction $`\varphi `$. At this point, we are only concerned with the occurrence of instabilities, rather than with the (average) amount of energy dissipated during an instability. Instabilities between incommensurate walls are shown in Fig. 1.

Unless $`\theta `$ is close to an integer multiple of 60, we find that the number of instabilities depends only weakly on $`\theta `$ and $`\varphi `$. If we chose the fundamental harmonics $`\stackrel{~}{V}(𝐠)`$ of both walls to be identical and the higher harmonics to be absent, then we find on average one instability each time the upper wall has been moved laterally with respect to the lower wall by a distance of $`200`$ $`d_{\mathrm{nn}}`$. Increasing the interaction strength for just one wall does not change the behavior until the ratio $`\stackrel{~}{V}_\mathrm{t}(𝐠)/\stackrel{~}{V}_\mathrm{b}(𝐠)`$ or its inverse exceeds about 4.7. Above this threshold value, the metastable positions and hence the particles follow the motion of just one wall and no instability occurs.

Like Daly et al. daly02cm , we note that the instabilities are possible due to transverse motion of the impurities, see also Fig. 1. One of the issues Daly et al. also discussed was the question above which value of $`\stackrel{~}{V}_\mathrm{t}(𝐠)/\stackrel{~}{V}_\mathrm{b}(𝐠)`$ the lubricant particle remains pinned to the (top) wall. They reported a value of 4.5, while we find a slightly higher value of 4.7, which essentially confirms the prediction. Furthermore, we also analyzed the effects of first higher harmonics, which were neglected in Ref. daly02cm, . Including the first higher harmonic $`𝐠_1`$ in addition to the fundamental harmonics $`𝐠_0`$ increases the number of instabilities, in particular for higher harmonics with positive sign. Thus, the occurrence of instability remains a robust feature of incommensurate walls. For a ratio $`V(𝐠_1)/V(𝐠_0)=0.1`$, the number of instabilities is increased by a factor of six. (At this ratio the absolute value of the second harmonic would be in the order of 0.01 $`V(𝐠_0)`$, see Eq. (17), and will thus be neglected.) One may argue that the observed increase in pops is related to an increase of incommensurability due to an additional (small) length scale.

### IV.3 3D Model

#### IV.3.1 Effect of commensurability on PDs

We now turn to the analysis of the full three-dimensional model, described in detail in Sect. III. Here, we also include the interaction between the lubricant atoms. However, as the coverage is only a quarter layer, the results remain almost identical to the ideal impurity limit. Despite these changes with respect to Section IV.2, all arguments discussed there remain valid under the new conditions. For instance, Fig. 3 shows the expected dynamical behavior of two commensurate walls separated by lubricant impurities in sliding motion, i.e., an alteration of hcp and fcc type configurations. Most importantly, the trajectories of lubricant atoms become continuous for commensurate walls. It is instructive to compare Figs. 1 and 3. For completeness, we mention that the simulations in Figs. 1 and 3 were both done at an external load of $`L=30`$ per top wall atom, a thermal energy of $`T=0.01`$, and relative sliding velocity of $`v_0=10^3`$.

The different trajectories of the mechanically stable states result in qualitatively different velocity distributions, even in the presence of thermal fluctuations, which is shown in Fig. 4. It can be seen that the velocity PDs of impurities between incommensurate walls can indeed be described with the non-equilibrium PD suggested in Eq. (9). It turns out that the PDs longitudinal ($`x`$) and transverse ($`y`$) to the sliding direction ($`x`$) are almost identical. We note in passing that the velocity PD normal to the interface ($`z`$ direction) is affected much less than the in-plane PDs.

#### IV.3.2 Effect of sliding velocity and temperature on PDs

As the relative sliding velocity between the walls is changed by a factor of 100, the prefactor of the exponential tail scales with the same factor, as suggested in Eq. (9). The commensurate walls behave differently. First, the non-equilibrium velocity distribution $`P(v)`$ deviates from equilibrium much less than for incommensurate walls and it does not obey Eq. (9) as well. More importantly, the tails of $`P(v)`$ behave differently from those of incommensurate surfaces under a change of sliding velocity. This difference is due to the absence of instabilities for the commensurate system. At $`v_0=10^4`$, the velocity PD for commensurate walls is almost identical to the equilibrium Maxwell Boltzmann PD, while at the same $`v_0`$, the PDs for incommensurate walls show distinct non-equilibrium tails. We employ a logarithmic scale for the PDs, because the tails can hardly be discerned on a linear scale.

Further examination of the distribution functions for incommensurate surfaces as shown in Fig. 5 reveals that the coefficients $`A`$ and $`B`$ in Eq. (11) are approximately constant for a wide range of velocities and temperatures. The parameters can be easily read off the graphs: The slope of the tails equals $`B`$ and the exponential of the $`y`$-axis intercept of a fitted line through the tails equals $`Av_0`$. The data for Fig. 5 were produced with load $`L=10.0`$ for temperatures $`T=10^3\mathrm{}10^1`$ and sliding speeds $`v_0=10^3\mathrm{}10^1`$.

The present discussion is valid when the non-equilibrium tails are clearly visible such as in Fig. 5. It becomes invalid when $`v_0`$ reaches extremely small values, i.e., when the tails are starting to disappear under the central Maxwell-Boltzmann peak. Eq. (9) then ceases to be a good description of the PDs in that limit and Eqs. (9) through (12) are no longer applicable. However, the equation describing the heat-flow balance between thermostat and confined system, Eq. (6), is unaffected by this argument and remains valid even in the limit $`v_00`$.

#### IV.3.3 Effect of load on PDs

The load dependence of the coefficients $`A`$ and $`B`$ was investigated as well. We show the effect of load on the PDs for one of our model systems exemplarily in Fig. 6. Many similar calculations were done for other loads, coverages, and sliding velocities with similar results for incommensurate surfaces. In all cases, we found that $`A`$ is roughly proportional to $`L^{0.8}`$, while $`B`$ is approximately proportional to $`L^{0.4}`$.

From the normalization factor of the central equilibrium peak in Eq. (12), one may infer that the ratio $`A/B^2`$ is a measure for the number of atoms far out of equilibrium and hence for the number of invoked instabilities. Given the proportionalities $`AL^{0.8}`$ and $`BL^{0.4}`$, this number remains constant when $`L`$ is increased. Inserting the proportionalities $`AL^{0.8}`$ and $`BL^{0.4}`$ into Eq. (12) results in a small deviation from Amontons’s law $`F_\mathrm{k}L`$ at the microscopic level.

Potential differences scale with $`L`$ in lowest order, thus we obtain for the energy dissipated in a pop $`\mathrm{\Delta }E_{\mathrm{diss}}Lv^2`$. Hence, for exact proportionality, the width of the non-equilibrium tails was $`L^{0.5}`$, respectively $`BL^{0.5}`$, yielding Amonton’s law. This shows that $`BL^\lambda ,\lambda 0.5`$ is to be expected, while the precise value of the exponent $`\lambda `$ will depend on the specific system potentials.

The deviation in our system is due to a shift of the relative significance of lower- and higher-order harmonics. This shift would presumably be smaller if the repulsive forces were modeled with (slightly more realistic) exponentially repulsive forces muser01prl .

Fig. 6 reveals that the exponential tails fall off less slowly when the pressure is increased. Thus, large pressures in sliding contacts can dramatically increase the probability of large velocities, even though the lubricant’s average kinetic energy $`T_{\mathrm{kin}}`$ (or effective thermal energy) may barely change. This favours the occurrence of rare events such as chemical bond breaking, as it becomes much more likely that a bond is hit quasi-simultaneously by two high-velocity atoms. As the non-equilibrium PDs fall off less slowly than the equilibrium PDs, bond breaking will occur more frequently in non-equilibrium than in equilibrium at a given thermal energy $`T_{\mathrm{kin}}`$. It will thus be difficult to assign a unique effective temperature that reflects at the same time the reactivity of the molecules in the junction and the energy contained in the vibrations.

#### IV.3.4 Comparison between calculated and measured friction coefficients

The fit of curves equivalent to those shown in Figs. 4 and 5 allows one to estimate the kinetic friction force $`F_\mathrm{k}`$ with the help of Eq. (12). This result can then be compared to the friction force that is measured directly in the simulation. It turns out that such a comparison typically leads to an agreement within approximately 25 % accuracy, which can be improved by also taking into account the effects of instabilities on the motion normal to the surfaces. The deviation between the ‘predicted’ $`F_\mathrm{k}`$’s and the directly measured $`F_\mathrm{k}`$’s is due to the fact that the tails are not exactly exponential. This is particular important when the temperature is large or $`v_0`$ extremely small. If we accumulate the correct $`P(v)`$’s in the simulation and use Eq. (6) to predict $`F_\mathrm{k}`$, the agreement between predicted and observed kinetic friction is almost perfect, also when $`v_0`$ tends to zero.

Fig. 7 shows the degree of agreement for one particular model system. One can see that the kinetic friction coefficients $`\mu _\mathrm{k}`$ as obtained from the full velocity PD, see Eqs. (6) and (7) agree perfectly quite well with the directly measured $`\mu _\mathrm{k}`$. Neglecting the contribution of the motion normal to the surface results in an $`𝒪(20\%)`$ underestimation of the friction force. Estimating $`\mu _\mathrm{k}`$ indirectly with the help of Eqs. (11) and (12) leads to an underestimation of about 25%.

#### IV.3.5 Effect of temperature

It was shown by He and Robbins he01tl  that the model system on which this study is based yields logarithmic velocity corrections to the friction force $`F_\mathrm{k}`$ for incommensurate surfaces, provided the temperature is positive and the sliding velocity is not too small, see also the discussion in Ref. muser02prl . Our simulation results of the $`v_0`$ corrections to $`F_\mathrm{k}`$ for incommensurate surfaces are not shown explicitly, however, they confirm the results by He and Robbins, The basic reason for a logarithmic-type correction had already been recognized by Prandtl prandtl28 . Due to thermal fluctuations, the embedded atoms can jump over local energy barriers and the instabilities will be ignited prematurely. This reduces the necessary external force to maintain sliding, because it does not need to move the embedded atom all the way to the top of the energy barrier.

For commensurate surfaces, discontinuous instabilities are absent and therefore the effect of thermal fluctuations must be different. This issue is investigated in Fig. 8. Due to the large loads and the small temperatures employed, the linear-response regime is not necessarily reached at the sliding velocities $`v_0`$ accessible to the simulations, i.e., $`v_0=10^5`$. Therefore, we obtain kinetic friction coefficients $`\mu _\mathrm{k}`$ that apparently vanish according to

$$\mu _k\stackrel{v_00}{}v_0^\beta ,$$
(18)
with exponents $`0.25\beta 1`$.

It is remarkable that a small change in temperature has a rather strong effect on $`F_\mathrm{k}`$. For the small load $`L=1`$, the exponent $`\beta `$ is approximately unity at temperature $`T=10^2`$ and one may argue that the corresponding $`F_\mathrm{k}(v_0)`$ reflects a linear response curve. As $`T`$ is lowered to $`T=10^3`$ , a different exponent $`\beta `$ is obtained, reflecting non-equilibrium behavior. When the load is now increased by a factor of ten, the energy barriers also increase approximately by a factor of ten. Therefore the $`F_\mathrm{k}(v_0)`$ curves belonging to masses $`L10`$ should be considered far from equilibrium, i.e athermal. This would favour exponents $`\beta `$ less than unity. However, this expectation is not true. Instead a Stokes-type friction is observed. The almost linear relation of $`F_\mathrm{k}`$ and $`v_0`$ for these largests loads ($`L=100`$) may thus be an effect due to higher harmonics in the lubricant-wall potential. As one can see in the 1D, incommensurate systems, i.e., Fig. 2, the friction-velocity relationship can change qualitatively at certain critical values of the higher-order harmonics.

## V Beyond the impurity limit

So far, we have neglected the direct interactions between the impurities or the coverages were small enough in order to render the direct interactions negligible. This approximation is reasonable when the coverage is small and when the lubricant particles are simple spherical units without inner degrees of freedom. When either condition is violated, the energy landscape and hence the detailed characteristics of the instabilities will change. This in turn might lead to a qualitative change in the tribological behavior of the junction. In this section, we will study the applicability, the limitations and corrections of the impurity limit model that are due to the interactions between lubricant atoms.

### V.1 Coverage effects

When the lubricant coverage is close to or greater than one monolayer and the junction is sheared, particles will have to move in a correlated fashion. In order for one atom to jump to another mechanically stable site, its neighbor has to jump as well, etc. A detailed description of the dynamics will be very complicated, i.e., it may involve sliding of correlated blocks along grain boundaries and the formation of dislocation-type structures persson00jcp . Yet, the argument persists that instabilities and sliding induced deviations from the equilibrium velocity distribution function lead to friction.

Besides the correlated motion, some more details change when the coverage is increased. For example, pops also occur in the direction normal to the interface with a similar magnitude as parallel to the interface. This is reflected in the probability distributions $`P(v)`$ for the in-plane velocity $`v_{}=\sqrt{v_x^2+v_y^2}`$ and the normal component $`v_{}=v_z`$ of the fluid particles, see Fig. 9. The system under consideration is incommensurate, the walls are separated by a double layer and the externally imposed load per wall atom is $`L=30`$. Although the detailed dynamics of the lubricant atoms must be very different from those in the impurity limit, Eqs. (9) and (11) provide again a reasonable description for respectively $`P(v_{})`$ and $`P(v_{})`$, i.e., a central Maxwell-Boltzmann peak and a non-equilibrium exponential tail. Similar curves, which are not shown explicitly, were obtained for a coverage up to 5 monolayers.

As before, the kinetic friction force $`F_\mathrm{k}`$ is the integral over the deviation of the $`P(v)`$’s from the Maxwell-Boltzmann distribution, as stated in Eq. (6). $`F_\mathrm{k}`$ is shown for various coverages and sliding velocities in Fig. 10. Both commensurate and incommensurate systems are investigated and again their behavior is strikingly different.

We start our discussion with the commensurate system. At a coverage of $`C=0.25`$, results are very close to the impurity limit. $`F_\mathrm{k}`$ decays to zero with a powerlaw $`v^\beta `$ where the exponent $`\beta `$ is less than one. As the coverage is increased to $`C=0.5`$ or even $`C=0.75`$, $`F_\mathrm{k}`$ decreases considerably less quickly with decreasing $`v_0`$ than in the impurity limit. The behavior remains strikingly different from Coulomb friction. This changes when the coverage reaches and exceeds one full monolayer. For coverages beyond double layers, the kinetic friction force even exceeds that of incommensurate systems. The prediction in Ref. muser02prl,  that commensurate systems should show smaller kinetic friction than incommensurate system must thus be limited to extreme boundary lubrication. Above one monolayer lubrication, this trends seemingly turns around. Experiments suggest that commensurability leads to enhanced friction between mica surfaces lubricated by a double layer or more ruths00 . Unfortunately, no study is known to the authors in which a monolayer of lubrication or less was used between two (smoothly) sliding commensurate walls.

At the smallest velocity investigated, $`\mu _\mathrm{k}`$ increases by a factor greater than 200 for the commensurate case, when we increase the coverage from $`C=0.5`$ to $`C=2`$. The same change in coverage for incommensurate surfaces only yields a factor of two. Hence, incommensurate surfaces show much weaker coverage dependence than commensurate interfaces. Overall, there is relatively little change of $`F_\mathrm{k}`$ with coverage for incommensurate walls with the exception of $`C=1.5`$. Due to the large load employed, the 1.5 moloyers are squeezed into a single layer, which then essentially acts like a solid. This situation would not occur - or at least occur only for a short period of time - if the lubricant could flow out of the junction. We conclude that the coverage dependence is weak for incommensurate walls.

### V.2 Effects due to molecular bonds

Most lubricant particles possess an inner structure. Here we will focus on the most simple generalization of the spherical molecules considered so far, namely dimers, and hexamers (6-mers). Dimers would represent small linear molecules such as C<sub>2</sub>H<sub>6</sub>, while hexamers are representative of short, linear alkane chains. The dynamics of the lubricant particles will change due to the additional internal degrees of freedom. Alternatively, one may argue that the dynamics of monomers is restricted because every monomer is constraint by at least one chemically bonded neighbor.

While monomers only have translational degrees of freedom, dimers also have rotational degrees of freedom. It is tempting to speculate that ’rotational’ instabilities can occur in addition to the ’translational’ instabilities. Therefore, one might expect $`F_\mathrm{k}`$ to be larger for dimers than for monomers. However, the rotational and translation motion will not be independent of each other and the coupling between them might reduce the effect of a ’translational’ instability. The question which effect dominates can only be answered by analytical calculations or by molecular dynamics simulations. Simulation results for the kinetic friction force in a boundary-lubricated interface are shown in Fig. 11.

In Fig. 11, one can learn that the ratio $`\rho =d_{\mathrm{nn}}/d_{\mathrm{mol}}`$ of the next neighbor spacing $`d_{\mathrm{nn}}`$ in the walls and the intra-molecular bond length $`d_{\mathrm{mol}}`$ plays an important role, particularly for commensurate surfaces. When the intra-molecular bond length is close to the next-neighbor distance of wall atoms $`d_{\mathrm{nn}}`$, $`F_\mathrm{k}`$ disappears as a power law with sliding velocity $`v_0`$. This means that the ’interference’ effects between commensurate walls persist and that no instabilities occur. Surprisingly, this is even observed for hexamers. However, if $`d_{\mathrm{mol}}`$ differs from $`d_{\mathrm{nn}}`$, instabilities also occur in boundary-lubricated systems, even for commensurate walls. These instabilities are invoked through the rotational degrees of freedom. While the misfit between $`d_{\mathrm{mol}}`$ and $`d_{\mathrm{nn}}`$ leads to Coulomb friction between commensurate walls, its value of $`F_\mathrm{k}`$ remains small as compared to the incommensurate case.

We now turn to the incommensurate walls. Interestingly, the smoother walls with $`d_{\mathrm{nn}}=1.0`$ produce higher kinetic friction than the walls with $`d_{\mathrm{nn}}=1.2`$, while the opposite is true for obvious reasons in the case of static friction $`F_\mathrm{s}`$. The reason for the effect in $`F_\mathrm{k}`$ is that the reduced nearest neighbor spacing leads to a higher rate of popping processes, as more atoms sliding past each other at a given $`v_0`$ while the energy gain in the pops is only slightly decreasing. This effect can be verified by comparing the distributions of the particle velocities. The effect remains stable for all degrees of polymerization. The high friction of a dimer is caused by the contribution of their fast rotations. This is revealed by the distribution of the bonds’ angular velocities. Despite these trends, incommensurate systems prove again to be less susceptible to quantitative changes in the parameters that determine the details of the model than commensurate systems.

## VI Conclusions

Kinetic friction requires the prevailance of instabilities (mechanical hysteresis) in a system. In the present paper, we have focused on instabilities in the trajectories of particles confined between two walls, which are sheared against each other. When an instability is reached, the particle does not find a local potential energy minimum in its vicinity anymore and is thus forced to “pop” into the next local minimum it sees. At small sliding velocity $`v_0`$, this will lead to a high velocity, which depends solely on the energy landscape. The kinetic energy is gradually dissipated, resulting in a frictional force. We derived a relationship between the (non-equilibrium) velocity distribution function $`P(v)`$ and the friction force $`F_\mathrm{k}`$. The characteristics of $`P(v)`$ and thus $`F_\mathrm{k}`$ depend only weakly on coverage, sliding velocity, load and other parameters for incommensurate surfaces.

In a generic setup, we first used two Steele potentials reflecting two two-dimensional, triangular walls, which could be rotated with respect to each other to achieve an incommensurate system. We then computed numerically the adiabatic trajectory of a test particle. It was found that instabilities were a robust feature of the incommensurate system. Different off-symmetry wall rotations and inclusion of higher-order contributions to the Steele potential as well asymmetric interaction strengths of the walls did not alter the occurrence of the instabilities, but only affected their frequency.

Including interactions between lubricant atoms does not change the existence of instabilities and hence the presence of Coulomb friction either. In contrast, the commensurability of the walls allowed for especially smooth trajectories of impurity atoms. The trajectories remain smooth when interactions between lubricant atoms are included up to a coverage of one monolayer. Above one monolayer, the lubricant atoms do not move coherently any longer and instabilities are starting to occur within the film. Kinetic friction rises dramatically as a consequence.

We speculate that coherent motion similar to the one just described may also be responsible for the behavior observed in a pioneering quartz crystal oscillator study by Krim and Chiarello krim91JVTB . They found that the friction between a solid monolayer and a smooth surface was much smaller than the friction between a fluid monolayer and the same surface. The reverse was reported for a rough surface. Of course, the motion of a layer adsorbed on a microbalance is different from that of a ”between”-sorbed layer, because in the first case, there is no confining top wall and sliding-induced wall interference effects cannot occur. This is an important qualitative difference, which prevents us from making a direct comparison of our simulations with the above-mentioned experiments. One may yet argue that commensurability can induce coherent motion of the film, be it adsorbed or ”between”-sorbed. This suppresses erratic pops, which ultimately lead to energy dissipation. Thus, if one assumes that film and smooth substrate were commensurate, the small values for the kinetic friction force would not necessarily be in contradiction to the supposedly large static friction force. On the other hand for rough, disordered surfaces, a solid monolayer would not be able to move coherently, which would be consistent with its large friction as compared to a fluid layer.

We turn back to the discussion of the non-equilibrium velocity distributions. For incommensurate walls, the distribution consists of a central peak, which is essentially identical to the equilibrium velocity distribution, and of an additional non-equilibrium tails. These tails fall off only exponentially with $`v`$, which is slower than the exponential decay with $`v^2`$ in equilibrium systems. This observation is rather generic for incommensurate systems and independent of the lubricant coverage. As the real velocity distribution function is qualitatively different from Gaussians, it seems futile to describe the interface in terms of an effective temperature. We argued that given a specific kinetic energy associated with a lubricant (which could be used to define an effective temperature), the non-equilibrium system would be more likely to invoke chemical bond breaking or other chemical reactions.

Overall, the impurity model provides a good description of the typical characteristics of a boundary-lubricated system. However, it is essential to study two-dimensional interfaces and incommensurate surfaces. One-dimensional and/or commensurate surfaces lead to untypical behavior, i.e., rather large sensitivity of the friction force with respect to small changes in the model (details of interaction potential) or in the external parameters (sliding velocity, load, temperature, etc.). This is unfortunate, because incommensurate walls are much more common than commensurate walls, which leaves us with fewer possibilities to control friction.

A surprising result of our study for incommensurate walls is that increasing the atomic scale roughness of the walls may actually sometimes reduce the kinetic friction force.

It would be interesting to compare our predictions concerning the velocity distributions to experimental data. While scattering data from small, confined volumes is certainly notoriously difficult to obtain, recent advances have been made. Using fluorescence correlation spectroscopy, Mukhopadhyay et. al. mukhopadhyay02  measured translational diffusion in molecularly thin liquids confined within a surface forces apparatus. In the future, it might be possible to extend these studies to sliding situations so that velocity distributions can be measured.

###### Acknowledgements.

M.A. is indebted to J. Baschnagel and M. Fuchs for numerous discussions and acknowledges support by the German Academic Exchange Service (DAAD, “Hochschulsonderprogramm III von Bund und Ländern”), grant no. D/00/07994. We are also grateful to the Bundesministerium für Bildung und Forschung (BMBF) for partial support under grants no. D.I.P. 352-101 and 03N6500.

## Appendix A Relation between velocity distribution and friction

Consider a system in steady-state equilibrium with the following underlying equation of motion

$$m\ddot{𝐱}+m\gamma \dot{𝐱}=𝐅_\mathrm{b}(𝐱)+𝐅_\mathrm{t}(𝐱𝐯_0t)+𝚪(t).$$
(19)
Here, we chose the same terminology as in Sect. II, i.e., $`𝐅_\mathrm{b}(𝐱)`$ denotes the force of the bottom wall on an impurity atom located at position $`𝐱`$ and $`𝐯_0`$ is the velocity of the upper wall with respect to the substrate. We multiply Eq. (19) with $`\dot{𝐱}`$ and average over a long time interval $`\tau `$. We then interpret the resulting individual terms. They can be associated with the (average) power dissipated within the system or the (average) power put into the system. First, the average change of kinetic energy with time equals zero, namely

$`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑tm\ddot{𝐱}\dot{𝐱}`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑t{\displaystyle \frac{d}{dt}}T_{\mathrm{kin}}`$ (20)
$`=`$ $`{\displaystyle \frac{1}{\tau }}[T_{\mathrm{kin}}(t=\tau )T_{\mathrm{kin}}(t=0)]`$
$``$ $`0\mathrm{for}\tau \mathrm{}`$
The second term is proportional to the time-averaged kinetic energy of the system with respect to the lower wall:

$`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑tm\gamma \dot{𝐱}\dot{𝐱}`$ $`=`$ $`\gamma m\dot{𝐱}^2`$ (21)
$`=`$ $`2\gamma T_{\mathrm{kin}},`$
$`T_{\mathrm{kin}}`$ being the time-averaged or ensemble-averaged (steady-state) kinetic energy of an impurity. Thermostating also parallel to the top wall, e.g. by choosing $`\gamma =\gamma _\mathrm{t}=\gamma _\mathrm{b}`$, requires a trivial modification of the reference system. The next term is the average work of the bottom wall on the impurity

$`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑t\dot{𝐱}𝐅_\mathrm{b}(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑t\left({\displaystyle \frac{d}{dt}}V_\mathrm{b}(𝐱)\right)`$ (22)
$`=`$ $`{\displaystyle \frac{1}{\tau }}\left\{V_\mathrm{b}[𝐱(\tau )]V_\mathrm{b}[𝐱(0)]\right\}`$
$``$ $`0,`$
which must vanish in any steady-state system. Of course, if the model was generalized such that (steady state) wear would occur, then the contribution discussed in Eq. (22) would indeed remain finite.

For the discussion of the next term in Eq. (19), it is necessary to keep in mind that

$$\frac{d}{dt}V_\mathrm{t}(𝐱v_0t)=F_\mathrm{t}(𝐱v_0t)\left[\dot{𝐱}v_0\right].$$
(23)
This and the same considerations invoked for Eq. (22) yield

$$\frac{1}{\tau }_0^\tau 𝑑t\dot{𝐱}𝐅_\mathrm{t}(𝐱𝐯_0t)=𝐅_\mathrm{t}𝐯_0,$$
(24)
where $`𝐅_\mathrm{t}`$ is the time- or ensemble-averaged force that the top wall exerts on an impurity. This force or depending on the definition its projection onto the sliding direction can be associated with the kinetic friction force $`F_\mathrm{k}`$.

The contribution due to the random force $`𝚪(t)`$ is the most difficult contribution to calculate. However, if the system is close to local equilibrium for most of the time, then the expectation value of $`𝚪(t)\dot{𝐱}`$ can be expected to be close to the value of this expression in thermal equilibrium. In equilibrium, it must compensate the expression discussed in Eq. (21), hence

$$\frac{1}{\tau }_0^\tau 𝑑t\dot{𝐱}𝚪(t)2\gamma T_{\mathrm{kin}}_{\mathrm{eq}},$$
(25)
where $`T_{\mathrm{kin}}_{\mathrm{eq}}`$ denotes the average kinetic energy in thermal equilibrium.

Assembling all necessary terms, yields

$$2\gamma (T_{\mathrm{kin}}T_{\mathrm{kin}}_{\mathrm{eq}})=F_\mathrm{k}v_0.$$
(26)
Note that Eq. (26) is equivalent to Eqs. (6) and (7).


format: text

text:
# Addendum to “Coherent radio pulses from GEANT generated electromagnetic showers in ice”

## Abstract

We reevaluate our published calculations of electromagnetic showers generated by GEANT 3.21 and the radio frequency pulses they produce in ice. We are prompted by a recent report showing that GEANT 3.21-modeled showers are sensitive to internal settings in the electron tracking subroutine. We report the shower and pulse characteristics obtained with different settings of GEANT 3.21 and with GEANT 4. The default setting of electron tracking in GEANT 3.21 we used in previous work speeds up the shower simulation at the cost of information near the end of the tracks. We find that settings tracking $`e^{}`$ and $`e^+`$ to lower energy yield a more accurate calculation, a more intense shower, and proportionately stronger radio pulses at low frequencies. At high frequencies the relation between shower tracking algorithm and pulse spectrum is more complex. We obtain radial distributions of shower particles and phase distributions of pulses from 100 GeV showers that are consistent with our published results.

In a recent paper soeb  we reported the results of an extensive study of the radio frequency pulse emitted by an electron-induced electromagnetic shower. Application to radio signal detection of ultrahigh energy electron neutrino interaction in ice provided our motivation radhep . We used the standard particle physics detector-simulation package GEANT 3.21 (G3 afterwards) geantman  with default settings to generate the showers. Our total track-lengths of electrons and number of electrons at shower maximum were 68% of those in the paper by Zas, Halzen and Stanev (ZHS) zhs . Our pulse height at the Cherenkov angle at 1 GHz was 72% of the ZHS result. Our study of sensitivity of shower development to reasonable changes in input cross sections showed no effects large enough to explain the difference between G3 and ZHS shower profiles and track lengths.

The source of the bulk of the discrepancy was recently presented by Alvarez-Mũniz, Marquéz, Vázquez and Zas jaime . They noticed that the default setting of a parameter called “ABAN” in the electron tracking subroutine of G3, the setting used in our studies, gives shorter total, projected and weighted track lengths and lower particle yield at shower maximum than GEANT 4 (G4 afterwards) or G3 with a choice of tracking rules different from the default. Their results for track-length and shower profile from the latter simulations are 90% of those of the original ZHS work for these quantities and in closer agreement than G3 showers with the default setting soeb ; jaime .

Track lengths are one indicator of the strength of electric field pulses to be expected from the shower. However the electric field at a given angle and frequency depends on a combination of many factors. We make a comparison among the shower and electric field pulse results from G3 (with several different combinations of electron tracking rules) and G4. We report results for shower profiles, total track lengths, energy loss ($`dE/dx`$) and radial distribution of the number of electrons minus the number of positrons around the shower maximum. We tabulate the G3, G4 and ZHS track-lengths and yields at shower maximum for ease of comparison of shower results. Regarding the electric field in the Fraunhoffer regime, we report radio pulse strength at the Cherenkov angle and the phase behavior of the radio pulse both on and off the Cherenkov peak. We find that the particle yield and charge imbalance are larger with improved tracking. At low frequencies the pulses emitted at the Cherenkov angle are correspondingly stronger as well.

Electron tracking in GEANT: The subroutine responsible for tracking electrons in G3 is called GTELEC. In its default setting, ABAN=1, if the electron range is too short to reach the boundary of a detector volume, and the distance to its next bremsstrahlung emission is greater than its range, the electron’s energy is deposited in the volume and it is stopped. This is done regardless of the value of the electron kinetic energy threshold parameter CUTELE set by a user in the input card file. This was presumably made the default to increase the calculation Al speed of high energy particle physics detector simulations in non-sensitive volumes, where detailed tracking was not needed aban . With the default setting, however, track lengths are reduced significantly, with a consequent reduction in the calculated radio Cherenkov emission. This is what we observed in our previous study soeb , where we compared results from G3 with default setting to the results of zhs .

In practice, it seems for applications such as ours the preferred michel1  setting is ABAN=0, which tracks all electrons down to the threshold set by the user. The setting ABAN=2 allows the user to track electrons precisely (similar to ABAN=0 setting) in sensitive volumes or not so precisely in the non-sensitive volumes isvol . “Sensitive” vs. “non-sensitive” settings are distinguished by the value of a step parameter called “STMIN” in the tracking medium parameter list, where “sensitive” produces many more steps per track than “non-sensitive”. G4 has range cuts that correspond to energy thresholds for the production of particles. A particle, once produced, is tracked down to rest regardless of its production threshold.

Results: We expand on the study reported in Ref. jaime , where showers produced with G3 ABAN=1, G3 ABAN=2 with “sensitive” volume, G4 and ZHS are compared. We run G3 for four different tracking choices, namely ABAN=0, ABAN=1, ABAN=2 sensitive volume (S) and ABAN=2 non-sensitive volume (NS) settings, to create showers in ice. We run G4 for comparison as well. We then repeat our study soeb  of the key features of the showers and the consequent key features of the radio pulses: pulse angular distributions, phase distributions and frequency spectra.

1) Shower Characteristics: First we investigate the important issue of total track lengths in Table 1. We note that the track lengths for ABAN=2 (NS) setting are similar to the default ABAN=1 setting; the ABAN=2 (S) setting used in Ref. jaime  gives results similar to the preferred ABAN=0 setting and both are similar to the results from G4. The agreement between track lengths produced by these latter three GEANT versions and those produced by the ZHS code improves as one compares total, total projected $`e+p`$ and total projected $`ep`$ track lengths. Comparing the projected $`e+p`$ and $`ep`$ track lengths, we see that the loss of track length for ABAN=1 is more significant for electrons than for positrons. We expect this to be the case, since a large part of the electron excess comes from low energy electrons swept into the shower by Compton scattering and $`\delta `$-ray production.

Despite the improved agreement on total track lengths, there are still significant differences between the G3 ABAN=0, G3 ABAN=2(S), and G4 simulations. For each code individually, however, the variance in total track length $`3640`$ cm is reasonable. For a given shower energy $`E_0`$, the total track length $`L_{tot}`$ is determined primarily by the rate of ionization energy losses, where the energy transfer is less than a low threshold $`E_t100`$ keV. Then $`L_{tot}E_0/(dE/dx)_{ion}`$, as expected under the conditions that a) soft ionization losses are essentially independent of particle energy and b) hard interactions do not “remove” energy from the shower. Hard $`\delta `$-rays create their own track length and hard $`\gamma `$’s return their energy to the shower in pair production and Compton scattering processes. Some energy is lost from the shower as soft particles fall below $`E_t`$. The fluctuations in track length are thus determined by the fluctuations in the number of soft particles multiplied by $`E_t/(dE/dx)_{ion}`$. Estimating the number of soft particles by $`N_t=E_0/E_t`$, the fluctuations in total track length are of order $`\delta L_{tot}E_t/(dE/dx)_{ion}\sqrt{N_t}=\sqrt{E_oE_t}/(dE/dx)_{ion}=50`$ cm, for $`E_0=100`$ GeV, $`E_t=100`$ keV, and $`(dE/dx)_{ion}=2`$ MeV/cm. Although this is a rough estimate, the observed values of $`\sigma `$ for $`L_{tot}`$ are consistent with expectations for the newer simulations.

This estimate leads to several observations. a) The deviations in $`L_{tot}`$ between G3 ABAN=0, G3 ABAN=2(S), and G4 are small but significant; this is under study. b) Due to the default decision in the tracking code to abandon some particles, the older G3 simulations have an effective value of $`E_t5`$ MeV, which implies $`\delta L3.5`$ m, generally consistent with the Table entries. c) For $`ep`$ track length, the fluctuations are dominated by the competition between pair production and Compton scattering, for which one may estimate $`E_t20`$ MeV and $`\delta L7`$ m.

Particle yields at the maxima of longitudinal shower profiles are also listed in Table 1. The results from the ZHS code, G3 with ABAN=0, ABAN=2 with sensitive volume, and G4 all agree reasonably well. Agreement is especially good for the excess charge, where results agree within errors. Again, the low yields from G3 ABAN=1 and ABAN=2 (NS) are consistent with the interpretation that lost track length is caused by a deficiency of low energy electrons being scattered into the shower.

Next we consider how individual particle tracks contribute to the radio pulse. For low frequencies it can be shown that the contributions are in phase and the pulse is proportional to the weighted, projected track length, labeled (e-p) in Table 1. As the frequency increases, this simple coherence breaks down as described by two parameters, which we called soeb  the Cherenkov Phase (CP) and the Translational Phase (TP). The CP is related to the phase coherence within a single track segment, whereas the TP describes the relative phase among track segments due to their different positions within the shower, relative to the observer.

In general, we find that the different simulations produce phase distributions similar to those shown in Fig. 17 of Ref. soeb . Specifically, the radial charge distribution at shower maximum, shown to peak at $`\rho 0.5`$ cm in Fig. 22 of Ref. soeb , is reproduced apart from an overall increase in normalization buniy . Moreover, the radial charge distribution early in the shower evolution has a core dominated by high energy particles, and the normalization of this core does not depend on the details of the electron tracking at low energy.

2) Radio frequency pulse: Next we examine how the changes in shower simulation affect the modeled radio Cherenkov pulses. At low frequencies, all particles contribute in phase, so the increase in the excess track length found in G3 ABAN=0 and G4 should be reflected by an increase in the amplitudes of the low frequency part of the spectrum as compared to showers produced with G3 ABAN=1. In Fig. 1 we see that this is indeed the case, where the spectrum $`R\times |\stackrel{}{E}|`$ is shown for an observer on the Cherenkov cone.

In Fig.2 we extend the frequency scale out to 10 GHz, which suffices to show clearly a transition from low to high frequency for the same sets of simulations. In this range, the spectra for the G3 ABAN=0 and G4 showers flatten out, but the G3 ABAN=1 spectrum grows slowly, so that the difference in their normalizations, relative to G4, is less than 20% for $`f=10`$ GHz. We infer that in this region the spectra are not simply described by the phase coherent sum of amplitudes from a large number of low energy electrons. We note two effects that could be important in this region. First, a relatively small number of high energy particles near the shower axis could contribute strongly to the spectrum in this region. Since this part of the shower evolution is not sensitive to low energy electron tracking, we may see the G3 and G4 results converging. This trend continues at higher frequency. Second, a large fraction of the track-segment amplitudes are no longer coherent. For that fraction, the magnitude of $`\stackrel{}{E}`$ increases as a random walk as the track length increases. In this case, one would expect the amplitude to increase as the square root of the number of track segments, rather than linearly, which reduces the ratio of the pulse from the G4 generated shower to that from ABAN=1 as the frequency increases. This is seen clearly in Figs. 1 and 2. Given that total and projected track lengths are proportional to energy soeb ; zhs , the effects just outlined should be relevant when comparing pulses from showers of different energies using the same simulation.

Figure 3 shows a comparison between pulse spectra from 1 TeV showers and 100 GeV showers, both generated by G4. The spectra are extended up to 50 GHz. For frequencies above 15 GHz, the results are fully consistent with Fig. 25 of Ref. soeb . At low frequencies, below a few GHz, the results are consistent as well, after accounting for the change in normalization between G4 and G3 ABAN=1. At 50 GHz, the TeV average, scaled down by a factor of 10, is smaller by roughly a factor $`3\sqrt{10}`$ than the 100 GeV average, possibly hinting that the “random walk” increase of $`\stackrel{}{E}`$ with energy, and therefore total track length, is setting in. At the same time, it seems clear that a combination of effects must be in play for intermediate frequencies.

Since the flat spectra observed above several GHz contrast with the sharp falloff expected above the 1-2 GHz scale zhs ; avz00 , we referred to this phenomenon as an “extended coherence regime” in . This is different from the use of coherence to refer to the rapid rise of the electric field strength as frequency increases in the region where frequency is low compared to the scale set by the peak at 0.5 cm in the shower radial distribution.

Summarizing, we see that the improved electron tracking performed by G4 or G3 ABAN=0 results in proportionately increased radio Cherenkov spectra at frequencies below a few GHz. At high frequencies the situation is more complicated. We are conducting further studies at a variety of energies and including showers initiated by $`\gamma `$-rays to identify the determining factors.

Summary and further discussion of results:

Motivated by the observation in jaime  that the default setting for electron tracking in G3 stops many electrons prematurely, we have reviewed the calculations reported insoeb . In soeb  we used G3 with default setting and found a persistently smaller shower and a weaker electromagnetic pulse at low frequencies than those reported in zhs . We find that this discrepancy largely disappears when electrons are tracked more completely in G3 or when using G4 generated showers, similar to results reported in jaime . Figure 1 shows that the strength of the pulses at low frequency clearly reflect the track lengths and particle yield at shower maximum. The high frequency behavior is more complex and is under further study.

The particles in the electromagnetic showers in a uniform medium, ice in the case of interest here, should be tracked as long as they still contribute to the Cherenkov signal. This is the case, for the ZHS, G4 and G3 simulations with the tracking option ABAN=0 or ABAN=2 with sensitive volume; the results from all these simulations agree at the 10% level. There is physics, namely significant emission, that is lost if the tracks are stopped prematurely and their energy dumped into the medium. Though choosing the default (ABAN=1) setting affects only the low energy tail of the shower, below the few MeV range, the population of particles (predominately electrons) is so high at the low energy end that the field contribution cannot be ignored.

For applications of the shower and pulse codes to neutrino detection experiments, the difference is obviously important. For example, below 1 GHz where currently deployed antennas are sensitive, the expected radio emission pulse is 1.3 to 1.4 times higher using the G3 or G4 code with particles tracked to threshold compared to the truncated showers produced by G3 with the ABAN=1 setting. The more correct treatments that follow the particles to lower energy produce significantly larger fields on the Cherenkov cone. This increases the effective volume estimates and leads to higher sensitivity estimates for cosmic neutrino detection systems like the Radio Ice Cherenkov Experiment (RICE) ilya  and the recently approved balloon-borne experiment ANITA anita , both using the Antarctic polar ice cap as a neutrino target.

Acknowledgements: We thank J. Le Vaillant and J. Adams for help running GEANT simulations. We acknowledge helpful communications with M. Maire, J. Alvarez-Muñiz, E. Marqués and E. Zas. Work was supported in part by the NSF, the DOE and the University of Kansas General Research Fund. S. Seunarine was supported by a Marsden grant. S. Razzaque was partially supported by NSF AST0098416 grant.