text:
# Resonant Scattering of Emission Lines in Coronal Loops: Effects on Image Morphology and Line Ratios

## 1 Introduction

Two effects of resonant line scattering on emission lines are to change the morphology of images and to change line ratios from what would be expected from optically thin analyses. Qualitatively, the effects can be split into three broad optical depth categories, low, moderate, and large. At low optical depths ($`\tau <5`$) resonant scattering primarily changes the photon direction and photons emerge into the directions of lowest optical depth. At moderate optical depths ($`\tau 10100`$) multiple scattering is important and there is a large probability of photon energy conversion, e.g., the conversion of Ly$`\beta `$ photons to H$`\alpha `$ photons. Large optical depths ($`\tau >>100`$) result in photon trapping and collisional de-excitation, eventually leading to thermalization. The goal of this paper is to investigate the effects of resonant scattering on images and intensity ratios in the low optical depth regime. Situations where our work is applicable include UV, EUV, and soft X-ray emission line imaging and spectroscopy of the Sun and stellar coronae. While we have chosen a geometry based on coronal magnetic loops, similar considerations would apply, for instance, to some of the strong X-ray emission lines in older supernova remnants. We have simulated an X-ray line of Fe xvii, but the predicted effects on image morphology should be visible in the narrow-band EUV images of EIT and TRACE.

Ultraviolet and soft X-ray emission line ratios are frequently used as diagnostics of temperature, density, elemental abundances, and the ionization state of astrophysical plasmas. Several assumptions are often implicitly made when emission line ratios are analyzed. In particular, the optical depth is generally assumed to be negligible. However, solar coronal observations of the Fe xvii line at 15.01Å show line ratios that are lower than theoretical predictions (Phillips et al. 1997; Schmelz et al. 1997; Saba et al. 1999; Waljeski et al. 1994) or laboratory measurements (Brown et al. 1998). Rugge & MacKenzie (1985) suggested that optically thick resonant scattering of Fe xvii 15.01Å photons out of the line of sight may be responsible for the small line ratios. Saba et al. (1999) employed the escape probability method of Kastner & Kastner (1990) in their analysis of various line ratios and found that up to 50% of the Fe xvii flux is unaccounted for. Bhatia & Kastner (1999) found similar results for planar and cylindrical geometries. This explanation for the intensity decrements is likely, given the large oscillator strength of the Fe xvii 15.01Å line (Table 1). For typical densities ($`n_e10^{10}`$) and sizescales ($`10^9`$cm) in the solar corona, the optical depth of the Fe xvii 15.01Å line is around 1. This is in the low optical depth regime in which the major effect is scattering of photons into optically thin directions.

In this paper we present three dimensional Monte Carlo radiation transfer simulations of the transfer of Fe xvii 15.01Å photons. Our analysis does not treat transfer within the line itself, instead considering all photons to be emitted and resonantly scattered at a single wavelength. This approximation is valid so long as the optical depth is not large enough that photons scatter out of the Doppler core. The advantage of our Monte Carlo technique is that it naturally accounts for arbitrary illumination and multiple scattering in complex geometries. The results of our Monte Carlo simulations are images that may be compared directly with optically thin simulations to determine differences in image morphology and optically thick to optically thin intensity ratios.

Previous simulations by Alexander & Katsev (1996) produced optically thin images of loops in which the density and temperature follow the scaling law prescription of Rosner, Tucker, & Vaiana (1978, hereafter RTV). The simulations of Alexander & Katsev showed that it is crucial to consider the geometry and temperature dependence of the emissivity along with the thermal response of the imaging telescope when analyzing solar emission line images. Such images may be in the light of a single line or an ensemble of lines as is the case with the broad wavelength response of Yohkoh SXT. In particular, a loop viewed face-on may appear uniformly bright, but line-of-sight effects can yield an apparent “looptop” source when the same loop is viewed close to edge-on (see Fig. 5 of Alexander & Katsev). Our work is a natural extension of the Alexander & Katsev analysis to optically thick, multiple scattering situations, and for continuity purposes we also adopt RTV loops for our simulations.

Section 2 presents the model ingredients (loop structure, Fe xvii emissivity and opacity, and radiation transfer technique). Results of our simulations (images and intensity ratios of optically thick/optically thin models) are presented in § 3, and in § 4 we discuss our results with respect to observational data.

## 2 Model Ingredients

In order to construct model images we must specify the emissivity, opacity, and scattering geometry within our Monte Carlo radiation transfer code. These model ingredients are now described.

### 2.1 Loop Structure

Following Alexander & Katsev, we consider loops with a density and temperature structure given by the RTV scaling law model. The RTV model is a hydrostatic loop model in which the maximum loop temperature, $`T_{\mathrm{max}}`$, and loop length, $`L`$, are related to the pressure, $`p`$, by $`T_{\mathrm{max}}1.4\times 10^3(pL)^{1/3}`$. The maximum temperature occurs at the loop apex, and the balance among thermal conduction, radiative cooling and mechanical heating determines the temperature distribution along the length of the loop. An RTV loop is therefore specified by two parameters, the maximum temperature, $`T_{\mathrm{max}}`$, and loop length, $`L`$, measured along the loop from a footpoint to the looptop. Throughout our simulations we assume the loops are semicircular tori with length $`L=10^{10}`$cm and radius $`r=L/10`$. We then construct simulations for a range of maximum temperatures and viewing angles. See Fig. 1 for a diagram showing the loop length, radius, and orientation.

In some cases we embed the emitting loop in a lower density background plasma. Alexander & Katsev modeled isolated loop structures. We wish to investigate situations where a background plasma may provide additional scattering opacity, thereby allowing for attenuation of loop emission and the formation of scattering halos around the loops. For the background we embed the RTV loop in a hemispherical isothermal constant density plasma. We adopt a temperature of $`T=2.5\times 10^6`$K and set the density such that the total emissivity from the hemispherical volume is half of that emanating from the RTV loop (see § 3).

We note here that an isothermal, constant density background plasma may be physically unrealistic. The precise form of the background is unimportant for the purposes of this work, as we wish to show the general effects of resonant scattering. For future modeling of images from TRACE, EIT, and SXT where resonant scattering effects are important, more realistic loop structures and backgrounds will be generated.

### 2.2 Fe xvii Emissivity and Opacity

The Fe xvii 15.01 Å line emissivity and opacity are functions of temperature within the loop. The Fe abundance was taken as the typical FIP (first ionization potential) enhanced value of three times the photospheric abundance ($`1.2\times 10^4`$; e.g. Meyer 1985; Feldman 1992). We have taken the Fe xvii ionic fraction as a function of temperature from the paper of Arnaud & Raymond (1992), which is similar to the more recent results of Mazzotta et al. (1998). Because of its stable Ne-like electron configuration, Fe xvii is the dominant ion over a broad temperature range, and its emission lines are among the brightest observed in astrophysical plasmas at temperatures around $`3\times 10^6`$ K. The excitation rate for the 15.01 Å line ($`2p^62p^53d^1P_1`$, Table 1) is taken from Smith et al. (1985), and it is very close to the results of other computations (e.g., Goldstein et al. 1989). The oscillator strength is from Bhatia & Doschek (1992). To obtain optical depths, we assume a line width of $`30\mathrm{km}\mathrm{s}^1`$, typical of the solar corona (Mariska, Feldman & Doschek 1979). The emissivity and opacity for the range of temperatures we consider are displayed in Fig. 2.

### 2.3 Radiation Transfer

We perform the radiation transfer with a Monte Carlo radiation transfer code that accounts for arbitrary sources of emission and multiple scattering within an arbitrary geometry. The scattering code is based on that described by Code & Whitney (1995), but has now been modified to run on a three dimensional linear Cartesian grid (Wood & Reynolds 1999) and includes forced first scattering (Witt 1977) and a “peeling off” procedure (Yusef-Zadeh, Morris, & White 1984). These modifications enable us to construct model images of three dimensional systems from specified viewing angles very efficiently.

