Amara Touil
email: amara.touil@univ-brest.fr

Jean Vareille
email: jean.vareille@univ-brest.fr

Fred Lherminier
email: fredl@terranov.com

Philippe Le Parc
email: philippe.le-parc@univ-brest.fr

Modeling and Analysing Ubiquitous Systems Using MDE Approach

Keywords: Domain Specific Language, Model Driven Engineering, Analysis, Telecontrol Systems

The growth of industrial activities during the last decades and the diversity of industrial products require standards and common methodologies for building and integrating different parts. It is also required that working groups use the same terminologies and concepts needed for each domain. The Model Driven Engineering approach aims to give an answer while using a high level method based on models and transformations. In this paper, we use this approach to model ubiquitous systems. Those systems are composed of devices interconnected through various kinds of network and offer to get and send information. We present a model for this class of system and its use in the field of energy while studying real cases from our industrial partner Terra Nova Energy. This company aims to give solutions to monitor energy use and to reduce consumption. First results, where the use of a Model Driven Engineering approach makes it possible for our partner to improve and to get other points of view of his systems, are presented.

I. INTRODUCTION

Ubiquity is often define as being able to be in several places at the same time. In the field of Information Technologies, this definition can be refined in, at least, two ways that may be apparently opposites. The first approach is to considered as in [START_REF] Weiser | The computer for the 21st century[END_REF] that users are surrounded with "intelligent" systems, that may delivers them needed information: in this case, computing facilities are used to locate users, to understand their environment, to anticipate their needs and to make it possible that all information they require is available anywherere at anytime. The second approach, is to offer people to be able to get information from a system located somewhere and to be able to act safely on it: in this case computing facilities are used to make it possible to be "virtually" present in several places at the same time. We place our work in the second approach in order to model and analyze telecontrol application as a kind of ubiquitous systems.

Terra Nova Energy (TNE) is an innovative company that provides solutions for data mining in the fields of electricity, hydraulic and pulse-energy. It integrates sensing, acting and communicating devices in telecontrol systems, for collecting data from different sites using various technologies. TNE's systems allow real-time operating and provide processes for handling different data treatments.

In order to improve its development, this company is facing two problems: how to conceive remote monitoring systems as automatically as possible and how to fasten their on-site deployment ?

To answer these questions, we propose to use a model based approach, relying on specific components for telecontrol domain, and libraries to integrate industrial parts coming from various providers. Our idea is to build a framework [START_REF] Le Parc | A model-driven approach for building ubiquitous applications[END_REF] following the Model Driven Engineering (MDE) methodology [START_REF] Kent | Model driven engineering[END_REF] [START_REF] Fleurey | Validation in model-driven engineering: testing model transformations[END_REF] to create a telecontrol ontology and proper transformations for building, generating, analyzing and deploying such ubiquitous telecontrol systems.

Our aim is to define a generic meta-model and to use it for various systems, as case study examples, one of them being the TNE's system. This paper is organized as follows: first of all, the MDE approach, Domain Specific Languages (DSL) and transformations that may be conducted are briefly introduced. Second, we introduce the generic meta-model proposal. Third, this general meta-model is applied to our case study, while defining an instance for TNE. Then, we explain how to carry a static analyzis on a given instance and show first results. Finally we discuss our method and we finish by further work.

II. RELATED WORKS

In this section, we introduce some basic notions about the MDE approach and DSL, and how to conduct transformations on models.

A. The MDE approach

A model is an abstracted view of a system that expresses related knowledge and information. It is defined by a set of rules, used to interpret the meaning of its components [START_REF] Gerard | Using the model paradigm for real-time systems development: Accord/uml[END_REF] [START_REF] Moody | Graphical Entity Relationship Models: Towards a More User Understandable Representation of Data[END_REF]. A model is defined according to a modeling language that can give a formal or a semi formal meaning description of a system depending on modeler's intention. Modeling languages can be textual or graphical.

The model paradigm has gained in importance in the field of systems engineering since the nineties. Its breakthrough was favoured by working groups like the Object Management Group (OMG) [START_REF] Cranefield | UML as an ontology modelling language[END_REF] that has normalized modeling languages such as Unified Modeling Language (UML) and its profiles (such as SysML and MARTE for real-time systems). This group also provides the Model-driven architecture (MDA) software design standard. The MDA is the main initiative for Model Driven Engineering (MDE) approach.

According to OMG, four abstraction levels have to be considered: a meta-meta-model that represents the modeling language, which is able to describe itself; a meta-model level that represents an abstraction of a domain, which is defined through the meta-meta-model; a model level that gives an abstraction of a system as an instance of the meta-model; finally, the last abstraction level is the concrete system.

Tools have been defined to implement OMG standards. Some have general and wide purpose like those for UML language, others have been defined for specific and reduced class of applications like in [START_REF] Losilla | Wireless sensor network application development: An architecturecentric mde approach[END_REF] that aims to specify wireless sensor network systems in order to generate code. Another approach is to use Domain Specific Language (DSL) to specify a special semantic and/or syntactic rules for a class of application. One of DSLs definition is given by [START_REF] Van Deursen | Domain-specific languages: an annotated bibliography[END_REF]: "A domain-specific language (DSL) is a programming language or executable specification language that offers, through appropriate notations and abstractions, expressive power focused on, and usually restricted to, a particular problem domain". In our case we define our own DSL for telecontrol, according to Eclipse Core (Ecore) meta-model [START_REF] Gronback | Eclipse Modeling Project: A Domain-Specific Language (DSL) Toolkit[END_REF].

B. Model transformation

In order to breath life in models [START_REF] Czarnecki | Classification of model transformation approaches[END_REF], model transformations aims to exceed the contemplative model view to a more productive approach, towards code generation, analysis and test, simulation, etc. Models are transformed into other models that may be manipulated by specific tools or that may be transformed again into other models. Two kinds of transformation are used, exogenous if source and target transformation do not have the same meta-model and endogenous if source and target have the same meta-model. We use the latter here in order to perform static analysis of the built models.

Generally, static analysis deals with observing system and checking some properties before real system creation, production or code implementation. As example, formal verification and model checking techniques [START_REF] Alenljung | Formal Verification of PLC Controlled Systems Using Sensor Graphs[END_REF][13] are used for verifying models. Behavioural analysis are not beyond in the scope of this article. Static analysis are only conducted in simple cases, but relevant ones for TNE. In this paper we will focus on placement but others studies have been conducted (cost, bandwith,...). As soon as the model is created, it may give answers to the engineer: Are my sensors, repeaters, routers placed in a proper location? Are there better configuration? etc.. Concerning the first question, there are several works treating these points like [START_REF] Ghosh | Sensor network design for smart highways[END_REF], which study Wireless Network Sensors (WSN) placement for smart highways; it gives a solution based on geometric resolution algorithms in outdoor and open space. In our case, indoor deployment environments, including various fix and mobile obstacles, have to be studied. Component's placement depends on monitored zone, on the chosen topology, on the network capacity, etc. The second question can be adressed by optimyzing the placement configuration.

Next section describes the meta-model proposed for telecontrol and its derivation for TNE's remote monitoring needs.

III. PROPOSED META-MODEL

In this section, we, describe the meta-model we are proposing to specify systems based on communicating objects. Only a high level view and the main concepts are presented. Starting from this meta-model, we derive a specific instance for TNE.

A. Generic meta-model

At a first abstraction level, we built our domain around systems containing entities that communicate with each other as stated in [START_REF] Committee | A consensus of the incose fellows[END_REF]: "A system is a construct or collection of different elements that together produce results not obtainable by the elements alone". In the meta-model proposal (Figure 1), a system, modeled by System, can contain other systems according to the ownedSystem reference, in order to provide the "system of the system" notion. Regarding entities, they are modeled by Entity and its system containment is referenced by the itsEntity relationship. Entity paradigm in our meta-model aims to be a generic and general concept formed by extracting common features from our specific telecontrol domain. It can be a logical or physical and can be also composed of other entities (ownedEntity).

Moreover, entities are described using several related concepts trying to keep information about their physical properties, communication facilities or behaviour. These concepts are modeled as follows:

• Structure element defines the interrelation or arrangement of different physical parts or the organization of elements that provides coherence, shape and rigidity to an entity. It may describe mechanical, chemical or biological aspects. • Actions performed by entities are described by the Task concept. They can be internal, such as making computation, external such as sending information and also connected to the real world such as sensing or acting on a real device. • ContactInterface element allows connection between entities. It can be a physical contact, or a logical interface depending on the specification level, on the refinement degree or on its own structure. • An entity may contain Data related to itself or to its environment.

• Message element provides data exchange between entities. To send (or to receive a message) from entity A to entity B, A and B must be connected through ContactInterfaces, directly or, it may exist a path of entities that may forward messages. • To each entity, a Behavior concept is added, in order to describe the different tasks an entity may realize. From this concept, code generation or simulation may be refined.

As the Entity concept is very general, in order to improve our generic meta-model, different sub-entities have been specialized using the concept of heritage:

1) Device represents a physical component such as a sensor, a router, a computing unit, etc. 2) Medium element represents an entity deployed and used between other entities to specify communication's medium. In classical approaches, medium is often ommitted and considers as perfect links. In our work, medium has its own structure, behavior, contact interfaces, datas, messages and tasks. Different mediums may be described and classified following their own specifications, wired or wireless for example. 3) User element is used to model a person interacting with other entities, it is a kind of human avatar. Like other entities User can have different presentation depending on model view. 4) The different entities composing the system are under some physical conditions and constraints that affect and influence the global behaviour. In our meta-model we introduce the Environment concept to describe such conditions. Because we are dealing with a complex system, the Environment is an entity of this system and not just an external element. In our modeling approach, Environment acts on internal elements in order to build a global circumstances of evolving system and gives context for ubiquitous entities and systems that are context-aware. 5) The concept of time is crucial for telecontrol systems that is why a TimeClock entity is added to our generic meta-model for measuring, recording or indicating time.

The presented generic meta-model intends to specify any kind of systems containing communicating objects. At this abstraction level, we present just a first facet (also named Basic package), without getting into all the details of the metamodel.

B. Terra Nova Energy meta-model

In the previous section, a generic meta-model, at a first level of abstraction with some inheriting elements from entity such as Device has been presented.

