Analysis of Chaotic Response of Frenkel-Kontorova-Tomlinson Model

: The Frenkel-Kontorova-Tomlinson (FKT) model represents mechanical systems in which the atomic smooth surfaces of two bodies slide against each other. The model is very sensitive to changes of the system parameters, and ranges from simple stable harmonic to chaotic solutions. The design of the model between two bodies for the dynamic problem, following the network method rules, is explained with precision and run on standard electrical circuit simulation software. It provides the phase diagrams of atom displacement for each atom and the total friction force by the summation of all the atom displacements. This article is focused on studying the effect of the selected time step on the result and in the lack of sensitivity of Lyapunov exponents to assess chaotic behaviour.


Introduction
The physical phenomenon of friction between sliding bodies is a relevant topic in mechanical engineering [1,2]. In the case of dry friction and atomic smooth surfaces, different and complex mechanisms have been considered in the literature. In dynamic friction, the aim is to know the friction force, which is strongly dependent on the relative sliding body velocity.
Historically, the first attempt to apply Coulomb's laws at atomic level is attributed to Tomlinson [3] considering the atoms of one surface as uncoupled oscillators. The coupling between atoms was later considered by the Frenkel-Kontorova model [4,5].
The combination of the former models provides the widely used Frenkel-Kontorova-Tomlinson (FKT, hereinafter) [6,7] model. With regard to dynamic friction, the solutions ranging from periodic solutions or limit cycles to chaotic orbits or bound non-periodic solutions.
A short summary of the used of FKT models could start with the FKT version proposed by Elmer and Weiss (1995) [6,7], which was improved by Gyalog and Thomas (1997), who developed a 2-D version of this model in order to represent the motion of an atomic force microscope tip over a surface [8]. Later, Elmer (1997) used a sticking-dependent static friction force and a velocity-dependent kinetic friction force [9]. Hölscher and Zwörner (1998) claimed this kind of models in the analysis of this type of microscopes [10]. Gnecco et al. (2000) analysed the friction in these devices for NaCl (100) [11]. Sang et al. (2001) corrected the obtained relation in the former study at not very low velocities, proposing new expression for the relation between the scanning velocity and the friction

The Governing Equations
The physical scheme of the FKT model, Figure 1, assumes two rigid sliding bodies. The interaction between the atoms of the surface layer of the upper sliding body and the neighbouring atoms and those in the bulk are represented by springs. The hypotheses currently assumed in the FKT model are summarised in the literature [11]. One of these considers that the interaction potential between the atoms of a fixed surface and the atoms of the surface, which has a relative displacement from the other one, can be modelled as a harmonic function. Considering an atom of the moving surface, the following equation can be written from the balance of the forces in this element [15] where ξ i is the displacement of the particle i from its equilibrium position, m is the mass of the particle, κ 1 and κ 2 is the stiffness between neighbouring atoms on the moving surface and substrate in the moving surface, respectively, v B is the body velocity, c is the ratio between the atoms number on each surface, β is the viscous damping associated with the moving surface, γ is the viscous damping associated with the still surface, and V v B , ξ j , t is the potential function. As mentioned above, the potential function is represented by: where b is the amplitude of potential function.
Each equation system in the FKT model is similar to the equation of an oscillatory system, which is described by where k represents the stiffness of the spring elements, c is the damping coefficient, m represents the mass, and u represents displacement variable. The function f can be a trigonometric function or represent the non-linear load-displacement curve of a spring, Figure 2. Let us assume that we have a system without energy dissipation and that the restoring force can be approximated by Thus, where ω 2 0 = k/m and µ represents dissipation factor. Equation (5) is called Duffing's equation [12,15,[20][21][22]. Another variant of this equation is: where a 1 , a 2 , a 3 , and a 4 represent coefficients, and ω is the angular frequency in the harmonic function. Therefore, Equation (6) will also be analysed in this work to compare the effects of the procedures that will be applied in the FKT model equations system.