For the current investigation we set up our grid with an RTV loop structure. The temperature dependent emissivity and opacity within the loop are determined as described above. We assume that all Fe xvii photons are emitted and resonantly scattered (with an isotropic scattering phase function) at a single wavelength. It is possible that when a 15.01 Å photon is absorbed, the excited level will decay to a 3p level, converting the energy into two EUV photons and a 3s-2p photon near 17 Å. However, the branching ratio for this conversion is only $`5.7\times 10^4`$, which can safely be ignored for the optical depths less than about 10 that we consider. We also assume that the photon does not leave the Doppler core. This is a good approximation for modest optical depths, though it would seriously underestimate the flux emerging in the line wings for large optical depths. We use the average scattering cross section over a Gaussian profile, rather than using the actual cross section as a function of velocity for computational speed.

## 3 Results

For each simulation of an RTV loop ($`L=10^{10}`$cm) we form images at different viewing angles. In § 3.1 we show some representative Monte Carlo images that illustrate the effects of optical depth and viewing angle. From these images we compute the total emergent intensity for the particular viewing angle for comparison with an optically thin simulation (§ 3.2).

### 3.1 Images

In Figs. 3, 4, and 5 we show a set of model images for loops of different $`T_{\mathrm{max}}`$ at different viewing angles. We show optically thin images (assuming zero optical depth), and optically thick images with and without background plasma. In this and the following section we refer to viewing angles $`\theta `$ and $`\varphi `$, where the polar angle, $`\theta `$, is measured from the $`z`$ axis and the azimuthal angle, $`\varphi `$, is measured anti-clockwise from the $`x`$ axis. In a Cartesian coordinate system, the loop is oriented in the $`zy`$ plane (see Fig. 1). There are several features of the images that we now describe.

The $`T_{\mathrm{max}}`$ dependence affects the distribution of the emissivity along the loop. For the lowest temperature, $`\mathrm{log}T_{\mathrm{max}}=6.4`$ (Fig. 3), the emissivity is concentrated towards the loop apex. This is because at $`\mathrm{log}T_{\mathrm{max}}=6.4`$ the Fe xvii 15.01Å line emissivity is almost zero (Fig. 2). This temperature occurs at the loop apex and the loop legs are cooler resulting in the concentration of emission at the apex and negligible emission from the loop legs. The emissivity distribution along the loop changes as we increase $`T_{\mathrm{max}}`$. As $`T_{\mathrm{max}}`$ increases ($`\mathrm{log}T_{\mathrm{max}}=6.6`$ in Fig. 4 and $`\mathrm{log}T_{\mathrm{max}}=6.8`$ in Fig. 5) the emissivity becomes more concentrated towards the legs and footpoints of the loop. For $`\mathrm{log}T_{\mathrm{max}}=6.8`$, the Fe xvii 15.01Å line emissivity is again almost zero (Fig. 2) and in this case the emission is concentrated towards the loop legs and there is negligible emission from the loop apex.

The geometrical line-of-sight effects for the optically thin images discussed by Alexander & Katsev are clearly seen in column A of Figs. 3, 4, and 5. When viewed from above, geometrical projection makes the footpoints appear bright (image A1 in Figs. 4 and 5), while the apex of the loop is brightened by projection if the loop is viewed edge-on (image A5 in Figs. 4 and 5). Loops that appear uniformly bright when seen face-on (image A3 in Figs. 3 and 4) appear to have bright looptops when viewed edge-on (image A5 Figs. 3 and 4). Similarly, loops that show brightening along the legs for face-on viewing (image A3 in Fig. 5), appear more uniformly illuminated when viewed edge-on (image A5 in Fig. 5). This is simply due to the increased path length towards the apex for edge-on loops.

Optical depth drastically reduces the geometrical projection effects because the long path length that enhances the brightness of an optically thin line implies a large optical depth for an optically thick line. The effect on image morphology is that isolated loops appear more uniformly bright (compare columns A and B in Figs. 3, 4, and 5). This arises because we see emission originating from the “optical depth one surface” which covers a progressively larger projected area on the loop as it becomes optically thicker. This effect is shown quantitatively in Fig. 6 where we show the intensity variation across two of the loop images from Fig. 4. Figure 6a shows the intensity variation across the images in Fig. 4 for the overhead view (images A1 and B1 in Fig. 4). The solid line shows the intensity variation across the optically thin image (image A1 in Fig. 4) with the intensity peak towards the legs of the loop evident. The dotted line shows the intensity variation for the same loop, but with the effects of resonant scattering included (image B1 in Fig. 4) — the overall intensity level is lower and the intensity variation is more uniform. Figure 6b shows the intensity variation across the images in Fig. 4 for the edge-on view (images A5 and B5 in Fig. 4). Again the solid line shows the intensity variation along the optically thin image (image A5 in Fig. 4) showing the intensity peak at the loop apex. The dotted line shows a much reduced intensity level and more uniform intensity for the loop model that includes resonant scattering effects (image B5 in Fig. 4).

Real solar active regions are complex structures, with loops of differing $`T_{\mathrm{max}}`$, length and pressure. The emissivity scales as $`n_e^2`$, while the opacity scales as $`n_e`$, so that dense, low-lying loops may dominate the emission, but a background plasma may provide additional Fe xvii scattering opacity. We simulate this effect by embedding the emitting loop in a lower density background (column C in Figs. 3, 4, and 5). The total emission from the background is half that from the RTV loop but the background is very diffuse (see Table 2), so the RTV loops are clearly visible.