According to the concept of inheritance, the Terra Nova Energy meta-model (Figure 2) is built as an extention and a generalization of the generic meta-model. The first element is the TNESystem that inherits from the generic element System. Other elements have been spread over two sets: generic devices and off the shelf devices that could be grouped in a specialized library.

1) Terra Nova Energy's generic devices: This first set of elements represents generic devices that allow to add a new industrial component to a library or to build a specific system. All these devices inherit from generic Device, imported from Basic package as shown in Figure 2. We present them successively as follows:

• BoxaNova models the main device of TNE systems. It collects information, sends orders and manages the entire flow of data between different devices and users. A Box-aNova device is generally created by some others devices depending on system configuration and deployment. • Router models a device that handles message transfers between different TNE elements. It handles also some other functions like controlling repeaters, sensors, actuators and some other server facilities. • Concentrator is used to model a collecting place of data coming from repeaters, sensors, actuators before sending them to routers depending on requests or preprogrammed sending tasks. • Sensor ActuatorDevice is used to capture measurement data and to send them or to act pursuant to an order. 2) Terra Nova Energy's off the shelf devices: This second set of elements represents real devices ready to use in the framework. They are another level of specialization of the first set. They may be grouped in a kind of library (blue box in Figure 2), composed with industrial devices coming from different manufacturers and in order to be integrated into TNE systems. In this paper we neither describe this library of devices nor the manner of their built because we focus on some other properties and aspects of model analysis.

In the next section, we propose a methodology of analysis showing the usefulness of these presented meta-models, and their examination and validation on real industrial cases.

IV. MODEL ANALYSIS

In this section, we present how to exploit the MDE approach to analyse placement for a given instance, and automatic placement of repeaters and concentrators for an uncompleted system, where only sensor, actuator positions are known.

The first part describes our methodological approach, based on processing properties and applied constraints for model's elements.

A. Presentation

From the TNE domain specific meta-model, presented in the previous section, an instance that describes a real system of telecontrol may be defined. According to the user needs, this instance captures a particular concern for a system and gives the necessary elements and their relevant properties. Such a model requires some processing stages to meet the final requirements in terms of analysis. These processing stages are done in response to expected system constraints applied to Fig. 3. A view of analyzer's engine model its various elements. Figure 3 shows a model overview of the various elements used to carry out those treatments for analysis methodology purpose.

The first analysis component (ModelIn) loads the instance that describes a telecontrol system according to its metamodel. This instance is composed of elements called Ele-mentIn. The Engine scans iteratively the input model, element by element, to find different properties and constraints. For each property the Analyzer checks if it fits well the associated constraint. When properties are verified, the element will be transformed by Transformation and treated by ElementOut. The ModelOut generates automatically output model formed by those transformed elements. Output model should be consistent to its meta-model that may be the input one or another depending on performed transformation.

Generally, input models have a tree shape structure and their components can be composed of other components. Properties and applied constraints can also have this composition criterion. In order to deal with this kind of structure and having a generic analyzer, possible compositions have been taken into account for analyzing and transforming tasks. In Figure 3, these composition aspects are represented by the Composed reference. The next part of this section deals with an application of this methodology to a system. First of all, analysis for Fig. 4. A screen shot of manual instance creation a completely described given instance are conducted to verify that elements are well placed and may communicate together properly. Second we tell that the Model Driven Engineering approach used, may also be helpful to build instances, starting with an uncomplete system, where only sensors and actuators are known, to add necessary repeaters and concentrators.

B. First case study with a complete TNE's instance

In this first case study, an instance of a TNE meta-model has been created with a tree editor as shown in Figure 4. It is an abstraction of a TNESystem named Telecontrol-System 1 composed of two BoxaNova that manage twenty Sensor ActuatorDevices. Each BoxaNova consists of two Concentrator elements and a Router. Communication between these components uses two Mediums; a coaxial cable that connects the first concentrator and the router and a RS232 cable that provides serial communication between the second concentrator and the router.

Sensor and actuator devices use a ZigBee protocol (IEEE802.15.4 standard) [START_REF] Baronti | Wireless sensor networks: A survey on the state of the art and the 802.15. 4 and ZigBee standards[END_REF] to communicate with the Box-aNova concentrators. To achieve this task, they use radioModules. Concentrators also have their own radio modules. In general, wireless propagation environments in TNE systems may have different attenuations. Thus, each communication between a sensor/actuator and its repeater or concentrator is defined by a medium with relevant properties and necessary (contact) interfaces. In our model, this information is kept by properties such as the attenuation value, speed of transmission, signal length, etc. Connection with the medium and other elements is ensured by the definition of two entities as ContactInterface; one for the medium and the other for the associated element.

The model also captures other information like the physical location (x, y coordinates) of the different parts of the system. In this way, Communication vs. placement analysis have been conducted. In fact, component's placement in TNE sytems is crucial to ensure good communication and proper information transfer. As the model is fully known, an analyzer can be created to verify that every Sensor ActuatorDevice element according to its coordinates may communicate to its corresponding concentrator. This verification takes the attenuation between communicating elements into account, specified as a property of the medium, which links them. Another analyzer can be defined and used for inspecting the number of sensor/actuator devices for each concentrator and its compliance to specified maxElement property.

From information analysis text file obtained (see transformation model to text in [START_REF] Czarnecki | Classification of model transformation approaches[END_REF]), manual corrections may be realized. This task is manageable for small systems but gets difficult and very expensive for large ones with many components. The next part gives an illustration of a methodology that can treat this difficulty.

C. Second case study with a partial TNE's instance

In this second case study, only the number and the position of sensor/actuator devices are given. The instance is partial, and the objective is to automatically find a solution for repeater's and concentrator's placement. This problem is a very complex one and several solutions may be found in literature [START_REF] Ghosh | Sensor network design for smart highways[END_REF][16] [START_REF] Cheng | Relay sensor placement in wireless sensor networks[END_REF].

In our case an input model is generated automatically with the required properties. As a first step, an instance of the TNE meta-model is created with a system TNESystem, a number of sensor ActuatorDevice and an environment that represents a factory hall. In our case, we suppose that radio measurements were made in the environment and have identified attenuation values.

In a second step, an analysis engine for placement, composed of three analysers, is used:

1) The primary analyzer performs an initial refinement of the input model. First, it divides the environment, specified in the input model, into zones where mediums have the same attenuation so that the range of sensoractuator radio modules, who are currently in, reach the middle. This range is calculated by applying attenuation value imposed in the input model. Then, depending on the number of sensor/actuator devices, it creates the necessary repeaters with their radio modules, and adds them to the input model with coordinates in the vicinity of the centre of the sub-environment. Finally, it creates connections (ContactInterface) between repeaters and sensor/actuator elements and adds them to the model. 2) A second pass of refinement with the same treatment approach is made by the second analyzer. Its objective is to connect repeaters with concentrators. The input model environment is split up into several zones with new radio modules ranges of repeaters. The number of added concentrators will depend on their capacity and also on the number of repeaters within reach.

3) The third analyzer performs an optimization on the obtained model that contains sensor/actuator devices that can communicate directly with concentrators without the use of a repeater. This analyzer checks for each sensor/actuator device if it reaches any concentrator and changes its connection from a repeater to a concentrator. Repeaters without connected sensor/actuator devices will be removed. At the end of these three analyzer passes, a new output model is obtained. It is composed of initial sensor/actuator devices and added components (repeaters, concentrators and connections). Figure 5 shows a simple placement layout of TNE system components obtained with the following initial model inputs: 200 sensors/actuator devices (location of these elements has been chosen randomly for the test), an environment that represents a factory with 800 meters length and 600 meters width and an attenuation value of 5 decibels by meter. Red rectangles represent repeaters, blue triangles represent concentrators and green circles represent sensor/actuator devices. The generated model allows us to have a first proposal for the positions of the repeaters and of the concentrators in the environment. This obtained placement is of course not optimal, as simple hypothesis have been used. Indeed, the presence of mobile obstacles in the environment and their influence on signal propagation should be considered. To optimize the placement and to ensure the highest data rate optimal algorithms should also be used. The orientation of elements should also be taken into account.

Nevertheless, our objective is to show that using model analysis in ubiquitous systems in MDE approach may support engineers in the design of their systems and verifying some properties in earliest stage without real implementation such as communicating object placement.

V. CONCLUSION AND FUTURE WORKS

In this paper, a generic meta-model for modeling ubiquitous telecontrol systems and, specially telecontrol systems is proposed. A variant of this metamodel is specified to fit the needs of Terra Nova Energy.

The introduction of the Model Driven Engineering approach in this field allows us to exceed the contemplative dimension of models towards productive models. The presented cases studies highlights this approach by proposing a possible model for placement analysis for an example of the TNE system, using transformations.

In the future, we would like to improve our generic metamodel in order to take the behaviour of the different components into account. We would like to move from static analysis to simulation with a complete integration of a system inside its environment.

Fig. 1 .

 1 Fig. 1. A view of generic meta-model components (Basic package)

Fig. 2 .

 2 Fig. 2. A view of Terra Nova Energy meta-model

Fig. 5 .

 5 Fig. 5. Placement layout

ACKNOWLEDGMENT

This work is supported by the city of Brest, "Brest Métropole Océane", and performed in partnership with Terra-Nova Energy. We thank them for their help.


O M Braun

Michel Peyrard
email: michel.peyrard@ens-lyon.fr

D V Stryzheus

Erio Tosatti

Collective effects at frictional interfaces

Keywords: boundary lubrication, nanotribology, viscosity, master equation, stick-slip PACS numbers: 81.40.Pq, 46.55.+d, 61.72.Hh

We discuss the role of the long-range elastic interaction between the contacts inside an inhomogeneous frictional interface. The interaction produces a characteristic elastic correlation length λc = a 2 E/kc (where a is the distance between the contacts, kc is the elastic constant of a contact, and E is the Young modulus of the sliding body), below which the slider may be considered as a rigid body. The strong inter-contact interaction leads to a narrowing of the effective threshold distribution for contact breaking and enhances the chances for an elastic instability to appear. Above the correlation length, r > λc, the interaction leads to screening of local perturbations in the interface, or to appearance of collective modes -frictional cracks propagating as solitary waves.

I. INTRODUCTION