The Convergence Criterion and the Lyapunov Exponent
As mentioned, we apply a software to the system of equations. Kirchhoff's laws for the circuit are applying at time steps selected in order to obtain a stable convergence and that numerical approximations of integrations are sufficiently accurate [13].
When the algorithm has reached convergence for any iteration k, we have |v k+1 where v n is the voltage in the node n, and These two parameters, RELTOL and VNTOL, which are the relative tolerance and absolute convergence, respectively, are related with the time step.
The ITL4 sets an upper limit on the number of iterations allowed to obtain a convergent solution at a timepoint.
The algorithm employed in this paper for determining Lyapunov exponents was proposed by Wolf et al. [14]. The algorithms allow the estimation of non-negative Lyapunov exponents from an time series. In this reference, the method is tested on model systems with known Lyapunov spectra and applied to data for the Belousov-Zhabotinskii reaction and Couette-Taylor flow [14].
The open-source free-software called Ngspice, which is used in this work, provides code models to support numerical simulation through a fast event-driven algorithm. This software is used for the simulation of large and complex circuits. In this work, we take advantage of its performance. Moreover, several tests were developed with equivalent MATLAB functions with acceptable results in case of simpler models, only with three atoms.

Simulation and Results
The model used in this work is based on the formal equivalence between the finite difference equations of the mathematical and network models.
Once the equivalence between electric and mechanical variables has been chosen, linear terms of the PDE are easily implemented by linear electrical devices, such as resistors, capacitors, and coils, while non-linear and coupled terms are implemented using auxiliary circuits or controlled current and voltage sources. The last are a special kind of source whose output can be defined-by software-as a function of the dependent or independent variables defined in any node or any component of the model. In addition, boundary and initial conditions-linear or not-are also immediately implemented by suitable electric components.
Once the network model has been designed, it is run with no need for other mathematical manipulations since the simulation code does this. The well-tried and powerful software for circuit simulation Ngspice, requires relatively short computing time thanks to the continuous adjusting of the internal time step required for the convergence. Solution simultaneously provides all the variables of interest: relative displacement and velocity of each atom as a function of time, phase diagrams, etc.
The Duffing equation will be solved for the following parameters: a 1 equal to −1, a 2 equal to 1, a 3 equal to 0.7, a 4 equal to 0.7, and ω equal to 1.4. The initial values are zero for the displacement and the velocity. For the calculation, the ITL4 has been set to 4000 iterations and DTMIN, the minimum internal time step, is set up in order to reduce the calculation time.
We solved the equation numerically for different convergence tolerances. Figure 3 shows the phase diagram for RELTOL and VNTOL equal to 10 −10 . The simulated time is 1000 s. Figure 4 represents only the interval between 900 and 1000 s, while Figure 5 represents only the interval between 0 and 100 s. Figure 6 shows the Lyapunov exponents of this system, which are negative and therefore prove the system stability. Nevertheless, it is necessary to take care of, since it is not immediate that in physical phenomena described by differential or difference equations, the existence of a trajectory with a positive Lyapunov exponent implies the system is instable and by opposite a trajectory with negative Lyapunov exponent is stable. As Balibrea and Caballero showed [23], the positivity or negativity are held uniformly in subsets of the phase space. In such cases positivity implies instability and negativity implies stability. A detailed analysis of Lyapunov exponents is applied to some examples in certain articles [23,24].   At this moment, we seek to analyse the effect of the convergence tolerance in the solution in the first seconds. Thus, Figure 7 represents the solution for RELTOL and VNTOL equal to 10 −15 . There are differences between Figures 5 and 7 even when using these small values for tolerances. It is also possible to find differences in the case of a courser convergence tolerance as in Figure 8, with RELTOL and VNTOL equal to 10 −6 . Finally, with RELTOL and VNTOL equal to 10 −20 , Figure 9, we have found the same solution as in Figure 7. Since the number of calculated points is very high with this convergence tolerance, only a selection of a few points corresponding to each second is recorded separately in order to obtain the last figure; and the calculation then starts again from zero to obtain the following second recorded.    After the previous analysis, the FKT model, a more complex case, will be solved for the following parameters: v B equal to 0.5, c equal to 144/233, b/m equal to 0.1, κ 1 /m equal to 1, κ 2 /m equal to 1.5, γ/m is equal to zero, and β/m is equal to 0.1. We solved the model numerically for different convergence tolerances. Figure 10 shows the phase diagram for RELTOL and VNTOL equal to 10 −6 . The represented time is range from 1450 to 1500 s. Figures 11-13 show the effect of reducing the value of these tolerances. Figures 14 and 15 show the results for the same tolerances as Figure 13 in different time intervals.         An additional calculation, represented in Figure 18, has been developed using different characterisation parameters for the element of the network. Thus, the initial estimation for zero and infinite resistance, −1e −8 and 1e +8 , respectively, has been change to −1e −20 and 1e +20 . The differences between Figures 10 and 18 shows the difficulty to obtain the real behaviour of this system. At this moment we seek to distinguish the transition between limit cycles and chaos by changing the value of k 2 to 1.4. Figures 19-24 show the effect of the convergence tolerance.          It is very difficult for the hardware available to obtain the Lyapunov exponent for this equations system during at least 300 s. We try to study the phenomenon with a reduced model with 21 equations. In this case, the variable c is 13/21. Figures 28 and 29 show the effect of κ 2 .   Figure 29, with a lower tolerance than Figure 30, is the most accurate. The phase diagram shows that the system has no determined trajectory.
It is possible to calculate the Lyapunov exponent for this equations system, Figure 31. There is no positive exponent; however, the system is chaotic.  In addition to the former calculations, a comparison is made between the results using Gear Method with MATLAB and Network Simulation methods with ngspice for FKT problem with only 3 atoms, Figures 32-37. Figures 32 and 33 represent the phase diagrams with the parameter κ 2 equal to 1.6 using MATLAB and ngspice, respectively. For this parameter, the system behaviour is periodical, as both programs show. Furthermore, the results in both programs are the same. Figures 34 and 35 represent the phase diagrams with the parameter κ 2 equal to 1.4 using MATLAB and ngspice, respectively. For this parameter, the system behaviour is chaotic, as both programs show. Furthermore, the results in both programs are similar but not equal. The Lyapunov exponents for the two values of parameter κ 2 are verified in the results shown in Figures 36 and 37.      The Lyapunov exponents for two values of parameter κ 2 and 233 atoms, considering only 10 s, verify these results, Figure 38. In this case, the Lyapunov exponents are not able to display the chaotic behaviour. Furthermore, Figure 39 shows a bifurcation diagrams and graph with large Lyapunov exponent (LLE) for κ 2 = 1.5 [25][26][27].

Final Comments and Conclusions
The non-linear mechanical system of equations of FKT is solved by an efficient and reliable model which was designed using the network method. From the solutions, the time step effect has strong influence and is able to remove the chaotic behaviour of the system. In this case, all the Lyapunov exponents are negative.
It is an extended practice particularly in experimental and applied dynamics, to associate orbits with positive Lyapunov exponents to their instability or chaotic behaviour and negative Lyapunov exponents with stability or non-chaoticity. However, such statement has no firm mathematical foundation if some restrictions on the maps describing the dynamical systems are not introduced. In fact, in Reference [23,24], several examples of system have exactly the opposite behaviour. We have systems with orbits in which Lyapunov exponents are all negatives but being chaotic in a neighbourhood of the orbit and others having positive Lyapunov exponents but non-chaotic in the neighbourhood. In Figure 13, we graphically show another example of such contradictory behaviour.
Author Contributions: J.S.R. Dry friction, numerical simulation, network simulation method. F.B.G. Dynamical system, chaos, partial differential equations. J.A.M.N. Numerical simulation, network simulation method, surface technology. F.M.G. Network simulation method. All authors have read and agreed to the published version of the manuscript.