Two main effects of the background opacity are easily seen in Fig. 5. The first is the scattered light halo present around the loops (e.g., image C5 in Fig. 5), arising from scattering of loop emission in the background plasma. Such halos might appear as fuzziness in EIT or TRACE images in EUV emission lines. The second effect is that the loop emission must pass through the background plasma and when the optical depth is large this results in scattering out of the line of sight and an asymmetry between the near and far side of the loops, resulting in an apparent one-sided loop (e.g., image C4 in Fig. 5). Asymmetric loops are often observed in SOHO/EIT and TRACE images (examples may be seen on the TRACE homepage: http://vestige.lmsal.com/TRACE/Public/Gallery/Images/). Most are probably asymmetric in reality, but optical depth may play a role in some observed morphologies.

We note that the 171 and 195 Å bands of EIT and TRACE and the 284 Å band of EIT are dominated by resonance lines of Fe ix, Fe x, Fe xii and Fe xv. While these ions do not occupy as large a temperature range as Fe xvii, the scattering cross section scales in proportion to wavelength, and images in these lines are likely to show similar effects to those predicted for Fe xvii (it is the line optical depth and not the wavelength that determines the image morphology in our simulations). Indeed, Schrijver & McMullen (2000) have inferred that scattering effects are important from the variation of intensity with height above the limb from TRACE and EIT images.

### 3.2 Unresolved Intensity

In the previous section we showed spatially resolved model images in the light of the Fe xvii line, as may be observed in the solar corona. Observations of stellar coronae cannot spatially resolve coronal loops and the CHANDRA satellite has recently observed Capella and other active late-type stars as part of the Emission Line Project. As our resolved images showed, optically thick resonant scattering yields image morphologies and intensity distributions along the loops that differ from optically thin simulations. In this section we investigate the effects of optically thick resonant scattering on the Fe xvii line intensity for different orientations of unresolved loops, and compare it with the intensity expected if the line were optically thin. For optically thin emission all emitted photons escape from the source isotropically, and all loops will yield the same intensity independent of viewing angle. Optical depth effects result in photons exiting the simulation region preferentially along optically thinner sightlines. Also, for optically thick loops, the intensity is proportional to the projected area of the loop, so face-on loops will yield larger intensities than edge-on loops.

Figure 7 shows the intensity ratio for an unresolved optically thick model compared to the corresponding model when optical depth effects are ignored. Each panel shows the intensity ratio for a range of $`T_{\mathrm{max}}`$, with the individual curves showing different viewing angles of the loop.

The left column shows intensity ratios for an isolated RTV loop (no background plasma). Face-on viewing ($`\varphi =0`$, $`\theta =\pi /2`$) results in intensities that are larger than the optically thin case, due to photons being scattered into our line of sight. For edge-on viewing ($`\varphi =\pi /2`$, $`\theta =\pi /2`$), the intensity is lower than the optically thin case. This is for two reasons; firstly the projected emitting area is smaller for edge-on loops (see the loop images in Figs. 3, 4, and 5). Secondly as the optical depth increases, photons from the far side of the loop have to traverse a greater optical depth through the loop and are more likely to be scattered into the optically thinner face-on directions. Depending on $`T_{\mathrm{max}}`$ and loop orientation, the intensity can range from being 20% brighter to 50% fainter than the intensity from an optically thin loop.

The second column in Fig. 7 shows intensity ratios for loops embedded within a background plasma. The background has a temperature of $`T_{\mathrm{back}}=2.5\times 10^6`$K and the density is chosen so that the total intensity from the background is half that from the loop (Table 2). The optically thick to optically thin intensity ratios for these simulations show that the intensity ranges from around 60% to 120% of the optically thin intensity depending on $`T_{\mathrm{max}}`$ and viewing angle. These simulations show that optical depth, geometry, and loop orientation can have a large effect on observed line ratios. Clearly great care is required in interpreting coronal images and line ratios due to the effects of projection and optically thick resonant scattering.

## 4 Summary

Our simulations have shown that optical depth effects, both within loops and in a background plasma, can change the amount of intensity emerging into a given direction and also the image morphology. Intensity ratios may vary in the range 50% to more than 120% of their optically thin values. If one observes an ensemble of loops, for instance spread over the surface of another star, the effects tend to average out. Equal numbers of photons are scattered into and out of the line of sight. However, a net effect may still arise if the emissivity and opacity vary in different ways. Given the electron density factor in the emissivity and its absence in the opacity, such differences are likely. To account for the general trend of reduced Fe xvii 15.01 Å intensity (e.g. Saba et al. 1999) by scattering rather than errors in atomic rates, the typical structure must preferentially scatter photons down to the solar surface. Such an example may be emission from dense, low-lying loops that is scattered in a more diffuse, larger scale region. This is expected from RTV models and is observed in solar active regions. Models similar to those in the righthand column of Fig. 5 can easily be constructed to reduce the intensities of lines having large oscillator strengths. They may be invoked when line ratio anomalies appear in Chandra or XMM data. Proving that optical depth is indeed the explanation for a given anomaly will undoubtably be a challenge.

KW acknowledges support from NASA’s Long Term Space Astrophysics Research Program (NAG5-6039), and JR acknowledges support from NASA grant NAG5-2845. We thank David Alexander, Carolus Schrijver, Rebecca McMullen, Julia Saba, and Danuta Dobrzycka for comments on this paper. We also acknowledge the referee, K.R. Phillips, for suggestions that improved the clarity of this paper.


format: text

text:
# Untitled Document

Figure Captions

Figure 1

The figure represents the dimensionless expression of the total energy $`E^{}=RE/GN^2m^2`$ of isothermal spheres as a function of a dimensionless inverse temperature $`b=GNm^2/RkT`$, where $`k`$ is Boltzmann’s constant, and $`G`$ the gravitational constant. Numbers along the line represent density contrasts $`=\rho (0)/\rho (R)`$, $`\rho (0)`$ is the density at the centre, $`\rho (R)`$ is that at the boundary. $``$ grows continuously in the direction of the arrow.

Figure 2

The full line represents the schematic behaviour of the linear series $`E(\beta )`$ near point D with corresponding points C and E of equal energy $`E_C`$. The dashed line is a schematic curve representing $`w(E_\mathrm{C};x^1)`$ as a function of $`x^1`$, $`x^i(i>1)=`$constant.


text:
# Numerical studies of the extended two-chain model of friction

## I. INTRODUCTION

There are several ways for studying friction phenomena. Experimental observations are very powerful to characterize the frictional behaviour and to formulate empirical friction laws, e.g., Coulomb-Amonton’s laws , from a phenomenological point of view. Macroscopically, everybody can usually observe that during two solids are slided against each other, wear take place. But wear is not necessary for the occurrence of friction. In the last few years the technological development makes it possible to study wearless dry friction between atomically flat contact interfaces experimentally (nanotribology) . These nano-scale experiments on systems with less complicity stimulated again the theoretical efforts to investigate friction phenomena on the base of microscopic models . The dissipative nature of friction is a typical non-equilibrium problem. From the theoretical point of view, wearless dry friction between two atomic flat solids could be explained by using a simple driven mechanical many-body system with a lattice structure and take into account interaction forces. Thus the question arises, how the structure of such a microscopic model of friction has to look.

One of the simplest structure on the microscopic level is given by the discrete version of the well-known Frenkel-Kontorova (FK) model , which is in general a model of an one-dimensional harmonic chain of particles in a periodic potential. The FK model is mainly used for describing an adsorbate atom monolayer on an atomically flat clean substrate surface. Using this picture, the driven FK model, where a constant external force is applied to each particle of the chain, can also be used as a simple model structure for a wearless friction situation, which has been studied by several works . Note, several physical phenomena such as the dislocations motion in crystals , commensurate-incommensurate transitions  and charge-density waves (CDW‘s)  could also be investigated by reducing on the FK model structure. Definitely, the simple FK model is a good starting point in order to construct more complex and realistic friction models on the atomic level. Thus, a natural extension of the FK model is the replacement of the rigid substrat (hard body) by a deformable substrat monolayer pinned on a bulk (soft body). Such a more complex model of friction was firstly introduced by Matsukawa and Fukuyama . They considered an one-dimensional model consisting of two deformable interacting harmonic chains, where each particle of the lower chain is harmonically pinned on a rigid bulk. The static and kinetic properties of this two-chain model of friction was analysed in great details by Kawaguchi and Matsukawa . So far, the FK model and the two-chain model are simple models of wearless friction between an atomic monolayer and a hard or soft body, respectively. For describing a more realistic friction situation, which should represent wearless friction between atomically flat bodies, it is necessary to replace the adsorbate monolayer by a substrat. For this reason, Weiss and Elmer  proposed an interesting approach, the so-called Frenkel-Kontorova-Tomlinson (FKT) model, which is an one-dimensional lattice model for a soft upper body sliding on a hard lower body. In contrast to the two-chain model, with the FKT model they extended the FK model structure in another direction by harmonically coupling of each particle of the adsorbat monolayer on an upper macroscopic sliding mass, whose position relative to the lower surface is characterized by the coordinate $`x_B`$. Note, all of these models are investigated as completely one-dimensional microscopic friction models, where the particles of the one-dimensional chain can only move parallel to the direction of the external force. Now, it is obviously clear, how you have to extend the previous model structures in order to write down a more complex microscopic model system, which is necessary to describe friction phenomena as simple but also realistic as possible. For that purpose we introduce an extended two-chain model consisting of two interacting harmonic atomic chains embedded in a two-dimensional space, where each particle of the upper or lower chain is harmonically pinned on an upper or a lower rigid bulk, respectively. This model describes two soft bodies with an one-dimensional contact interface, which can carry out a genuine two-dimensional sliding motion against each other. In our studies we consider the extended two chain model as a simple relaxation model by taking energy dissipation via a phenomenological damping term in the mechanical equation of motion into account. Hence, all model, where the energy dissipation is modeled by hand through an explicit damping term, cannot explain the origin of friction in complex physical systems refering to the sliding state. Well, that is also not the aim of our investigations, because with our approach we want to establish the typical macroscopic friction situation (tension experiment) on a microscopic level. Nevertheless, the introduced model structure can be used to discuss the mechanism of the occurrence of static friction, which is one characteristic feature of dry friction, based on an energy concept. The hope is that the properties of the ground state - defined as the minimal potential energy state of the mechanical model system in absence of an external force - describe completely the static friction behaviour as function of the interaction between the two deformable chains. In general the ground state structure can exhibit two situations as a driven force increased from zero. On the one hand, a zero static friction arises exactly then if the maximum static friction force vanished, i.e., the model system is frictionless. On the other hand, a finite static frition corresponds always with the existence of a finite maximum static friction force, which is the maximum potential energy state of the system above which a sliding state appears. In other words, if we evolve the particle configuration from their ground state configuration, then the maximum static friction force is defined by the largest depinning force, which is equivalent to the smallest driven force above which no stable stationary state exists.

For the analysis of static properties of the extended two chain model we follow the philosophy of Aubry by using strangely enough an incommensurate ground state structure, which has the only purpose to introduce a special kind of mathematical functions called as hull-functions. On the basis of the specific hull-function concept for the visualisation of the ground state structure of the discrete FK-model in the incommensurate case Aubry  has been shown that a new type of phase transition (called as breaking of analyticity) occurs when the amplitude of the periodic potential increased beyond a critical value. That means, below a critical amplitude the hull-functions have a continuous (analytic) structure whereas above a critical amplitude the hull-functions are discrete (non analytic). In our numerical investigations of the hull-functions of the extended two chain model we have also found the transition by breaking of analyticity (Aubry transition). In order to distinguish, whether the Aubry transition is really about a second order phase transition, we have specifically studied this type of phase transition like Coppersmith and Fisher (but which have been investigate the behaviour near the maximum static friction force ) from the viewpoint of critical phenomena. For that purpose we have introduced a suitable order parameter, which shows the typical critical behaviour for second order phase transitions and critical exponents can be calculated via scaling arguments.

Here, it should mentioned, that from the physical intuition a finite maximum static friction as function of the interaction strength $`ϵ`$ has to exist in all previous introduced models equal whether one choose a commensurate or an incommensurate ground state structure. The main differences are the following. In the commensurate case the static friction is always larger than zero for any $`ϵ>0`$. By way of contrast, in the incommensurate case whether the static friction is equal or unequal zero depends strongly on $`ϵ`$ below or above a critical value $`ϵ_c`$ due to the preculiar breaking of analyticity. Hence the commensurate and the incommensurate ground state structure reveals two types of physical system classes characterizing by different critical behaviour in the sliding state near the threshold force $`F_c`$. The first class is given by model systems with a commensurate ground state structurce. Effectively, commensurate models are characterized by finitely many internal degrees of freedom, which lead to trivial critical exponents in the sliding state for the velocity-force characteristic near the threshold $`F_c`$. An incommensurate ground state structure is equivalent to a second class of models with infintely many internal degrees of freedom. Well, both model classes emphasize the fact that the nature of static and kinetic friction is deeply rooted with internal degrees of freedom generated by a complex atomic interaction behaviour, which is the key for unterstanding dissipation in general. Commensurate model have one main disadvantage, because an uniquely visualisation concept for the ground state structure like the hull-function for the incommensurate case is missing. Moreover, commensurate models exist always in a pinned ground state structure, that means, models imply static friction for any value of the amplitude of the interaction strength larger then zero intrinsicly. Therefore, for understanding the physical fundamental mechanism of the arising of static friction it seems to be most interesting to study incommensurate models, because the ground state shows two different incommensurate phases related to amazing properties. That means, an unpinned incommensurate phase corresponds to zero static friction whereas the pinned incommensurate phase is equivalent to a finite static friction.

In the present paper we introduce the extended two chain model and the numerical calculation methods. The detailed investigations of the static properties, especially the ground state structure of the extended two chain model are discussed by defining an order parameter and proposing scaling relations.

## II. THE MODEL SYSTEM

Following Matsukawa and Fukuyama , we introduce the extended two-chain model, which is a two-dimensional microscopic lattice model system describing wearless dry friction between two atomically flat solids. Therefore, we consider two deformable atomic chains, so called an upper and a lower chain. The atoms in each chain are classical point particles and are able to realize a two-dimensional movement. We assumed a harmonic particle interaction in each chain. The interaction between the two chains is of the Lennard-Jones type. Moreover, both chains are coupled harmonically on an upper or a lower substrat respectively. Like Ref.  the effects of energy dissipation in the upper and lower chain are phenomenologically approached to be proportional to the difference between the velocity of the $`i`$-th particle and that of the center of gravity. The external force $`\underset{¯}{F}^{ex}=(F^{ex},0)`$ acts on the upper chain in $`x`$ direction without $`y`$ component. In order to simplify the situation, we consider overdamped motion. The equations of motion for the particles in the upper and the lower chain are given for each direction as follows (see also FIG. 1)

$`m_a\gamma _a\left(\dot{x}_i^a\dot{x}_i^a\right)`$ $`=`$ $`k_a^x\left[x_{i1}^a+x_{i+1}^a2x_i^a\right]k_{sub}^{a_x}\left[x_i^aic_a\right]+{\displaystyle \underset{j=1}{\overset{N_b}{}}}F_{int}^x\left(x_i^ax_j^b\right)+F^{ex}`$ (1)
$`m_a\gamma _a\left(\dot{y}_i^a\dot{y}_i^a\right)`$ $`=`$ $`k_a^y\left[y_{i1}^a+y_{i+1}^a2y_i^a\right]k_{sub}^{a_y}\left[y_i^a(b_b+L_0)\right]+{\displaystyle \underset{j=1}{\overset{N_b}{}}}F_{int}^y\left(y_i^ay_j^b\right)`$ (2)
$`m_b\gamma _b\left(\dot{x}_i^b\dot{x}_i^b\right)`$ $`=`$ $`k_b^x\left[x_{i1}^b+x_{i+1}^b2x_i^b\right]k_{sub}^{b_x}\left[x_i^bic_b\right]+{\displaystyle \underset{j=1}{\overset{N_a}{}}}F_{int}^x\left(x_i^bx_j^a\right)`$ (4)
$`m_b\gamma _b\left(\dot{y}_i^b\dot{y}_i^b\right)`$ $`=`$ $`k_b^y\left[y_{i1}^b+y_{i+1}^b2y_i^b\right]k_{sub}^{b_y}\left[y_i^bb_b\right]+{\displaystyle \underset{j=1}{\overset{N_a}{}}}F_{int}^y\left(y_i^by_j^a\right),`$ (5)
where $`\underset{¯}{r}_i^a=(x_i^a,y_i^a)`$ and $`\underset{¯}{r}_i^b=(x_i^b,y_i^b)`$ are the immediately position vectors of the $`i`$-th particle in the upper and lower chain. Furthermore, $`N_a`$ ($`N_b`$), $`m_a`$ ($`m_b`$) is the particle number and the particle mass, $`\gamma _a`$ ($`\gamma _b`$) denotes the friction coefficient and $`c_a`$ ($`c_b`$) is the mean atomic spacing in $`x`$-direction of the upper (lower) chain. It should be mentioned, that the particle mass and the friction coefficient are for all particles in each chain the same. For a compact representation one can introduce the chain index $`I=(a,b)`$ relating to the upper and lower chain. The velocity of the center of gravity is given by the average $`\dot{X}_i^I=\frac{1}{N_I}_{i=1}^{N_I}\dot{X}_i^I`$, where $`\dot{X}_i^I`$ stands for $`\dot{x}_i^I`$ and $`\dot{y}_i^I`$ respectively. $`k_I^x`$ and $`k_I^y`$ are the strength of the harmonic interaction force in the chains for the $`x`$ and $`y`$ direction. The strenght of the coupling of each chain on an upper and a lower substrat is represented by $`k_{sub}^{I_x}`$ and $`k_{sub}^{I_y}`$, which determine the rigiditiy of the chains. A further system parameter is the substrat-substrat distance $`L=b_a+b_b+L_0`$, where $`b_a`$ ($`b_b`$) is the mean atomic spacings in $`y`$-direction between the upper (lower) chain and upper (lower) substrat. $`L_0`$ denotes the distance between the two harmonical chains and is a function of their interaction strength. The interaction force between the upper and lower chain, projected on the $`x`$ and $`y`$ axes, can be written as

$`{\displaystyle \underset{j=1}{\overset{N_J}{}}}F_{int}^x\left(x_i^Ix_j^J\right)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N_J}{}}}F(r_{ij}){\displaystyle \frac{\left(x_i^Ix_j^J\right)}{r_{ij}}}`$ (6)
$`{\displaystyle \underset{j=1}{\overset{N_J}{}}}F_{int}^y\left(y_i^Iy_j^J\right)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N_J}{}}}F(r_{ij}){\displaystyle \frac{\left(y_i^Iy_j^J\right)}{r_{ij}}}`$ (8)
with the chain indices $`I=(a,b)`$, $`J=(b,a)`$ and the distance $`r_{ij}=\left|\underset{¯}{r}_i^a\underset{¯}{r}_j^b\right|\sqrt{(x_i^bx_j^a)^2+(y_i^by_j^a)^2}`$ relating to the $`i`$-th atom of the upper (a) and the $`j`$-th atom of the lower (b) chain. $`F(r_{ij})`$ can be expressed by the potential $`V(r_{ij})`$ through $`F(r_{ij})=\frac{V(r_{ij})}{r_{ij}}`$. As interaction potential we choose specifically the Lennard-Jones potential

$$V(r_{ij})=4ϵ\left[\left(\frac{\sigma }{r_{ij}}\right)^{12}\left(\frac{\sigma }{r_{ij}}\right)^6\right]$$
(9)
representing the character of the friction due to the intermolecular forces in the best way. Here $`ϵ`$ is the strength of the interaction and $`\sigma `$ is a characteristic length parameter, which parameterized the minimum of the potential $`V^{min}=V(r_{ij})|_{r_{ij}=r_0}`$ via the distance from equilibrium $`r_0=\sigma \mathrm{\hspace{0.17em}2}^{1/6}`$.

### III. NUMERICAL METHOD

For solving the two-dimensional extended two-chain lattice model we choose periodic boundary conditions in both chains

$`x_i^I`$ $`=`$ $`x_{i+N_I}^IL_I`$ (10)
$`y_i^I`$ $`=`$ $`y_{i+N_I}^I`$ (11)
with the chain index $`I=a,b`$ and the chain lengths $`L_I=N_Ic_I`$. This implies that the system size of the upper and lower chain are the same $`LN_ac_a=N_bc_b`$. Hence, the ratio or the misfit of the two lattice constants could be expressed by the particle numbers

$$\alpha =\frac{c_a}{c_b}=\frac{N_b}{N_a}.$$
(12)
This means, that the ratio $`\alpha `$ of the two lattice constants is always a rational number in the numerical realisation. Using the concept of the hull-function, which was introduced by Aubry and coworkers  for describing the ground state of the underlying model system, we need an irrational ratio $`\alpha `$ between $`c_a`$ and $`c_b`$. If we have a rational $`\alpha `$, i.e., in the most commensurable case $`\alpha =1`$, we cannot apply the concept of the hull-function. Because the whole ground state configuration of each chain, consisted of $`N_I`$ particle positions, would be represented by one value of the argument of the hull function and hence the hull-function plot lose his mathematical power of proposition. For this reason, it is necessary to investigate an incommensurate ground state structure relating to the ratio $`\alpha `$ of the two lattice constants. Unfortunately, in the numerical calculation we could only approximate an irrational number through optimal rational ratios of two natural numbers using a finite continued-fraction expansion up to a certain order

$$\alpha =a_0+\frac{1}{a_1+{\displaystyle \frac{1}{a_2+{\displaystyle \frac{1}{\mathrm{}+{\displaystyle \frac{1}{a_{n1}+{\displaystyle \frac{1}{a_n}}}}}}}}}$$
(13)
Following the literature, we choose the simplest continued-fraction expansion with $`a_0=1`$ and $`a_1=a_2=\mathrm{}=a_n=1`$, so that $`\alpha `$ could be expressed by the ratios of the Fibonacci numbers $`F_{n+1}=F_n+F_{n1}`$, where $`n=1,2,3,\mathrm{}`$ as well as $`F_0=1`$ and $`F_1=1`$. Hence,

$$\alpha =\frac{F_{n+2}}{F_{n+1}}=\frac{3}{2},\frac{5}{3},\frac{8}{5},\frac{13}{8},\frac{21}{13},\frac{34}{21},\frac{55}{34},\frac{89}{55},\frac{144}{89},\frac{233}{144},\frac{377}{233},\frac{610}{377},\frac{987}{610},\frac{1597}{987},\mathrm{}$$
(14)
and the irrational limit $`\alpha lim_n\mathrm{}\frac{F_{n+2}}{F_{n+1}}=\frac{1}{2}\left(\sqrt{5}+1\right)`$ of the approximation series (14) is the inverse of the so-called golden mean. From the point of number theory  , this is the most incommensurable case of the underlying model system, which could be emulated in a systematical way. In our studies we set the lattice constant of the lower chain equal one ($`c_b=1`$). Accordingly, the upper mean lattice spacing $`c_a`$ is equal $`\alpha `$ and we calculated with the system sizes

$$c_a=\frac{N_b}{N_a}=\frac{89}{55},\frac{233}{144},\frac{377}{233},\frac{610}{377},\frac{1597}{987}.$$
(15)
It should be mentioned, that not all ratios oft the approximants are ”good” irrational number approximations in the same manner because for certain particle number ratios there exist some artifacts in the hull function.
For the numerical calculation of the equations of motion we used a fourth-order Runge-Kutta algorithm with an adaptive stepsize control. In order to discuss the essential features of the present model we restrict the wide range of model parameters by setting $`m_I=1`$, $`\gamma _I=1`$, $`k_I^x=k_I^y=1`$, $`k_{sub}^{I_x}=k_{sub}^{I_y}=1`$ and $`\sigma =1`$. As initial conditions we assumed at $`t=0`$ a particle configuration in the x-y-plane, where the $`N_a`$ and $`N_b`$ chain particles are arranged at regular lattice sites periodically. Besides, the distance $`L_0`$ is set to an arbitrary initial value, i.e. $`L_0=1`$. For calculating the stable stationary states we used the same criterion like Kawaguchi and Matsukawa . That means, we used a velocitiy condition

$$\sqrt{\frac{\left(_{i=1}^{N_a}(\dot{r}_i^a)^2+_{j=1}^{N_b}(\dot{r}_j^b)^2\right)}{(N_a+N_b)}}<\mathrm{\hspace{0.17em}10}^5$$
(16)
for stopping the Runge-Kutta calculation. Hence, the system could be considered to be in a static state. The next step consists in the check up of the total sum of the coupling forces in the (a)-chain and (b)-chain in $`x`$-direction

$$F_{sub}^{a_x}=\underset{i=1}{\overset{N_a}{}}k_{sub}^{a_x}\left[x_i^aic_a\right],F_{sub}^{b_x}=\underset{j=1}{\overset{N_b}{}}k_{sub}^{b_x}\left[x_j^bjc_b\right]$$
(17)
and $`y`$-direction

$$F_{sub}^{a_y}=\underset{i=1}{\overset{N_a}{}}k_{sub}^{a_y}\left[y_i^a(b_b+L_0)\right],F_{sub}^{b_y}=\underset{j=1}{\overset{N_b}{}}k_{sub}^{b_y}\left[y_j^bb_b\right]$$
(18)
respectively. In absence of a driven force ($`F^{ex}=0`$) it is sufficient to consider only the value of the total sum of the coupling forces (18) in $`y`$-direction, because for $`t=0`$, after the Runge-Kutta steps and in the final state the $`x`$-direction forces (17) are equal to zero. The regular periodical arrangement of the particles for the initial time $`t=0`$ leads to $`F_{sub}^{a_y}=F_{sub}^{b_y}=0`$. After the relaxation in the final stationary state relating to the velocity condition (16) the total sum of coupling forces (18) must always be zero, which is in general not true for an arbitrary chosen initial value $`L_0`$. Accordingly, during the numerical calculation every final stable state in respect to $`L_0`$ must be proofed and if necessary $`L_0`$ must be changed like $`Ł_0=L_0\pm \delta L_0`$. Then the chains relax again and the procedure would be repeated until the system reached a final stationary state, which fulfil the condititon $`F_{sub}^{a_y}=F_{sub}^{b_y}<510^5(0)`$. In this way we could find numerically the stationary configuration of the extended two-chain model and for $`F^{ex}=0`$ one can assumed that this is approximately the ground state, defined as having an energy which cannot be lowered by moving of some system particles. The properties of the ground state could be analysed by changing the interaction strength $`ϵ`$, which corresponds to set the correct $`L_0=L_0(ϵ)`$ numerically.
In the presence of an external force $`F^{ex}>0`$ acting in $`x`$-direction one must slightly change the numerical calculation procedure. Instead of $`F_{sub}^{a_x}=F_{sub}^{b_x}=0`$ the $`x`$-direction forces (17) have for $`t=0`$ the values $`F_{sub}^{a_x}=F^{ex}`$ and $`F_{sub}^{b_x}=0`$ respectively. After each Runge-Kutta step we have to guarantee that the control parameters (17) have again these values. Hence, we are able to evolve the model system from the ground state into the last possible static solution by continuously increasing of the applied external force. In this sense we could numerically define the threshold, where no longer a static solution of the underlying system exist. This threshold $`F^{ex}=F_c`$ is the critical force of static friction, which is also called maximum static frictional force. Below $`F^{ex}<F_c`$ the chains are in a steady state, such that the afterwards space-time averaged velocity $`v\dot{x}_i^a_{i,t}=0`$. Above $`F^{ex}>F_c`$ the upper chain, which is after an initial transient decay in a steady sliding state, slides with a constant velocity $`v`$. Note, in an easy way one can investigate the present model in a gravitational field via the constraints (18).

## IV. NUMERICAL RESULTS

### A. Ground state

First, we investigate the stationary states of the extended two-chain model in the absence of an external force, i.e., $`F^{ex}=0`$. The stationary states are given by the solution of equation system (1-5) with vanishing right-hand side

$$0=m_a\gamma _a\left(\underset{¯}{\overset{\dot{}}{r}}_i^a\underset{¯}{\overset{\dot{}}{r}}_i^a\right)=m_b\gamma _b\left(\underset{¯}{\overset{\dot{}}{r}}_i^b\underset{¯}{\overset{\dot{}}{r}}_i^b\right)$$
(19)
with consideration of the conditions (17) and (18), where $`F_{sub}^{a_x}=F_{sub}^{a_y}=F_{sub}^{b_x}=F_{sub}^{b_x}=0`$. Hence, the ground state of the extended two chain model is the stationary particle configuration having the minimal potential energy for $`F^{ex}=0`$. The knowledge of the ground state determine a physical system, such as the static friction, completely.

#### 1. The hull-function

Similarly to the FK model  or FKT model , the lattice structure of the particle positions in the incommensurate ground state can be uniformly described by a hull-function. For examining the lattice structure of the present two-dimensional extended two-chain model, it is necessary to introduce four hull-function due to two degrees of freedom of the particles in each chain. The two hull-functions $`h_a`$ and $`h_c`$ of the upper chain are defined as

$`x_i^a`$ $`=`$ $`ic_a+\varphi +h_a(ic_a+\varphi )`$ (20)
$`y_i^a`$ $`=`$ $`b_b+L_0+\varphi ^{}+h_c(ic_a+\varphi ^{})`$ (21)
and the two hull-function $`h_b`$ and $`h_c`$ of the lower chain are given by

$`x_i^b`$ $`=`$ $`ic_b+\psi +h_b(ic_b+\psi )`$ (22)
$`y_i^b`$ $`=`$ $`b_b+\psi ^{}+h_d(ic_b+\psi ^{}),`$ (23)
where $`\varphi `$, $`\varphi ^{}`$, $`\psi `$ and $`\psi ^{}`$ are constant phases. For our investigations, it is conviencent to set $`\varphi =\varphi ^{}=\psi =\psi ^{}=0`$. The hull-functions are periodic and even, which can usually be expressed by

$`h_a(z+c_b)`$ $`=`$ $`h_a(z)=h_a(z),h_c(z+c_b)=h_c(z)=h_c(z)`$ (24)
$`h_b(z+c_a)`$ $`=`$ $`h_b(z)=h_b(z),h_d(z+c_a)=h_d(z)=h_d(z).`$ (25)
The argument $`z`$ of $`h_a`$ and $`h_c`$ is defined in the range from $`0`$ to $`c_b`$, whereas of $`h_c`$ and $`h_d`$ it is given from $`0`$ to $`c_a`$. Note, the visualisation concept of the hull-function makes it necessary that the argument of the just defined hull-functions $`h_a`$, $`h_c`$ or $`h_b`$, $`h_d`$ must be the same, respectively. Because, the displacements in $`x`$\- resp. $`y`$-direction of all particle in each chain can be assigned to different values of the argument $`z`$ due to the irrational ratio $`\alpha `$ of the mean lattice spacing $`c_a`$ and $`c_b`$.

First, we have numerically calculated the ground state of the extended two-chain model for a fixed ratio $`\alpha =N_b/N_a`$ of the particle numbers given by $`N_a=233`$ and $`N_b=377`$. Several investigations of the FK model , FKT model  such as the two-chain model  showed, that the ground state strongly depends on the interaction strength $`ϵ`$. Therefore, we started the computation for different values of $`ϵ`$ in the wide range from $`0.05`$ to $`0.4`$. For small values of $`ϵ`$ up to a certain value $`ϵ_c`$ we get the result that the periodic hull-functions $`h_a`$, $`h_b`$, $`h_c`$ and $`h_d`$ are analytic, i.e, smooth and continuous. This behaviour is illustrated in Fig. 2, where the interaction strength has the value $`ϵ=0.22`$. Furthermore, the stationary state for $`ϵ=0.22`$ is characterized by the distance $`L_0=0.986273`$ between the upper and lower chain. Figure 3 shows another numerical calculation of the hull-functions for the interaction strength $`ϵ=0.33`$ with the corresponding distance $`L_0=0.959584`$. We can see that the hull-functions in Fig. 3 are discrete, which means that they are no longer analytic. This discrete structure is related to jumps, which occur at certain points of the argument $`z`$ of the hull-functions. Now, the comparison of Fig. 2 with Fig. 3 shows the phenomenon of breaking of analyticity, which is due to two significant aspects of the hull-functions $`h_a`$, $`h_b`$, $`h_c`$ and $`h_d`$ for the interaction strength $`ϵ<ϵ_c`$ and $`ϵ>ϵ_c`$, which defined a certain transition value $`ϵ_c`$ above the hull-functions are no longer analytic. For a detailed analysis of the threshold $`ϵ_c`$ we introduced a suitable order parameter in the following part.

Well, Aubry has shown  that the breaking of analyticity of the incommensurate ground state of the FK model is characterized by the existence of a largest central gap at the half of the period of the hull-function. Whereas the situation in the two-chain model of Matsukawa and Fukuyama  is quite different for several elastic parameters and three simple cases are observed: both hull-functions $`h_a`$ and $`h_b`$ have a largest gap at the half of the period (FK limes), whether $`h_a`$ or $`h_b`$ has only a largest central gap and finally both hull-functions do not have a largest central gap but two symmetrical gaps with the same size refering to the half of the period. In Fig. 3 we can also observe that no largest central gap for all hull-functions of the extended two-chain model exists. Furthermore, the hull-functions show also two symmetrical gaps of the same size refering to the half of the period of the hull-function in the whole analysed parameter range. Up to now we have only investigated the ground state structure for a fixed system size ($`N_a=233`$ and $`N_b=377`$) and a certain range of interaction strength $`ϵ`$. Now, we vary the system size of the underlying model which corresponds the different approximation orders of the irrational number through a rational ratio of the particle numbers. Here, we calculated especially the ground state structure for the ratios $`\alpha =89/55`$, $`233/144`$, $`610/377`$, and $`1597/987`$ in the same range of the interaction strength $`ϵ`$ like above. In particular, we proofed whether there occurs an anormalous behaviour in the structure of the hull-functions for certain values of $`ϵ`$ due to the using of a rational ratio instead of an irrational number in the numerical computation. Indeed, Fig. 4 shows such a anormalous behaviour of the hull-functions for $`ϵ=0.33`$ and for the particle numbers $`N_a=377`$ and $`610`$. In the comparison with Fig. 3 we can recognize that the even particle number $`N_b=610`$ generates a wrong particle position at the half of the period $`x=0.5`$ of $`h_b`$ in Fig. 4 due to the periodical boundary condition. Consequently, in term of a suitable rational particle ratio Fig. 3 shows the right picture of the hull-functions for $`ϵ=0.33`$, whereas Fig. 4 shows some artifacts in the hull-function structure for the same interaction strength $`ϵ=0.33`$ as a result of bad choice $`N_a=377`$ and $`N_b=610`$ of the system size for the interaction strength $`ϵ=0.33`$. Beyond it, there appear much more artifacts in the hull-function structure different ranges of $`ϵ`$ als consequence of a bad choice for the particle number ratios. For $`N_a=144`$ and $`N_b=233`$ the occurence of such artifacts destroyed the order parameter analysis of the transition area. Kawaguchi and Matsukawa interpreted in  these artifacts, which also occur in the two-chain model for certain values of the interchain strength and $`N_a=144`$, $`N_b=233`$, as information of the appearance of a new gap structure at the elastic parametes $`K_b`$. This interpretation is wrong, because, if they calculated the model system with the same interaction strength for different particle numbers these artifacts will be vanished. The appearence of artifacts in the hull-function for the ground state defines the range of validity of the hull-function concept for several $`ϵ`$ in the numerical calculation, which makes it necessary to proof the hull-function structure for the same $`ϵ`$ for several system size. Relating to stationary state which is given by Eq.(19), Fig. 5 shows the distance $`L_0`$ between the upper and the lower chain as a function of the strength $`ϵ`$ of the interaction force, which was calculated for the particle numbers $`N_a=233`$ and $`N_b=377`$. The plot in Fig. 5(a) shows the computation of $`L_0`$ for $`ϵ`$ in range from $`0.05`$ to $`0.45`$ and Fig.5(b) shows $`L_0`$ in the transition area. Note, the distance $`L_0`$ is remarkably insensitiv to the system size. For the system sizes $`\alpha =377/233`$, $`610/377`$, and $`1597/987`$ we found a deviation for $`L_0=L_0(\alpha )`$ smaller than $`10^5`$. In the occurence of local artifacts in the hull-functions for certain particle numbers and certain interaction strength the deviation for $`L_0(\alpha )`$ lies in the order $`10^4`$. Note, that in the ground state above the transition point $`ϵ_c`$ also many local rearrangements take place which can be identified as first order phase transitions.

#### 2. The order parameter

Now, we analyse the transition of breaking the analyticity near the threshold $`ϵ_c`$ in terms of second order phase transition. For this it is convenient to introduce an order parameter, which is an indicator for the occurence of a phase transition at a critical point. The order parameter distinguishs both phases which are involved in the phase transition, because it is defined to be zero in one phase. Refering to the transition of breaking of analyticity the order parameter is defined in the discrete phase, which measures the discrete jumps in the hull-functions. Therefore the order parameter is in the continuous phase for $`ϵ<ϵ_c`$ equal to zero and in the discrete phase above $`ϵ_c`$ unlike zero. Here, we use the definition

$$OP_m=\frac{{\displaystyle _0^1}\left|{\displaystyle \frac{^mf(z)}{z^m}}\right|^2𝑑z}{\delta (0)}.$$
(26)
Here, $`\delta (0)`$ behaves like $`\delta (0)\mathrm{\Delta }z^1|_{\mathrm{\Delta }z0}`$ . If the hull-function is smooth in the interval $`[0,1]`$, the integral becomes finite and the order parameter vanished due to the prefactor $`\delta (0)^1`$. On the other hand, if $`f(z)`$ has a discrete step at $`z_0`$, we obtain $`_zf(z)\delta (zz_0)`$ and therefore $`_0^1[_zf(z)]^2𝑑z_0^1\delta (zz_0)^2𝑑z\delta (0)_0^1\delta (zz_0)\delta (0)`$, i.e. the order parameter becomes finite. Due to the finite system size, which devides the region of integration in equidistant subintervals $`\mathrm{\Delta }z=1/N_a`$ ($`\mathrm{\Delta }z=1/N_b`$) for the upper (lower) chain, the integral in Eq.(26) reduced to a sum over the particle numbers. Therefore, we have intrinsic discontinuities in the numerical computation, because the hull-functions are given by a discrete set of points at $`\mathrm{\Delta }z`$ at all, which define the systematical error of the calculation. In order to smeared out the systematical intrinsic discontinuities, we choose for our numerical investigation $`OP_4`$ as order parameter. Hence, we get the following expression

$$OPOP_4=\frac{{\displaystyle _0^1}\left|{\displaystyle \frac{^4f(z)}{z^4}}\right|^2𝑑z}{\delta (0)}\frac{1}{(\mathrm{\Delta }z)^6}\underset{k=2}{\overset{\genfrac{}{}{0pt}{}{(N_b2)}{N_a2}}{}}\left(\left|f_{k+2}4f_{k+1}+6f_k4f_{k1}+f_{k2}\right|\right)^2$$
(27)
It should mentioned, that $`OP`$ shows necessary for all hull-functions $`h_a`$, $`h_c`$, $`h_c`$ and $`h_d`$ the same qualitativ behaviour. Therefore, in order to show all basic features of the phase transition (transition of breaking of analyticity) it is enough to calculate $`OP`$ (via Eq.(27)) only for the hull-function $`h_a=h_a(z)`$.

#### 3. The scaling behaviour

Next we have calculated the order parameter given by Eq.(27) for different system sizes $`\alpha =89/55`$, $`233/377`$, $`610/377`$ and $`1597/987`$. Figure 6(a) shows the order parameter as function of the interaction strength $`ϵ`$ in the transition region. We clearly see that the order parameter $`OP`$ is changed from zero below a certain critical threshold $`ϵ_c`$ to a value unlike zero above $`ϵ_c`$ for all particle number ratios. This is a certain hint that the transition of breaking of analyticity is a second order phase transition. Using a finite-size scaling  for the data in Fig. 6(a) given by the scaling function

$$OP(ϵϵ_c,L)=L^{\frac{\beta }{\delta }}OP\left(L^{\frac{1}{\delta }}(ϵϵ_c)\right)$$
(28)
we find that the order parameter $`OP`$ exhibits a scaling behaviour in the transition region shown in Fig. 6(b). That means, all curves of Fig. 6(a) could mastered on one resulting curve (Fig. 6(b)). The system size $`L`$ is determined by the particle numbers $`55`$, $`233`$, $`377`$ and $`987`$. In the transition region, which can be characterized by the reduced interaction strength $`t=(ϵϵ_c)/ϵ_c`$, we find that the data of Fig. 6(b) are consistent with the following power law for $`ϵ>ϵ_c`$

$$OP(ϵϵ_c)L^{\frac{\beta }{\delta }}\left(ϵϵ_c\right)^\beta $$
(29)
with the critical exponents $`\beta =0.239\pm 0.05`$, $`\delta =0.687\pm 0.05`$ defining the ratio $`\mu =\beta /\delta 0.347\pm 0.05`$. Due to the finite-size scaling we obtain for the threshold interaction strength of the transition of breaking of analyticity the critical value $`ϵ_c=0.25575\pm 0.001`$.

### B. The static friction

In the following we want to calculate the maximum static friction force $`F_c`$ of the extended two-chain model refering to a pinning structure which is given in the incommensurate case for $`ϵ>ϵ_c`$. It is clear that in the presence of an applied forces $`F^{ex}>0`$ ground states are no longer defined. From the physical intuition, we hope that one can adiabatically evolve the underlying model system from the ground state by increasing $`F^{ex}`$ from zero, because the threshold force $`F_c`$ is the largest depinning force $`F^{ex}`$ above no a longer stationary state exist, i.e., Eq.(19) is no longer valid. Like Ref.  the static friction is given by

$$F_R=\underset{i=1}{\overset{N_a}{}}\underset{j=1}{\overset{N_b}{}}F_{int}\left(x_i^ax_j^b\right)_tN_aF^{ex}_t$$
(30)
which can be obtained from the equation of motion (1)-(5) by summing over $`i`$ and averaging with respect to $`t`$. It should remarked, that Eq.(30) is also valid for the calculation of the kinetic friction in a sliding state, i.e., for $`F^{ex}>F_c`$. For the determination of the threshold force $`F_c`$, we choose especially the interaction force $`ϵ=0.33`$ and the system size $`N_a=233`$ and $`N_b=377`$. Note, the numerical calculations of $`F_c`$ for the extended two-chain model are very hard and the computation time is very large. In a first calculation, we find that the maximum static friction force has approximately the value $`F_c0.4`$ but the numerical deviation is very large. Therefore, it takes much more efforts in order to calculate the accurate value of the threshold force.

## V. CONCLUSIONS

We have investigated the ground state for a two-dimensional extended two-chain model with an incommensurate lattice structure. We have found that the structure of the hull-functions depends strongly on the strength $`ϵ`$ of the interaction between the upper and the lower body. It was numerically observed the breaking of analyticity in the ground state for the critical value $`ϵ_c=0.25575`$. Furthermore, it was shown that this type of phase transition could be uniform described by an order parameter. This order parameter representation allows a systematic analysis of the critical behaviour for $`ϵ>ϵ_c`$ due to the calculation of the critical exponents. Consequently, the peculiar breaking of analyticity shows all features of a second order phase transition.

A second interesting result of the ground state structure occurs in associtation with the numerical computation for different interaction strengths as function of the rational system sizes. It has been shown that artifacts in the hull-functions structure appear for certain rational particle numbers of a fixed interaction strength $`ϵ`$. This fact pointed out that the correct ground state structure depends strongly on the choose of sufficient rational ratio of the particle numbers as a irrational number approximation for a certain range of the interaction force.

Concerning the static properties of the extended two-chain model in the presence of an applied forces $`F^{ex}>0`$ we have calculated the maximum static friction force $`F_c`$.

Finally, it should remarked, that the extended two-chain model is a simple model structure , which established the macroscopic dry friction situation on the microscopic level between atomically flat surfaces by using a simple driven mechanical many-body system. It is convenient to define the model as simple as possible in order to dicuss some basic features like the ground state structure or the maximum static friction behaviour. Therefore we took in our investigation a simple relaxation dynamic into account and neglected thermal fluctuation ($`T=0K`$). Surely, such kind of microscopic models could be always extended by additional thermal fluctuations representing by stochastical forces which leads to creep in the presence of applied forces. Note, the vality of the underlying model system is given by the assumption of simple harmonic particle interaction. It is clear that this assumption in many real systems where large distortions occur.

Note, there are one possible and very interesting extention of the extended two-chain model by the introduction of quenched randomess which corresponds surface impurities. It is an open question, whether randomess destroy the frictionless state in the incommensurate case for $`ϵ<ϵ_c`$ in the extended two-chain model like in the FK-model . Well, it is an interesting question, whether randomess, which exists in any real system, leads in general to a finite maximum static friction force for any interaction strength $`ϵ`$. Refering to the compare whith experiments, another natural extention of the present model is to proceed to a three-dimensional system, which is able to describe wearless dry friction between two atomically flat plans pinning in each case on a substrat and embeding in a three-dimensional space.

In later works we want to investigate the threshold behaviour near $`F_c`$ especially for $`F^{ex}>F_c`$. It should be proofed whether the velocity-force characteristics shows a critical behaviour or not. Well, for a complete understanding of the behaviour near the threshold $`F_c`$ a renormalization group calculation is needed, but surely very complicated for a deterministic mechanical system.

Note in the present paper, we have considered a simple overdamped motion. Hence, it is also interesting to discuss the influence of inertia terms (underdamped dynamic) on the behaviour near the threshold $`F_c`$ and on the sliding state.

## ACKNOWLEDGMENTS

This work has been supported by the Graduiertenkolleg ”Defektstrukturbestimmte physikalische Eigenschaften” and the Deutsche Forschungsgemeinschaft DFG (SFB 418).

FIGURES

FIG. 1. A schematic picture of the two-dimensional extended two-chain model. The external force $`\stackrel{~}{F}^{ex}`$ is given by the applied force $`F^{Fex}`$ which acts on each particle of the upper chain times the particle number $`N_a`$ of the upper chain where $`N_a`$, ($`N_b`$) is the particle number of upper (lower) chain.

FIG. 2. The hull-functions of the ground state plotted for the particle numbers $`N_a=233`$ and $`N_b=377`$ with the interaction strength $`ϵ=0.22`$. For an uniform representation all hull-functions are renormalized by the mean lattice spacing $`c_b`$ and $`c_a`$, respectively. Therefore, the substitution for the hull-function $`h_a`$ and $`h_c`$ of the (a)-chain is $`zz^{}=z/c_b`$, whereas for the hull-functions $`h_b`$ and $`h_d`$ of the (b)-chain it is $`zz^{}=z/c_a`$.

FIG. 3. Hull-functions of the ground state for $`ϵ=0.33`$ and the particle numbers $`N_a=233`$ and $`N_b=377`$.

FIG. 4. Artifacts in the hull-functions for a bad choice of the particle numbers, e.g., given by $`N_a=377`$ and $`N_b=610`$ and the interaction strength $`ϵ=0.33`$.

FIG. 5. The interaction distance $`L_0`$ in $`y`$-direction between the upper and lower chain as function of the interaction strength $`ϵ`$. The distance $`L_0`$ is plotted in different ranges (a) from $`ϵ=0.05,\mathrm{},0.4`$ and (b) from $`ϵ=0.24,\mathrm{},0.31`$ (transition range).

FIG. 6. In the computation it is convenient to rescale the order parameter by $`OP`$ $``$ $`OP(\mathrm{\Delta }z)^6`$. The order parameter $`OP=OP(ϵϵ_c)`$ describes the critical behaviour near the threshold $`ϵ_c`$. The graph (a) shows the behaviour of the order parameter for different particel numbers and (b) is the finite-size scaling plot of the mastered order parameter, where $`OP(ϵ)L^{\frac{\beta }{\delta }}`$ is plotted versus $`(ϵϵ_c)L^{\frac{1}{\delta }}`$.