Studies of sliding friction, a subject with great practical importance and with rich physics, attracted an increased interest during last two decades [START_REF] Persson | Sliding Friction: Physical Principles and Applications[END_REF][START_REF] Braun | Nanotribology: Microscopic mechanisms of friction[END_REF]. Tip-based experimental techniques as well as atomistic molecular dynamics (MD) computer simulations describe with considerable success the processes and mechanisms operating in atomic-scale friction. Much less is known when going to meso-and macro-scale friction, where one has to take into account that the frictional interface is inhomogeneous and generally complex. An immediate example is dry friction between rough surfaces. Even when the sliding surfaces are ideally flat but for example the substrates are not monocrystalline, or there is an interposed solid lubricant film consisting of misoriented domains, the frictional interface is again inhomogeneous. The same may be true even for liquid lubrication, if under applied load the lubricant solidifies making bridges due to Lifshitz-Slözov coalescence. In these cases, the so-called earthquake-like (EQ) type models can be successfully applied [START_REF] Braun | Nanotribology: Microscopic mechanisms of friction[END_REF][START_REF] Olami | Selforganized criticality in a continuous, nonconservative cellular automaton modeling earthquakes[END_REF][START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF][START_REF] Braun | Transition from stick-slip to smooth sliding: An earthquakelike model[END_REF][START_REF] Filippov | Friction through Dynamical Formation and Rupture of Molecular Bonds[END_REF][START_REF] Farkas | Static versus dynamic friction: the role of coherence[END_REF][START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Kinetics of stick-slip friction in boundary lubrication[END_REF][START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF][START_REF] Braun | Dependence of kinetic friction on velocity: Master equation approach[END_REF][START_REF] Braun | Kinetics and dynamics of frictional stick-slip in mesoscopic boundary lubrication[END_REF]. In the EQ model, the two (top and bottom) mutually sliding surfaces are coupled by a set of contacts, representing, e.g., asperities, patches of lubricant, or 2D crystalline domains. A contact is assumed to behave as a spring of elastic constant k c so long as its length is shorter than a critical value x s = f s /k c ; above this length the contact breaks, to be subsequently restored with lower stress. The sliding kinetics of this model may be reduced to a master equation (ME), which allows an analytical study [START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF][START_REF] Braun | Dependence of kinetic friction on velocity: Master equation approach[END_REF].

In the simplest approach, the slider is treated as a rigid body. Due to the non-rigidity of the substrates, however, several length scales naturally appear in the problem. First, different regions of the interface will exhibit different displacements. The length λ L such that for distances r ≫ λ L the displacements are independent, is known as the Larkin-Ovchinnikov length [START_REF] Larkin | Pinning in Type-II superconductors[END_REF]. It was shown [START_REF] Persson | Theory of friction: elastic coherence length and earthquake dynamics[END_REF] that for the contact of stiff rough solid surfaces λ L may reach unphysically large values ∼ 10 100 000 m. Second, deformation of the solid substrates leads to the elastic interaction between the contacts. Elasticity will correlate variations of forces on nearest contacts over some length λ c known as the elastic correlation length [START_REF] Caroli | Hysteresis and elastic interactions of microasperities in dry friction[END_REF]. Third, displacements in one region of the slider will be felt in other regions on the distance scale set by of a screening length λ s . Finally, the breaking of one contact may stimulate neighboring contacts to break too (the socalled concerted, or cascade jumps), following which an avalanche-like collective motion of different domains of the interface may appear [START_REF] Braun | Transition from stick-slip to smooth sliding: An earthquakelike model[END_REF].

In this paper we discuss collective effects in the frictional interface and propose approaches to treat them from different viewpoints. In particular, our aim is to clarify the following questions: (i) what is the law of interaction between the contacts; (ii) at which scale can the slider be considered as a rigid body, or what is the coherence distance λ c within which the motion of contacts is strongly correlated; (iii) whether the interaction effects can be incorporated in the master equation approach; (iv ) how does the interaction modifies the interface dynamics; (v ) what is the screening length λ s ; (vi) when do avalanche motion of contacts (a self-healing crack) appear, and what is the avalanche velocity?

The paper is organized as follows. The earthquakelike model, its description with the ME approach, and the elastic instability responsible for the stick-slip motion are introduced in Sec. II. The interaction between contacts is studied in Sec. III. An approach to incorporate the interaction between contacts into the ME approach in a mean-field fashion is described in Sec. IV. The role of interaction at the meso/macro-scale is considered in Sec. V. Finally, discussions in Sec. VI conclude the paper.

II. EARTHQUAKE-LIKE MODEL, MASTER EQUATION AND ELASTIC INSTABILITY

A. The earthquake-like model

In the EQ model the sliding interface is treated as a set of N contacts which deform elastically with the average rigidity k c . The ith contact connects the slider and the substrate through a spring of shear elastic constant k i . When the slider is moved, the position of each contact point changes, the contact spring elongates (or shortens) so that the slider experiences a force -F = f i from the interface, where f i = k i x i and x i (t) is the shift of the ith junction from its unstressed position. The contacts are assumed to be coupled "frictionally" to the slider. As long as the force |f i | is below a certain threshold f si , the ith contact moves together with the slider. When the force exceeds the threshold, the contact breaks and a rapid local slip takes place, during which the local stress drops. Subsequently the junction is pinned again in a less-stressed state with f bi , and the whole process repeats itself. Thus, with every contact we associate the threshold value f si and the backward value f bi , which take random values from the distributions P c (f ) and R(f ) correspondingly. When a contact is formed again (re-attached to the slider), new values for its parameters are assigned. The EQ model was studied numerically in a number of works [START_REF] Olami | Selforganized criticality in a continuous, nonconservative cellular automaton modeling earthquakes[END_REF][START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF][START_REF] Braun | Transition from stick-slip to smooth sliding: An earthquakelike model[END_REF][START_REF] Filippov | Friction through Dynamical Formation and Rupture of Molecular Bonds[END_REF][START_REF] Farkas | Static versus dynamic friction: the role of coherence[END_REF][START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Kinetics of stick-slip friction in boundary lubrication[END_REF][START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF], typically with the help of the cellular automaton numerical algorithm.

B. The master equation approach

Rather than studying the evolution of the EQ model by numerical simulation, it is possible to describe it analytically [START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF][START_REF] Braun | Dependence of kinetic friction on velocity: Master equation approach[END_REF]. Let P c (x) be the normalized probability distribution of values of the stretching thresh-olds x si at which contacts break; it is coupled with the distribution of threshold forces by the relationship P c (x) dx = P c (f ) df , i.e., the corresponding distributions are coupled by the relationship P c (x) ∝ x P c [f (x)], where f ∝ x 2 [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]. We assume that the distribution P c (x) has a dispersion ∆x s centered at x = x c .

To describe the evolution of the model, we introduce the distribution Q(x; X) of the contact stretchings x i when the sliding block is at position X. Evolution of the system is described by the integro-differential equation (known as the master equation, or the kinetic equation, or the Boltzmann equation) [START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]

∂ ∂X + ∂ ∂x + P (x) Q(x; X) = R(x) Γ(X) , (1)
where

Γ(X) = ∞ -∞ dξ P (ξ) Q(ξ; X) (2)
and

P (x) = P c (x)/J c (x) , J c (x) = ∞ x dξ P c (ξ) . (3)
Then, the friction force (the total force experiences by the slider from the interface) is given by (

k c = k i ) F (X) = N k c ∞ -∞ dx x Q(x; X) . (4)
In the steady state corresponding to smooth sliding, the ME reduces to

dQ(x)/dx + P (x) Q(x) = R(x) Γ , (5)
which has the solution

Q s (x) = N E P (x) 1 + Γ x 0 + dξ R(ξ)/E P (ξ) , ( 6
)
where N is the normalization constant,

∞ 0 dx Q s (x) = 1
, and

E P (x) = exp [-U (x)] , U (x) = x 0 dξ P (ξ) . (7)

C. Elastic instability

The solution of the ME [START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF] shows that when a rigid slider begins to move adiabatically, Ẋ > 0, it experiences from the interface a friction force F ∞ (X) < 0. Initially |F ∞ | grows roughly linearly with X, |F ∞ | ≈ K s X (here K s = N k c is the total elastic constant ("rigidity") of the interface), until it reaches a value ∼ F s -∆F s , where F s ≈ K s x c and ∆F s ≈ K s ∆x s . Gradually however contacts begin to break and reform, slowing down the increase of |F ∞ | and then inverting the slope through a displacement ∆x s until almost all contacts have been reborn. Successively the process repeats itself with a smaller amplitude until, due to increasing dispersion of breaking and reforming processes, the force asymptotically levels off and attains a position independent steady state kinetic friction value with smooth sliding.

According to Newton's third law, the external driving force F d = K(vt -X) which causes the displacement X (here K is the slider rigidity and v is the driving velocity), is compensated by the force from the interface, F d = F (X). Smooth sliding is always attained with a rigid slider. It persists for a nonrigid slider as well, so long as the pulling spring stiffness is large enough, K > K * , where

K * = max F ′ ∞ (X) , F ′ ∞ (X) ≡ dF ∞ (X)/dX . (8)
When conversely the slider, or the pulling spring elastic constant are soft enough (K < K * ) there is a mechanical instability. The driving force F d cannot be compensated by the force from the interface, and the slider motion becomes unstable at X c , where X c is the (lowest) solution of F ′ ∞ (X) = K (for details see Refs. [START_REF] Braun | Modeling friction on a mesoscale: Master equation for the earthquakelike model[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]). The mechanical instability yields stick-slip frictional motion of the slider.

Thus, the regime of motion -either stick-slip for K ≪ K * or smooth sliding for K ≫ K * -is controlled by the effective stiffness parameter K * ∼ K s x c /∆x s . When all contacts are identical, ∆x s = 0 so that K * = ∞, one always obtains a stick-slip motion.

D. Material parameters

It is useful here, before proceeding with the analytical and numerical developments necessary to answer the questions posed in the Introduction, to review the practical significance and magnitude of the model parameters.

Elastic constant of the slider. The slider (shear) elas-

tic constant K is equal to K = [E/2 (1 + σ)][L x L y /H],
where L x , L y and H are the slider dimensions, E and σ are the substrate Young modulus and Poisson ratio, respectively [START_REF] Landau | Theory of Elasticity[END_REF]. For example, for a steel slider of Young's modulus E = 2 × 10 11 N/m 2 , Poisson's ratio σ = 0.3 and the size

L x × L y × H = 1 cm × 1 cm × 1 cm, we obtain K ∼ 10 9 N/m.
Rigidity of the interface contacts. Here we characterize the typical magnitudes of the contact stretching length x c and stiffness k c . Assume the slider and the substrate to be coupled by N = L x L y /a 2 contacts, and that the contacts have a cylindrical shape of (average) radius r c with a distance a between the contacts. It is useful to introduce the dimensionless parameter γ 2 = r c /a, which may be estimated as follows [START_REF] Persson | Sliding Friction: Physical Principles and Applications[END_REF]. Consider a cube of linear size L on a table. The weight of the cube F l = ρL 3 g (ρ is the mass density and g = 9.8 m/s 2 ) must be compensated by forces from the contacts, F l = N r 2 c σ c , where σ c is the plastic yield stress. Then, γ 2 2 = (N r 2 c )/(N a 2 ) = (ρL 3 g)/(σ c L 2 ), or γ 2 = (ρLg/σ c ) 1/2 . Taking L = 1 cm, ρ = 10 g/cm 3 and σ c = 10 9 N/m 2 (steel), we obtain γ 2 ≈ 10 -3 which should be typical for a contact of rough stiff surfaces. For softer materials, and especially for a lubricated interface, the values of γ 2 would be much larger, e.g., γ 2 ∼ 0.1.

The second dimensionless parameter γ 1 = k c /Ea characterizes the stiffness of the contacts. To estimate γ 1 , assume again contacts with the shape of a cylinder of radius r c and length h (h is the thickness of the interface). Suppose in addition that one end of a contact ("column") is fixed, and a shear force f is applied to the free end. This force will lead to the displacement x = f /k c of the end, where k c = 3E c I/h 3 , E c is the Young modulus of the contact material and I = πr 4 c /4 is the moment of inertia of the cylinder [START_REF] Landau | Theory of Elasticity[END_REF]. In this way we obtain 3 . For the contact of rough surfaces, where E c = E and a > ∼ h, we have γ 0 > ∼ 1, while for lubricated interfaces where E c ≪ E, one would expect

k c = (3π/4)(E c r c )(r c /h) 3 , so that γ 1 = (3π/4) E c a 3 /Eh 3 (r c /a) 4 = γ 0 γ 4 2 with γ 0 = (3π/4)(E c /E)(a/h)
γ 0 < ∼ 1.
An estimate of characteristic values [START_REF] Persson | Sliding Friction: Physical Principles and Applications[END_REF] leads to r c ∼ (10 -3 ÷ 10 -2 ) a. Thus, for the steel slider considered above, taking r c = h = 1 µm and intercontact spacing a = 3 × 10 2 r c , we obtain N ∼ 10 3 and k c ∼ 5 × 10 5 N/m, so that the global stiffness of the interface is K s ∼ 5 × 10 8 N/m. Stick-slip versus smooth sliding. As mentioned above in Sec. II C, the regime of motion (either stick-slip or smooth sliding) is controlled by the parameter K * ∼ K s x c /∆x s . For the steel slider considered above, estimates gave K ∼ 10 9 N/m and K s ∼ 5 × 10 8 N/m. Thus, if the surfaces are rough so that ∆x s ∼ x c , then K > K * and one should typically get smooth sliding. Stick slip appears further disfavored if we consider a realistic P c (x) distribution. For all cases mentioned in Introductionthe contact of rough surfaces (both for elastic or plastic asperities), the contact of polycrystal (flat) substrates, and the case of lubricated interface, when the lubricant, melted during a slip, solidifies and forms bridges at stick, -the distribution P c (x) is rather wide with a large concentration of small-threshold contacts [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF], which makes the value of K * very small. Thus, the theory predicts that most systems do not undergo an elastic instability and should not therefore exhibit stick-slip. This conclusion contradicts everyday experience as well as careful experiments, where stick-slip is pervasive. As suggested by EQ simulations [START_REF] Braun | Kinetics and dynamics of frictional stick-slip in mesoscopic boundary lubrication[END_REF], the discrepancy is most likely caused by ignoring the elastic interaction between the contacts.

The role of interaction is considered in the next sections. First, however, we need to define the form and parameters of the interaction between contacts.

III. INTERACTION BETWEEN CONTACTS

Friction is not a simple sum of individual contact properties. The collective behavior of the contacts is important. Recently Persson [START_REF] Persson | Theory of rubber friction and contact mechanics[END_REF][START_REF] Persson | Elastic contact between randomly rough surfaces: Comparison of theory with numerical results[END_REF][START_REF] Persson | On the elastic energy and stress correlation in the contact between elastic solids with randomly rough surfaces[END_REF][START_REF] Yang | Contact mechanics: contact area and interfacial separation from small contact to full contact[END_REF][START_REF] Almqvista | Interfacial separation between elastic solids with randomly rough surfaces: comparison between theory and numerical techniques[END_REF] developed a contact mechanics theory based on the fractal structure of surfaces in order to determine the actual contact area at all length scales, which determines the friction coefficient. This approach includes the presence of multiple contacts and leads to the correct low-threshold limit P c (f → 0) = 0. Persson found that the distribution of normal stresses σ (σ > 0) at the interface may approximately be described by the expression P σ (σ) ∝ exp -(σσ) 2 /∆σ 2exp -(σ + σ) 2 /∆σ 2 , where σ is the nominal squeezing pressure, the distribution width is given by ∆σ = E * R 1/2 (E * is the combined Young modulus of the substrates, E * -1 = E -1 1 +E -1 2 where E 1,2 are the Young modula of the two substrates), and the parameter R is determined by the roughness of the contacting surfaces, R = (4π) -1 dq q 3 d 2 x h(x)h(0) e -iqx . Assuming that a local shear threshold is directly proportional to the local normal stress, f ∝ σ, we finally obtain the distribution, which is characterized by a low concentration of small shear thresholds, P c (f ) ∝ f at f → 0, and a fast decaying tail, P c (f ) ∝ exp(-f 2 /f * 2 ) at f → ∞, i.e., now the peaked structure of the distribution is much more pronounced.

However an important aspect which has to be included is the redistribution of the forces when some contacts deform or break. A concerted motion of contacts may emerge only due to interaction between the contacts which occurs through the deformation of the bulk in the directions parallel to the average contact plane. It is this aspect that we want to consider here. For the elastic interaction, a qualitative picture is presented in Fig. 1 (left). When a contact acts on the surface at r = 0 with a force f , it produces a displacement field u(r) ∝ r -1 which affects other contacts (Fig. 1a) -similar to the Coulomb potential for a point charge [START_REF] Landau | Theory of Elasticity[END_REF]. However, if there are two surfaces, then the same contact acts on the second surface with the opposite force -f and, if the two surfaces are in contact, the resulting displacement field should fall as u(r) ∝ r -3 (Fig. 1b) -similar to the dipole-dipole potential for a screened point charge near a metal surface [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF]. The question thus is the form of the interaction for the multi-contact interface (Fig. 1c). We will show that the interaction between the contacts has a crossover from the r -1 slow Coulomb decay at short distances to the faster dipole-dipole one at large distances.

A. Analytics

Let us consider an array of N elastic contacts (springs) with coordinates r i ≡ {x i , y i , 0}, i = 1, . . . , N , between the two (top and bottom) substrates. If the interface is in a stressed state, the contacts act on the top substrate with forces f i ≡ {f ix , f iy , f iz }. These contact forces produce displacements u (top) i of the (bottom) surface of the top substrate. The 3N -dimensional vectors

U (top) ≡ {u (top) i } and F t ≡ {f i } are coupled by the linear relationship U (top) = G (top) F t .
Elements of the elastic matrix G (top) (known also as the elastic Green tensor) for a semi-infinite isotropic substrate were given by Landau and Lifshitz [START_REF] Landau | Theory of Elasticity[END_REF]:

G ix,jx = g(r ij )[2(1 -σ) + 2σx 2 ij /r 2 ij ] G ix,jy = 2g(r ij ) σx ij y ij /r 2 ij G ix,jz = -g(r ij )(1 -2σ) x ij /r ij G iz,jx = -G ix,jz G iz,jz = 2g(r ij )(1 -σ) , (9)
where x ij = x ix j , g(r) = (1 + σ)/(2πEr), and σ and E are the Poisson ratio and Young modulus of the top substrate, respectively.

In the equilibrium state, the forces that act from the contacts on the bottom substrate, must be equal to F b = -F t according to Newton's third law. These forces lead to displacements of the (top) surface of the bottom substrate, U (bottom) = -G (bottom) F t . The elements of the bottom Green tensor G (bottom) are defined by the same expressions (9) except the xz elements for which

G (bottom) ix,jz = -G (top)
ix,jz (if the substrates are identical, the z displacements are irrelevant). Thus, the relative displacements at the interface due to elastic interaction between the contacts are determined by

U ≡ U (top) -U (bottom) = -GF , (10)
where

F = -F t and G = G (top) + G (bottom) .
On the other hand, the forces and displacements are coupled by the diagonal matrix (the contacts' elastic matrix) K, K iα, jβ = k iα δ ij δ αβ (α, β = x, y, z):

F = K (U 0 + U) , (11)
where U 0 defines a given stressed state (because of linearity of the elastic response, final results should not depend of U 0 ). The total force at the interface, f = i f i , must be compensated by external forces applied to the substrates, e.g., by the force f (ext) = f applied to the top surface of the top substrate if the bottom surface of the bottom substrate is fixed. Combining Eqs. ( 10) and ( 11), we obtain F = K (U 0 -GF), or

F = BKU 0 , where B = (1 + KG) -1 . ( 12
)
If one changes the contact elastic matrix, K → K+δK, then the interface forces should change as well, F → F + δF. From Eq. ( 12) we have δF = (δB)KU 0 + B(δK)U 0 . Then, δB may be found from the equation

δ[B (1 + KG)] = (δB)(1 + KG) + B(δK)G = 0.
Therefore, finally we obtain:

δF = B δK (1 -GBK)U 0 . ( 13
)
Above we have assumed that δK is small. If it is not small, we have to use the expression δF = B δK (1 -GB K)U 0 , where K = (1 + δKGB) -1 (K + δK). Now, if we remove the i * th contact by putting δk iα = -k iα δ ii * and then calculate the resulting change of forces on other contacts, we can find a response of the interface to the breaking of a single contact as a function of the distance r = r ir i * from the broken contact.

B. Numerics

Equation ( 13) may be solved numerically by standard methods of matrix algebra. We explore an idealized array of identical contacts, k iα = k c and (U 0 ) iα = u 0 δ αx for all i, organized in a square 89 × 89 lattice with spacing a = 1, with the broken contact i * at the center of the lattice. For singular terms of the Green function [START_REF] Braun | Kinetics of stick-slip friction in boundary lubrication[END_REF] we apply a cutoff at r ii = r c . Numerical results depend on two dimensionless parameters. The first is γ 1 = k c /E * a, which determines the stiffness of the array of contacts relative the substrates (here

E -1 * = E -1 top + E -1 bottom ).
The second parameter γ 2 = r c /a characterizes a single contact (or the density of asperities). For the Poisson ratio we took a typical value σ = 0.3. A typical distribution of breaking induced force changes is shown in Fig. 1

(right).

The numerical results for the x-component of dimensionless force δf = δF x /(k c u 0 ) are presented in Fig. 2. The function δf (r) exhibits a crossover from a slow Coulomb like decay δf (r) ∝ r -1 at short distances r ≪ λ c to the fast dipole-dipole like decay δf (r) ∝ r -3 at large distances r ≫ λ c . The near and far zones are separates by the elastic correlation length λ c first introduced by Caroli and Nozieres [START_REF] Caroli | Hysteresis and elastic interactions of microasperities in dry friction[END_REF]. It may be estimated in the following way: the stiffness of the "rigid block" K ∼ Eλ c should be compensated by that of the interface, K ∼ k c (λ c /a) 2 (stiffness of one contact times the number of contacts). This leads to The rigid slider corresponds to the limit E → ∞, or γ 1 → 0. Therefore, the slider may be considered as a rigid body (e.g., in MD simulation), if its size is smaller than λ c . For the steel slider considered in Sec. II D, estimation gives λ c /a ∼ 10 2 . Up to distance λ c the contacts strongly interact. If the ith contact breaks and its stretching changes on |δx i | ≈ x c , then the force on the jth contact at a distance r ij < λ c away, changes by δf j ≈ κk c a δx i /r ij , where the dimensionless parameter κ < 1 characterizes the strength of interaction (numerics gives κ ∼ 10 -3 ). In the near zone r ≪ λ c the interaction between the contacts may be accounted for within the master equation approach in a mean-field fashion as described in the next Sec. IV. At larger distances, different regions of the slider will undergo different displacements. Therefore, in the far zone, r ≫ λ c , we must take into account the elastic deformation of the slider.

λ c ≈ a/γ 1 = a 2 E/k c . (14)

IV. NEARBY CONTACTS: MEAN FIELD APPROACH

EQ model with interaction between the contacts. Let us now include the dynamical interaction between the contacts. When a contact breaks, the now unsustained shear stress must be redistributed among the neighboring contacts. We assume that because of elastic interaction between the contacts i and j, the forces acting on these contacts have to be corrected as f i → f i -∆f ij and f j → f j + ∆f ij , where ∆f ij = k ij (x jx i ) in linear approximation. For example, let at the beginning the contacts be relaxed, x j (0) = x i (0) = 0. Due to sliding motion, all stretchings grow together, so that still ∆f ij = 0. At some instant t let the jth contact break, x j (t) → 0, with the ith contact still stretched, x i (t) > 0. Clearly, as the jth contact breaks, the force on the ith

contact increases, ∆f ij (t) = -k ij x i (t) < 0.
The amplitude of interaction decreases with the distance r from the broken contact as ∆f ∝ r -1 at short distances r < λ c . Neglecting the anisotropy of interaction, we assume that

k ij = f /|r ij |,
where f is a parameter. We simulated a triangular lattice of N = 60 × 68 = 4080 contacts with periodic boundary conditions and lattice constant a = 1, with an average contact spring constant k c = 1 and radius of interaction λ c = 3a or λ c = 5a. We assumed f bi = 0 and a rectangular shape of the distribution P c (x), i.e., P c (x) = P c0 (x) = (2∆x s ) -1 for |xx s | < ∆x s and 0 otherwise, which admits an exact solution for noninteracting contacts [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF] (more realistic distributions give the same results).

Figures 3 and4 show the result of simulations for different values of the dimensionless strength of the interaction

κ = f /(k c x c ) , (15)
where x c = dx xP c0 (x) is the average stretching of the initial threshold distribution (for the rectangular distribution x c = x s ). These results yield the following conclusions. First, in the steady state, the interaction causes a narrowing of the final distribution Q s (x). At high interaction strength κ, the distribution approaches a narrow Gaussian. Second, the drop of frictional force F (X) at the onset of sliding (at X ∼ x c ) gets steeper and steeper as κ grows. Therefore, contact interactions reinforce elastic instability. Third, above a critical interaction strength, κ ≥ κ c ∼ 0.1, a multiplicity of contacts break simultaneously at the onset of sliding, and there is an avalanche, where the force F (X) drops abruptly. average avalanche size may be estimated similarly as done in Ref. [START_REF] Braun | Transition from stick-slip to smooth sliding: An earthquakelike model[END_REF].

While the full EQ model may be only studied numerically, it is always useful to have analytical results, even if only of qualitative level. In what follows we show that the main EQ results may be reproduced within the ME approach by using "effective" P c (x) and R(x) distributions defined in a mean-field fashion. In this section the ME equation is only used to reproduce the EQ results. This is however useful because it provides an additional understanding of the results as the effective distributions, obtained in this analysis, provide a description of the collective effects affecting the contacts in terms of simple functions.

Smooth sliding. Using the steady state solution of the ME, Eqs. ( 6) and [START_REF] Farkas | Static versus dynamic friction: the role of coherence[END_REF], one may approximately recover the functions P c (x) and R(x) if the stationary distribution Q s (x) is known. Indeed, for small x, where P (x) is close to zero, the left-hand side of Q s (x) allows us to find R(x) as R(x) ∝ Q ′ s (x) (see Eq. ( 5)), while the righthand side of Q s (x), where x ∼ x c and the contribution of R(x) to the shape of the steady state distribution is negligible, gives us [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]

P c (x) ∝ P (x) Q s (x) ∝ -Q ′ s (x)
. Thus, differentiating the function Q s (x) obtained in the EQ simulation, we may guess shapes of the effective distributions P c (x) and R(x) which, when substituted in the ME, would produce a solution Q s (x) close to that obtained in the EQ simulation.

Using the simulation results, let us suppose that the detached contacts form again with nonzero stretchings, i.e., that the distribution R(x) is shifted to positive stretching values,

R(x) = G(x -αx c , γx c ) , (16)
where G(x, σ) is the Gaussian distribution with zero mean and standard deviation σ,

G(x, σ) = 1 σ √ 2π exp - x 2 2σ 2 . ( 17
)
At the same time, we suppose that the effective threshold distribution P c (x) shrinks and shifts with respect to the original ("noninteracting") one,

P h (x) = βP c0 [β(x -αx c )] . (18)
Let us moreover take its convolution with the Gaussian function [START_REF] Landau | Theory of Elasticity[END_REF],

P c (x) = P h ⊗ G ≡ dξ P h (x - ξ) G(ξ, γ √ 2x c
). The results of this procedure for the rectangular distribution P c0 (x) are shown in Fig. 3. We see that with a proper choice of the parameters α, β and γ, the ME solutions Q s (x) perfectly fits the numerical EQ results (for the parameters α, β and γ in Fig. 3 Results of similar quality were also obtained for other simulated cases, e.g., for larger radius of the interaction or for wider threshold distribution P c0 (x).

The dependences of the fitting parameters α, β and γ on the dimensionless strength of interaction κ may be found in the following way. To begin with, for noninteracting contacts initially α = γ = 0 and β = 1. It is reasonable to expect that in the lowest approximation α, γ ∝ κ and β -1 ∝ κ. Indeed, because the shift of the effective distributionP c (x) appears because of the interaction, αf c = j ∆f ij , at small κ we have approximately

α ∼ 0.5 a -2 λc 0 d 2 r κx c /|r| = πκλ c x c /a 2 . ( 19
)
At large κ, however, α has to saturate, e.g., as α ∝ κ/β, because the shift cannot be larger than x c , i.e., α < 1. Then, because the distribution P c (x) shrinks from both sides, we have b 1 ∼ 2 b 2 . Thus, the interaction makes the threshold distribution P c (x) narrower by a factor β and shifts its center to the left-hand side, x c → νx c , where ν = α + β -1 changes from 1 to 0.5 as the interaction strength κ increases from zero to infinity.

Onset of sliding. The beginning of motion when started from the relaxed configuration, Q(x; 0) = δ(x), cannot be explained by the approach used above, because the effective distribution P c (x) is "self-generated" during smooth sliding, i.e., it can be applied only when the process of contacts breaking-reattachment is continuously operating. Nevertheless, the initial part of the F (X) dependence may still be described by the effective ME approach, but with the modified "forward" threshold distribution given by the expression

P ci (x) = N x ǫ0 P c0 [β 0 (x -α 0 x c )] , (20)
where N is a normalization factor, ∞ 0 dx P ci (x) = 1. The parameter α 0 is now defined so as to keep the lowest boundary unshifted, β 0 (x fix 0 -α 0 x c ) = x fix 0 with x fix 0 = x L = x s -∆x s , so that α 0 = (x fix 0 /x c )(1β -1 0 ). The "backward" distribution R(x) is still defined by Eq. ( 16) with the same parameters as above. Numerics shows that with a proper choice of the fitting parameters β 0 and ǫ 0 for a given value of κ, the initial part of the function F (X) may be reproduced with quite high accuracy. Moreover, for a rather wide range of κ values, the EQ simulation results may be reproduced by the ME approach with a reasonable accuracy using only three fitting parameter c 1 , c 2 and κ c , if the parameters β 0 and ǫ 0 in Eq. ( 20) are given by the expressions

β 0 = 1 + c 1 κ/(1 -κ/κ c ) and ǫ 0 = c 2 (β 0 -1)
, where the parameter κ c corresponds to the critical "breakdown" interaction strength when many contacts begin to break simultaneously. For κ > κ c , the drop of F (X) becomes jump-like, so that K * = ∞ and stick-slip will appear for any stiffness of the slider K < ∞. Note that the value of κ c may be estimated from the equation αx c ∼ ∆x s .

For the rectangular shape of the distribution P c0 (x) the result of this procedure is demonstrated in Fig. 4 (the fitting parameters are c 1 = 33.9, c 2 = 3.0 and κ c = 0.074).

Of course, the P ci (x) function, Eq. ( 20), can describe only the initial part of the F (X) dependence, when F (X) grows, reaches the first maximum and then decreases. To simulate the whole dependence F (X), one would have to involve the evolution of P c (x) with sliding distance, e.g., as some "aging" process P ci (x) → P c (x) (see Ref. [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]) with the initial distribution P c, ini (x) = P ci (x) and the final one P c, fin (x) = P c (x).

Stick-slip versus smooth sliding. As was mentioned above, stick-slip appears as a result of elastic instability which is controlled by the relation between the slider stiffness K and the effective interface stiffness K * . For noninteracting contacts K * ≈ K s x c /∆x s ; because typically ∆x s ∼ x c , estimates give K * < ∼ K so that stickslip should never appear. The interaction between contacts strongly enhances the elastic instability thus making stick-slip much more probable. Indeed, because of the effective shrinking of the threshold distribution, the parameter K * increases roughly as K * → K * eff ∼ β 0 K * , i.e., the effective interface stiffness K * eff grows with the strength of interaction κ, and the elastic instability can now appear. For example, for a realistic threshold distribution the dependence of K * eff on the strength of interaction κ is shown in Fig. 5.

The strength of interaction between the contacts may be found as κ ≈ κa/x c , where realistic values of the dimensionless parameter κ are of the order κ ∼ 10 -3 ; taking a ∼ (10 2 ÷ 10 3 ) r c and x c ∼ r c , we obtain κ ∼ 0.1 ÷ 1 which gives β 0 ∼ 3 ÷ 13 according to Fig. 5.

V. FAR ZONE: MESO/MACROSCALE FRICTION

At the mesoscopic scale, i.e. on distances r ≫ λ c , the substrate must be considered as deformable. Let us split the frictional area into (rigid) blocks of size λ c . In a general 3D model of the elastic slider, the nth λ c -block is characterized by a coordinate X n , and its dynamics is described by the ME for the distribution functions Q n (u n ; X n ). A solution of these MEs gives the interface forces F n (X n ). Then, the transition from the discrete numbering of blocks to a continuum interface coordinate r is trivial:

n → r, Q n (u n ; X n ) → Q[u; X(r); r], P n (u) → P (u; r), Γ n (X n ) → Γ[X(r); r], F n (X n ) → F [X(r); r]
(here r is a two-dimensional vector at the interface), and the master equation now takes the form:

∂Q[u; X(r); r] ∂X(r) + ∂Q[u; X(r); r] ∂u + P (u) Q[u; X(r); r] = δ(u) Γ[X(r); r] , (21)
where we assumed that, for the sake of simplicity, the contacts are reborn with zero stretching, R(u) = δ(u), and

Γ[X(r); r] = dξ P (ξ) Q[ξ; X(r); r] . (22)
Equations [START_REF] Yang | Contact mechanics: contact area and interfacial separation from small contact to full contact[END_REF][START_REF] Almqvista | Interfacial separation between elastic solids with randomly rough surfaces: comparison between theory and numerical techniques[END_REF] should be completed with the elastic equation of motion for the sliding body (we assume isotropic slider)

ü + η u = G 1 ∇ 2 u + G 2 ∇(∇ • u) , (23)
where u(R) is the 3D displacement vector in the slider (R = {x, y, z}), η is the intrinsic damping in the slider,

G 1 = E/2(1 + σ)ρ = c 2 t and G 2 = G 1 /(1 -2σ)ρ = c 2 l -c 2
t , E, σ and ρ are the Young modulus, Poisson ratio and mass density of the slider correspondingly, and c l (c t ) is the longitudinal (transverse) sound speed. Equation ( 23) should be solved with corresponding boundary and initial conditions. In particular, at the interface (the bottom plane of the slider, where z = 0 and {x, y} = r) we must have u x = X(r), u y = u z = 0, and the shear stress should equal F [X(r); r]/λ 2 c , where the friction force acting on the λ c -block from the interface,

F [X(r); r] = N λ k c du u Q[u; X(r); r] , (24)
should be obtained from the solution of Eq. ( 21) [here

N λ = (λ c /a) 2 ].
Equations [START_REF] Yang | Contact mechanics: contact area and interfacial separation from small contact to full contact[END_REF][START_REF] Almqvista | Interfacial separation between elastic solids with randomly rough surfaces: comparison between theory and numerical techniques[END_REF][START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF][START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF] form the complete set of equations which describes evolution of the large scale tribological system; in a general case it has to be solved numerically. However, a qualitative picture may be obtained analytically. The interface dynamics depends on whether or not the λ c -blocks undergo the elastic instability, i.e., on the ratio of the stiffness of the λ c -block K λ ≈ (2c 2 l + 3c 2 t )ρλ c [as follows from the discretized version of Eq. ( 23)] and the effective critical stiffness parameter of the interface K * λ eff = β 0 K * λ , where K * λ ∼ K λs x c /∆x s and K λs = N λ k c . If the elastic instability does not appear, then a local perturbation at the interface relaxes, spreading over an area of size λ s -the screening length considered below in Sec. V A. In the opposite case, when the elastic instability does emerge (locally), in may propagate through the interface. Below in Sec. V B we consider a simplified one-dimensional version, which allows us to get some analytical results and a rather simple simulation approach (such a model is also supported by the fact that the largest forces near the broken contact are just ahead/behind it according to Fig. 1). Recall that the interaction between the λ c -blocks is weaker than in the short-range zone, it follows the law δf ∝ r -3 which determines, e.g., the block-block interaction strength κ λ in Eq. ( 26) below (although the interaction is power-law, we may consider nearest neighbors only, because excitations at the interface, such as "kinks" introduced in Sec. V B, are localized excitations, and the role of long-range character of the interaction reduces to modification of their parameters [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF]).

A. Elastic screening length

Let us assume that the slider is split in λ c -blocks (rigid blocks) and consider the block-block interaction in a mean-field fashion (analogously to methods used in soft matter, see Refs. [START_REF] Sollich | Rheological constitutive equation for a model of soft glassy materials[END_REF][START_REF] Hébraud | Mode-Coupling Theory for the Pasty Rheology of Soft Glassy Materials[END_REF][START_REF] Bocquet | Kinetic Theory of Plastic Flow in Soft Glassy Materials[END_REF]). Due to sliding of neighboring blocks, the forces acting on contacts in the nth λ c -block get an addition shift. This effect may be accounted with the help of a substitution

f n → f n + ∆f n , ∆f n = m =n f m × Prob(m → broken) × Π mn ≈ x c m =n f m Γ m Π mn (
recall that the sum is over the λ cblocks here), or approximately

f n →   1 + x c m =n Γ m (X m ) Π mn   f n , (25)
where

Γ m (X m ) = du P m (u) Q m (u; X m ) so that N λ Γ m (X m
) is the number of broken contacts in the mth λ c -block per its unit displacement, and

Π mn ≈ N λ κ λ (λ c /r mn ) 3 (26)
describes the dimensionless (i.e., normalized on f s ) elastic interaction between the λ c -blocks separated by the distance r mn . In this way the force is given by f s Π; the numerical constant κ λ ∼ κa/λ c depends on the substrate and interface parameters.

Let us introduce the dimensionless variable ε n = x c m =n Γ m (X m ) Π mn . The shift of forces in the nth block due to broken contacts in the neighboring blocks may be accounted by a renormalization of the rate:

P n (u) → P n [(1 + ε n )u] . ( 27
)
Indeed, when contacts in the neighboring blocks break, then the forces in the given block increase, ε n > 0, and the contacts in the given block should start to break earlier, i.e., their threshold distribution effectively shifts to lower values. Making the transition from discrete sliding blocks to a continuum sliding interface, Π mn → Π(r ′r) and ε n → ε(r), we obtain a master equation of the form:

∂Q[u; X(r); r] ∂X(r) + ∂Q[u; X(r); r] ∂u + P ([1 + ε(r)]u) Q[u; X(r); r] = δ(u) Γ[X(r); r] , (28)
where we again assumed that the contacts are reborn with zero stretchings, R(u) = δ(u),

ε(r) = x c λ -2 c |r ′ -r|≥λc d 2 r ′ Γ[X(r ′ ); r ′ ] Π(r ′ -r) (29)
and

Γ[X(r); r] = dξ P ([1 + ε(r)]ξ) Q[ξ; X(r); r] . (30)
In the long-wave limit, when |dε(r)/dr| ≪ ε(r)/λ c , we may assume that the interface is locally equilibrated, i.e., the distribution of forces on contacts is close to the steady-state solution of the master equation, Q[u; X(r); r] ≈ Q s (u; r), which depends parametrically on the coordinate r through the function ε(r) entered into the expression for the rate P ([(1 + ε(r)]u). The stationary solution of the ME is known analytically [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF], and we may find the function [START_REF] Gerde | Friction and Fracture[END_REF], Γ(r) = [1 + ε(r)]/x c . Together with Eq. ( 29) this gives a self-consistent equation on the function ε(r):

ε(r) = λ -2 c |r ′ -r|≥λc d 2 r ′ [(1 + ε(r ′ )] Π(r ′ -r) . (31)
Taking into account the interaction of nearest neighboring λ c -blocks only and expanding ε(r) in Taylor series, we obtain the equation

ε(r) = Π 0 1 + ε(r) + 1 2 λ 2 c ε ′′ (r) , (32)
where Π 0 = νΠ(λ c ) = νN λ κ λ ∼ νκλ c /a and ν = 2 ÷ 4 is the number of nearest neighbors. Writing ε(r) = ε 0 + ∆ε(r), where ε 0 = Π 0 /(1 -Π 0 ), Eq. ( 32) may be rewritten as

λ 2 s ∆ε ′′ (r) = ∆ε(r) , (33)
where λ s = λ c (ε 0 /2) 1/2 is the characteristic screening length in the sliding interface.

From the known analytical steady state solution of the ME [START_REF] Braun | Dependence of kinetic friction on velocity: Master equation approach[END_REF], we may predict the dependence of screening length on temperature and sliding velocity. In particular, if T > 0, then λ s ∝ v -1/2 → ∞ as v → 0 in agreement with the results of Ref. [START_REF] Lemaître | Rate-Dependent Avalanche Size in Athermally Sheared Amorphous Solids[END_REF].

B. Frictional crack as a solitary wave

In the frictional interface, sliding begins at some weak place and then expands throughout the interface. Such a situation is close to the one known in fracture mechanics as the mode II crack, when the shear is applied along the fracture plane. In friction, a crack first opens, evolves (propagates, grows, extends) during some "delay" time τ , but then it either expands throughout the whole interface, or it will close because of the load. Below we consider the latter scenario, when one solid slips over another due to motion of the so-called self-healing crack [START_REF] Caroli | Slip pulses at a sheared frictional viscoelastic/nondeformable interface[END_REF][START_REF] Gerde | Friction and Fracture[END_REF][START_REF] Greenwood | The theory of viscoelastic crack propagation and healing[END_REF][START_REF] Greenwood | Viscoelastic crack propagation and closing with Lennard-Jones surface forces[END_REF] -a wave or "bubble" of separation moving like a crease on rug [START_REF] Vella | Statics and Inertial Dynamics of a Ruck in a Rug[END_REF]. Our plan is to adopt ideas from fracture mechanics, adapt them to the friction problem and then reduce it to the Frenkel-Kontorova (FK) model [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF] in order to describe collective motion of contacts in the frictional interface.

When one of the "collective contacts" (the λ c -block) breaks, it may initiate a chain reaction, with contacts breaking domino-like one after another. This scenario may be described accurately by reducing the system of contacts to a Frenkel-Kontorova-like model. Recall that the FK model describes a chain of harmonically interacting atoms subjected to the external periodic potential V sub (x) of the substrate. If the atoms are additionally driven by an external force f , then the equations of motion for the atomic coordinates u n take the form

mü n + mη un -g(u n+1 + u n-1 -2u n ) + V ′ sub (u i ) = f
, where m is the atomic mass, g is the strength of elastic interaction between the atoms, and η is an effective damping coefficient which describes dissipation phenomena such as the excitation of phonons etc. in the substrate. The main advantage of using the FK model is that its dynamics is well documented [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF]. Mass transport along the chain is carried by kinks (antikinks)local compressions (extensions) of otherwise commensurate structure. The kink is a well-defined topologically stable excitation (quasiparticle) characterized by an effective mass m k which depends on the kink velocity v k ,

m k = m k0 (1 -v 2 k /c 2 ) -1/2
(the relativistic Lorentz contraction of the kink width when its velocity approaches the sound speed c). Therefore, the maximal kink velocity v k max = c. In the discrete chain, kinks move in the so-called Peierls-Nabarro (PN) potential, whose amplitude is much lower than that of the primary potential V sub (x). Therefore, the kink motion is activated over these barriers, and its minimal velocity v k min is nonzero. The steady-state kink motion is determined by the energy balance: the incoming energy (because of action of the external driving force f ) should go to creation of new "surfaces" (determined by the amplitude of the substrate potential) plus excitation of phonons by the moving kink (described by the phenomenological damping coefficient η), so that

v k (f ) = f /(m k η).
FK-ME model. Thus, let us consider a chain of λ ccontacts ("atoms" of mass m = ρλ 3 c ), coupled harmonically with an elastic constant g, driven externally through a spring of elastic constant K with the end moving with a velocity v. Using the discretized version of Eq. ( 23), the elastic constants may be estimated as g ≈ 2λ c ρc 2 l and K ≈ λ c ρc 2 t . The λ c -contacts are coupled "frictionally" with the bottom substrate; the latter is described by the nonlinear force F s (u). The equation of motion of The interaction with the substrate is modeled by the function Fs(u) = kc[tanh(u) + 1.5e -u sin(3u)] with kc = 1 defined for 0 ≤ u < uc = 1 and periodically prolonged for other values of u. All contacts are driven through the springs of the elastic constant K = 0.07, their ends moving with the velocity v = 10 -4 . The motion is overdamped (m = 1, η = 100).

To initiate the breaking, two central contact interact with the substrate with smaller values of the elastic constant, kc = 0.5.

the discrete chain is

mü n +mη un -g(u n+1 +u n-1 -2u n )+F s (u n )+Ku n = f , (34)
where the driving force is given by f (t) = Kvt. The substrate force F s (u) is found from the solution of the ME for the rigid λ c -block. A typical evolution of the chain is shown in Fig. 6.

The general case may only be investigated numerically. Let us first consider a simplified case, when F s (u) has the sawtooth shape, i.e. it is defined as

F s (u) = k c u for 0 ≤ u < u c (35)
and periodically prolonged for other values of u. We assume that f is approximately constant during kink motion (otherwise, the kink will accelerate during its motion along the chain); this is correct if the change of the driving force ∆f = Kv ∆t during kink motion through the chain, ∆t = L/v k (L is the chain length and v k is kink velocity), is much lower than k c u c , or

K/k c ≪ (v k /v)(u c /L).
Let us define the function F (u) = F s (u)+Ku-f . The degenerate ground states of the chain are determined by the equation F (u) = 0. Let the right-hand side (n → ∞) of the chain be unrelaxed,

k c u R + Ku R = f , or u R = f /(k c + K) , (36)
while the left-hand side (n → -∞) already undergone relaxation, k c (u Lu c ) + Ku L = f , or

u L = (f + k c u c )/(k c + K) . (37)
Thus, the FK-like model of friction (the FK-ME model) is described by Eqs. [START_REF] Braun | Crack in the frictional interface as a solitary wave[END_REF] and [START_REF] Burridge | Model and theoretical seismicity[END_REF] with the boundary conditions given by Eqs. [START_REF] Schmittbuhl | Propagative Macrodislocation Modes in an Earthquake Fault Model[END_REF] and [START_REF] Kierfeld | Slow Crack Propagation in Heterogeneous Materials[END_REF].

Continuum-limit approximation. Let the system be overdamped (ü = 0); later on we shall remove this restriction. In the continuum-limit approximation, n → x = na (a = 1), the motion equation takes the form

mηu t -a 2 gu xx + F (u) = 0, F (u)| x→±∞ = 0 . (38)
We look for a solution in the form of a wave of stationary profile (the solitary wave), u(x, t) = u(xv k t), so that u t = -v k u ′ and u xx = u ′′ . In this case Eq. [START_REF] Braun | Hysteresis in the underdamped driven Frenkel-Kontorova model[END_REF] takes the form

mηv k u ′ + a 2 gu ′′ = F (u) , (39)
which may be solved analytically by standard methods [START_REF] Braun | Crack in the frictional interface as a solitary wave[END_REF].

A solution of Eq. ( 39) with these boundary conditions exists only for a certain value of the kink velocity v k , defined by the equation

(mηv k ) 2 = ga 2 (k c + K)(2 -β) 2 /(β -1) , (40)
where β = k c /(k * -K) and k * = f /u c . The solitarywave solution exists for forces f min < f < f max only.

The minimal force which supports the kink motionthe Griffith threshold -is given by

f min = 1 2 k c + K u c . (41)
The maximal force, for which a kink may exist, is given by

f max = (k c + K) u c ; (42)
at higher forces, the barriers of F s (u) are degraded, the stationary ground states disappear, and the whole chain must switch to the sliding state. From Eq. ( 40) we can find the kink velocity as a function of the driving force. At low velocities

v k ≈ (f -f min )/m k η , (43)
where we introduced the effective kink (crack) mass

m k = m 4a u c g k c 1 + K k c , (44)
while at f → f max the velocity tends to infinity,

mηv k ≈ gk c (k c + K)a 2 u c (f max -f ) . ( 45
)
The latter limit should be corrected by taking into account inertia effects. The term mü in Eq. ( 34) gives mv 2 k u ′′ for the solitary-wave solution, so it can be incorporated if we substitute in the above equations g →

g eff = g(1 -v 2 k /c 2 0 )
, where c 0 = (ga 2 /m) 1/2 is the sound speed along the chain. The high-velocity limit now takes the form

v k ≈ c 0 1 + mη 2 (f max -f ) k c (k c + K)u c . (46)
Simulations. The continuum-limit approximate is accurate for the case of strong interaction between the contacts, g ≫ 1; in the opposite limit one has to resort to computer simulation. We solved Eq. ( 34) by the Runge-Kutta method. As the initial state, we took the chain of length N (typically N = 3 × 10 3 or 3 × 10 4 ) with periodic boundary conditions and all contacts relaxed, but the threshold breaking value for two central contacts was taken lower than for the other contacts. Then the driving force increases because of stage motion, two central contacts break first and initiate two solitary waves of subsequent contact breaking which propagate in the opposite directions through the chain. The value k ′ c of the lower threshold of the central contacts determines the driving force and therefore the kink velocity; the lower this threshold, the lower the threshold force for the motion to start [START_REF] Braun | Crack in the frictional interface as a solitary wave[END_REF]. As soon as the kink motion is initiated, the k c -values of the central contacts are restored to the same value as for other contacts (otherwise these contacts will act as a source of creation of new pairs of kinks), and we begin to move the stage in the opposite direction, v > 0 → v b < 0, so that the driving force linearly decreases with time (see Fig. 7b), the average chain velocity ui = N -1 i ui decreases as well (Fig. 7c) until the motion stops (Fig. 7a). Also, such an algorithm allows us to find the dependence of the kink velocity determined as

v k = n -1 k N ( ui -v) , (47)
where n k = 2 is the number of moving kinks in the chain and v = uL,R = v b K/(k c +K) is the background velocity, on the driving force f . These dependences are presented in Fig. 8; they agree well with that predicted by Eqs. ( 40) and [START_REF] Guozden | Supersonic Crack Propagation in a Class of Lattice Models of Mode III Brittle Fracture[END_REF]. Contrary to the continuum-limit approximation, in the discrete chain of contacts the kink oscillates during motion (see Fig. 7d) -the well-known discreteness effect of the FK model due to existence of the PN barriers f PN . The stronger the elastic interaction between the contacts, the larger the kink "width" and the smaller the kink oscillations (compare Figs. 8a and8b). The amplitude of oscillations also depends on the shape of the "substrate potential" [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF] -it is larger for a sawtooth potential F s (u), but smaller for a smoother shapes. Recall that the λ ccontacts are characterized by a smooth dependence F s (u) as follows from the master equation. The PN oscillations determine the lowest average kink velocity. Therefore, the lowest velocity allowed for the frictional crack propagation, v k min , is determined by the parameters g and λ c -the larger are g and λ c , the smaller is v k min .

Discussion. The FK-ME model used here is rather close to the well-known 1D Burridge-Knopoff (BK) model of earthquakes with a velocity-weakening friction law [START_REF] Burridge | Model and theoretical seismicity[END_REF]. The difference is in the interface force F s (u): we use the function derived from the ME-EQ model (with well-defined parameters which may be extracted from experiments or calculated from first principles), whereas the BK model adopts a phenomenological velocity-dependent function for F s . Nevertheless, the qualitative behavior of the two models is similar, the BK model also exhibits solitary-wave dynamics as was demonstrated numerically in Ref. [START_REF] Schmittbuhl | Propagative Macrodislocation Modes in an Earthquake Fault Model[END_REF]. In our case, however, by reducing the model to the FK-ME one, we can describe the solitary waves analytically and rigorously.

In the simulation we started from the well-defined initial configuration, when all contacts are relaxed except the one or two where kink's motion is initiated. If one starts from a random initial configuration, we expect that kinks will emerge at random places, so that several kinks may propagate through the system simultaneously, as was observed in simulation of the BK model [START_REF] Schmittbuhl | Propagative Macrodislocation Modes in an Earthquake Fault Model[END_REF].

Also we assumed that all λ c -contacts are characterized by the same F s (u) dependence and thus have the same threshold values F th . This is correct if the number of original contacts within a single λ c -contact, N λ = (λ c /a c ) 2 , is infinite. Otherwise, different λ c -contacts will have different threshold values F n however the distribution of their thresholds is narrower that the distribution of thresholds of single asperities by a factor √ N λ . A narrow distribution of thresholds will nevertheless have a qualitative effect because rupture fronts may stop when they meet λ ccontacts with a threshold above the driving force. When the interface is disordered, the avalanches will have finite lengths and may become short for forces near f ini , for which the rupture fronts propagate at the minimal velocity.

Our approach may also incorporate the existence of disorder and defects always present in real materials. On the one hand, defects may nucleate kinks (cracks); on the other hand, the kink propagation may be slowed down up to its complete arrest due to pinning by the defects. For example, the slowing down of the 1D crack propagating through a 2D system with quenched randomly distributed defects was considered in Ref. [START_REF] Kierfeld | Slow Crack Propagation in Heterogeneous Materials[END_REF].

Thus, reducing the EQ-ME model of friction to the FK-ME one, we described avalanche-like dynamics of the frictional interface -the solitary wave of contacts breaking. If the force F s (u) has a sawtooth shape, then the interface dynamics may be described analytically; otherwise one has to use numerics. The analogy with the FK model may be extended even further:

• The driven FK model exhibits hysteresis when the force increases and then decreases [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF][START_REF] Braun | Hysteresis in the underdamped driven Frenkel-Kontorova model[END_REF]. The same effect was observed in the large-scale crack simulation [START_REF] Holland | Ideal Brittle Fracture of Silicon Studied with Molecular Dynamics[END_REF], thus could be observed in the frictional interface too.

• Effects of nonzero temperature may be considered. One may predict that at T > 0 the sliding kinks will experience an additional damping, while the immobile (e.g., arrested) kinks will slowly move (creep) due to thermally activated jumps.

• As shown in Refs. [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF][START_REF] Braun | Driven kink in the Frenkel-Kontorova model[END_REF], a fast driven kink begins to oscillate due to excitation of its shape mode, and then, with the further increase of driving, the kink is destroyed. This effect is similar to what is observed in fracture mechanics, where cracks begin to oscillate and then branch [START_REF] Fineberg | Instability in dynamics fracture[END_REF].

• If the interaction between the atoms is nonlinear and stiff enough, the FK model admits the existence of supersonic kinks [START_REF] Braun | Supersonic and multiple topological excitations in the driven Frenkel-Kontorova model with exponential interaction[END_REF] which are similar to solitons of the Toda chain. It would be interesting to study if similar waves may appear in the frictional interface, as was predicted in crack propagation [START_REF] Guozden | Supersonic Crack Propagation in a Class of Lattice Models of Mode III Brittle Fracture[END_REF].

• One may suppose that the damping coefficient η in the equation of motion [START_REF] Braun | Crack in the frictional interface as a solitary wave[END_REF] depends on the kink velocity, η(v). In fracture mechanics, this coefficient defines the rate at which the energy is removed from the crack edge, thus it plays a crucial role.

• A large number of works is devoted to different generalizations of the FK model to 2D system (e.g., see [START_REF] Braun | The Frenkel-Kontorova Model: Concepts, Methods, and Applications[END_REF]). For example, if kinks attract one another in the y (transverse) direction, they unite into a line (dislocation) which moves as a whole (or due to secondary kinks).

VI. CONCLUSION

We discussed the crucial role in sliding friction of the elastic interaction between the contacts at the inhomogeneous frictional interface and proposed various approaches to treat this problem from different viewpoints. The interaction produces a characteristic elastic correlation length λ c = a 2 E/k c . At distance r < λ c the slider may be considered as a rigid body but with a strong contacts' interaction, which leads to shrinking of the effective contact breaking threshold distribution and an enhanced possibility for a mechanical elastic instability to appear, which is conducive to stick slip. At large distances r > λ c , the contact-contact interaction leads to screening of local perturbations in the interface, or to appearance of collective modes (frictional cracks) propagating as solitary waves.

In our work we assumed that the external stress (the driving force) is uniform across the system. In a general case, however, stress is nonuniform and may moreover change with (adjust itself to) interface dynamics, so that the problem should be considered self-consistently. For given boundary conditions, determined by the experimental setup, one should calculate the stress field, e.g., by finite element technique, which provides the driving force f (r) in the FK-ME model. The latter defines the displacement field at the interface through the solution of the FK-ME master equations. The displacement field in turn is to be used as the boundary condition for the elastic-theory equations at the frictional interface (from other sides of the slider, the boundary conditions should correspond to a given experimental setup).

~r 3 -FIG. 1 .

 31 FIG. 1. Left: decaying of the displacement field at the interface (schematic): (a) for a single contact u(r) ∝ r -1 , (b) for a single hole u(r) ∝ r -3 , and (c) for the array of contacts. Right: change of forces on contacts when the central contact is removed (γ1 = 0.06).

3 FIG. 2 .

 32 FIG. 2. (Color online): Dependence of the change of forces δf (r) on the distance x from the broken contact for three values of the interface stiffness: γ1 = 0.003 (blue down triangles, dashed line), 0.06 (red solid circles, dotted line) and 0.8 (black up triangles, solid line) at fixed value of γ2 = 0.3 (σ = 0.3). The lines show the corresponding power laws.

FIG. 4 .

 4 FIG. 3. (color online):The steady state distribution Qs(x) for the rectangular threshold distribution Pc0(x) with xs = 1 and ∆xs = 0.25 and different values of the interaction strength κ = 0, 0.02, 0.06, 0.1, and 0.5. The EQ simulations (dotted) are compared with the ME results (solid curves).

  we used expressions β = 1 + b 1 κ, α = b 2 κ/β and γ = b 3 αb 4 α 2 with the coefficients b 1 = 18, b 2 = 9.6, b 3 = 0.142 and b 4 = 0.232).

FIG. 5 .

 5 FIG. 5. The effective interface stiffness K * eff (normalized on the noninteracting value) as a function of the strength of interaction κ for the realistic threshold distribution Pc0(x) = (2/xs) u 3 e -u 2 , u ≡ x/xs with xs = 1, when K * /Ks = 0.179.

FIG. 6 .

 6 FIG. 6. (color online): Color map of atomic velocities for a typical evolution of the chain of contacts. The nearest neighboring contacts interact elastically with the constant g = 25.The interaction with the substrate is modeled by the function Fs(u) = kc[tanh(u) + 1.5e -u sin(3u)] with kc = 1 defined for 0 ≤ u < uc = 1 and periodically prolonged for other values of u. All contacts are driven through the springs of the elastic constant K = 0.07, their ends moving with the velocity v = 10 -4 . The motion is overdamped (m = 1, η = 100). To initiate the breaking, two central contact interact with the substrate with smaller values of the elastic constant, kc = 0.5.

FIG. 7 .

 7 FIG. 7. (color online):Evolution of the chain of N = 3000 contacts. The nearest neighboring contacts interact elastically with the constant g = 25, the interaction with the substrate is modeled by the sawtooth function[START_REF] Burridge | Model and theoretical seismicity[END_REF] with kc = 1 and uc = 1. All contacts are driven through the springs of the elastic constant K = 0.07, their ends moving with the velocity v = 10 -4 . The motion is overdamped (m = 1, η = 100). To initiate the breaking, two central contacts interact with the substrate with smaller spring constants, k ′ c = 0.5. When the kinks motion begins, the elastic constants of the central contacts restore their values to kc = 1, and the driving velocity changes its sign, v → v b = -2 × 10 -4 . (a) shows the kinks centers (defined as places where the atomic velocity is maximal), (b) shows the driving force f (t), (c) shows the average chain velocity ui = N -1 i ui, and (d) demonstrates oscillation of the velocity due to PN barriers.

FIG. 8 .

 8 FIG. 8. (color online): Kink velocity versus the driving force for (a) g = 5 (v b = -4 × 10 -5 ) and (b) g = 25 (v b = -2 × 10 -4 ); N = 3 × 10 4 , other parameters as in Fig. 7. Blue solid and red dashed lines correspond to Eqs. (40) and (43), correspondingly.

ACKNOWLEDGMENTS

We wish to express our gratitude to E.A. Jagla, B.N.J. Persson, M. Urbakh, and S. Zapperi for helpful discussions. This work was supported in part by CNRS-Ukraine PICS grant No. 5421, by ESF Eurocore FANAS AFRI through CNR-Italy, by PRIN/COFIN 20087NX9Y7, and by the SNF Sinergia Project NPA1617. O.B. acknowledges hospitality at SISSA and ICTP Trieste.