Simulating the friction between atomic layers by using a two- q model: Analysis of the relative motion and coherence

Studying friction between atomic layers is not only of great interest for the fundamental aspect of the tribology but also important for many applications such as the layer adhesion in wearable technologies and energy saving. The previous theoretical study has used the modified Prandtl-Tomlinson model to describe the motion of the tip above a two-dimensional atomic layer in an atomic force microscopy experiment. Here the degree of freedom for the substrate has been further explicitly included in the simulation, which is significant because the coherence between the sensing and the substrate layers can be explicitly addressed by computing their relative motion. For both layers, graphene has been chosen as an example for the simulations. Based on the simulations reported here, which agree with the previous relevant theoretical and atomic-force-microscopy experimental results, the motions between the sensing sheet and the substrate can be clearly distinguished. The dependence of motion and force on the parameters for the mechanical properties of the individual layers and the interaction potential between the layers has been carefully studied. For the relatively large values of the parameters for the me- chanical properties, the relative motions between the sensing sheet and the substrate show that there would be coherence between the layers, which is beneficial for the adhesion between them. However, many other parameter spaces can be studied further in the future. Similar to the simulations of the motions of the atomic layers, the computed force of the atomic-force-microscopy tip can also indicate the stability of the layers. The theoretical work reported can be used to identify explicitly the relative motions between the sensing sheet and the substrate, providing a substantial improvement for the understanding of the friction between atomic layers. Moreover, in principles, the modeling methodology proposed can be generalized to describe any number of layers in the thin-film devices, by adding a q -parameter for each layer.


Introduction
Wearable technology (WT) has attracted much attention recently owing to the rapid development of thin-film technologies and organic electronics [1]. In WT, multimedia, sensors, and wireless communications can be integrated into clothes, supporting human gestures, eye movements and the other interactive ways to help us realize functionalities such as the real-time monitoring of a human body. In particular, sensors are the important components for WT to collect and process information [2]. At present, three kinds of sensors have been developed for WT, based on fibers, composites and thin films, respectively, such as ZnO-based light sensors [3]. Among these, the thin-film sensors have many advantages, including light weight, great mechanical flexibility, and high sensitivity [4][5][6]. The thin films deposited onto the flexible substrate can also withstand pressures, bending and vibrations with a minimum impact on their performances, making them outstanding in the field of the smart WT [7]. Most of today's common sensors can respond to piezoelectric, piezoresistive, and capacitive pressures [8]. The mechanical couplings between the sensing sheet and the substrate have a vital influence on the performance of the sensing sheets. Among these mechanical couplings, frictions (or adhesions) are one of the most important interactions between the sensing sheet and the substrate [9]. To achieve a great device performance, the understanding of the adhesion between the sensing sheet and the substrate is crucial to the design of the robust sensing device functionality.
Friction is an indispensable counterforce against the relative motion, which fundamentally stems from electromagnetic interactions between the particles in the touching surfaces at the atomic scale [10]. Approximately one quarter of the world's energy is lost or dissipated to the thermal energy due to frictions. Finding an effective method to reduce the energy waste owing to frictions is still very challenging. The two-dimension layered materials, such as graphene, have exhibited excellent properties for lubrications and even super-lubrications [11]. Atomic force microscopy (AFM) experiments have recently suggested that the friction force was strongly dependent on the number of layers [12]. Especially, the friction between the sensing sheet and the substrate is interesting, requiring a proper understanding at the atomic scale. The finite element (FE) simulation has become a very effective method to understand the time evolution of the motions and the frictions at the atomic scale, implemented in the numerical solver for the differential equation derived by using Newton's laws of motion [13,14]. Moreover, the FE analysis has been widely applied to many research fields, which regards the solution domain as composed of many finite elements, finds approximate solutions, and then deduces the solutions of the problems by interpolations. The FE method not only has high computational accuracy, but also can be adapted to various complicated situations [15] such as solving differential equations with many variables.
The interactions between the AFM tips and the atomic layers have been simulated previously based on FE methods, by introducing a parameter q to represent the motion of the sensing sheets that could be bending, stretching, and shearing, etc [16]. Therein, the frictional process has been simulated based on the Newton's laws of motion and the reasonable assumption of the potential energy of the system. The Prandtl-Tomlinson (P-T) model has been widely used to describe the motion of the AFM tip at the atomic scale on a surface with a periodic potential with an amplitude V 0 . In the P-T model, a parameter η = 4π 2 V0 ka 2 (k is the spring constant for the tip cantilever and a is the periodicity of the surface lattice) is introduced to determine whether there is a continuous sliding (η < 1) or an abrupt jump (η > 1). The former region is called super lubricity region while the latter the stick-slip region. Andersson, et al., have shown that the frictional forces between two atomic layers could be simulated properly by using a simple yet fundamental model, the modified P-T model (the potential of the sensing sheet and the kinetic energy of the tip have been included explicitly) [17], compared with the relevant experiment [12]. Therein the motions of the layers and the forces of a tip (dragged over an atomically thin layer that lies on a substrate) can be simulated by solving the differential equations following the Newton's second law [16]. The modified P-T model [17] has also been improved in this work by introducing the lattice mismatching between the sensing sheet and the substrate. However, it would be good if the motions of the substrate can also be clearly described separately apart from the potential terms. It would be appropriate to introduce an additional degree of freedom (another q parameter) for the substrate such that the motions of the sensing sheet and the substrate can be described explicitly.
In this work, the previous computational results in Ref. [16] has been benchmarked. Then the potential energies for the substrate and the interaction energies between the sensing sheet and the substrate has been designed. An additional degree of freedom has been introduced to explicitly describe the motion of the substrate [16] such that we can have a comprehensive picture for the frictions between the sensing sheet and the substrate by computing the relative motions. Note that here the simulation of the motions and forces has been carried out for the thin films that are already formed during the AFM experiment after the fabrication process, rather than the tension or compression between the films during the deposition process [18,19]. The workflow for our modeling process is shown in Fig. 1. The calculations presented here can be properly validated by the previous theoretical [16] and experimental results [11,12].

Computational methods
Our simulations have been performed based on the Newton's laws of motion and FE methods [20]. A set of differential equations have been established to describe both the motions of the sensing sheet and the substrate, in addition to the motion of the AFM tip. Here FE methods were used to solve the simultaneous differential equations numerically as the explicit analytical solution for such complicated partial differential equation can rarely be obtained. Then the resulting equations with the variables for the motions of the AFM tip, the sensing sheet, and the substrate, will be solved by using the numerical differential equation solver, 'NDSolve' functionality in Mathematica [20], which is fundamentally based on the FE methods. The commands option "automatic" has been used such that the most appropriate method was chosen to solve the equations. The methods used include forward Euler, midpoint The single-q model has been first benchmarked (not shown here). Then the appropriate potential energies for the sensing sheet and the substrate and the interaction between them have been designed, which will be discussed in detail in the next section. And another independent q-parameter for the substrate has been included. rule, linearly implicit Euler, and numerical approximation to locally exact symbolic solution methods [20].
As shown in Fig. 2, the graphene multilayer structure (in addition to the AFM tip) was chosen for simulations. To study the friction between the sensing sheet and the substrate, two independent parameters have been introduced as discussed in the Introduction section, q 1 and q s , to represent the internal motional dynamics of the sensing sheet and the substrate, respectively. The thicknesses of the film and the substrate will be described by the parameters defined later. The total potential energy U of the system can be expressed as in Eq.1. [16].
Here the first term on the right-hand side represents the flexibility of the AFM tip (the position of the tip is labeled by x) and cantilever, consisting of a spring with a spring constant of k between the tip (sliding over the periodic sheet) and the support, moving at a constant velocity of v (here t represents time). The energies for the sensing sheet and the substrate can be expressed, respectively, as below, Since the corrugation should be at a minimum of the potential energy for an undistorted sheet, this dependence is chosen to be a combination of a quadratic term and a quartic term. Here ν 2 and ν 2s have been set to be zero to separate the sensing sheet (q 1 ) and the substrate (q s ) [16]. As the tip must overcome the energy barrier to slide over the sheet, the corrugation would depend on q 1 and q s . The ν 4 and ν 4s characterize the rigidities (equivalently the thicknesses) of the sensing sheet and the substrate, respectively. Here it is assumed that ν 4 is slightly larger than ν 4s , which has been varied with a range of values to understand the behaviors of the relative motions and the force. The interactions between the tip and the sensing sheet/substrate are described in Eq.4 and Eq.5 [16] as below.    Table 1. ν 4s ranges between 0 M atom nm − 2 ps − 2 and 100 M atom nm − 2 ps − 2 .
Here a 1 and a s are the lattice constants for the sensing sheet and the substrate, respectively. V 1 and κ 1 (V s and κ s ) describe the behavior of the tip-sheet (tip-substrate) potential. The periodicities of the lattices on the sheet and the substrate are characterized by cosine functions in the second bracket of the Eq.4 and Eq.5, respectively. In Eq.5, f can be set to 0 (1) to suggest a fixed (movable) substrate. For all the studies presented here, it is assumed that there is a movable substrate (f =1). The most important part is the interaction energy between the sensing sheet and the substrate. Through the comparison with the previous work [16] on the tip-sheet/substrate interactions, the form of the interaction energy between the sensing sheet and the substrate is assumed as Here the lattice mismatch with a dimensionless parameter c has been taken into account, with the assumption that the relative motions  between films will influence the potential energy in the form of a parabolic function with a periodic feature characterized by a cosine function. According to this methodology, in principles the situation of thin films with more than one sensing sheets can be simulated by introducing additional degrees of freedom. Following Eqs.(1-6), the dynamical equation can be written as below.
Here m x , m q1 , and m qs are the masses for the atoms for the tip, the sensing sheet, and the substrate, respectively. η x , η q 1 , and η q s are the damping parameters. Through these simulations, the motional trends of all the layers can be analyzed when the stick-slip-friction phenomenon occurs. In principles, this model can be generalized to simulate any number of layers, which will lead to the corresponding simultaneous differential equations similar to the Eqs. (7)(8)(9). Below, the initial parameters used for our simulations have been summarized in Table 1, which will be varied later in the simulations. Here the unit and some of the values used in Ref. [16] for graphene have been adopted, which can be varied as needed when the experimental conditions change. It is also assumed that the substrate is crystalline. Notice that smaller values for the substrate than the sensing sheet were chosen to start with, i.e., κ s and V 1s are half of κ 1 and V 1 , respectively. a s is kept the same as a 1 , assuming the film and the substrate have the same lattice parameter. In addition, it is assumed that the atoms in the sensing sheet and the substrate have the same mass and the damping parameters. Here the chosen parameters can be varied while fixing the other parameters to the values as shown in Table 1. The boundary conditions include x(0) = 0;ẋ(0) = 0; q 1 (0) = 0;q 1 (0) = 0; q s (0) = 0;q s (0) = 0. Here M atom is the atomic mass of a hydrogen atom.

Results and discussions
The single-q model has been benchmarked (not shown here), which is consistent with the results presented in Ref. [16]. In Figs. 3-5, the dependence of the motion and forces on the ν 4s has been studied. As shown in Fig. 3, the behaviors of the individual motions of the sensing sheet and the substrate can be clearly observed. When ν 4s is small (this implies the potential of the substrate is small), the motions of the sensing sheet and the substrate are not coherent, as shown by the behaviors of the q 1 and q s in Figs. 3a and 3b. In other words, the strengthening process in terms of the friction force is not stabilized. By contrast, when the ν 4s is large, the motions between the sensing sheet and the substrate are more coherent (Figs. 3c and 3d), suggesting the adhesion effect is more robust. The slip-stick friction indicated by the motions of the films is illustrated by the blue arrows in Fig. 3a. In Fig. 4, the time evolution of the difference (the relative motion) between q 1 and q s , Δq = |q 1 − q s | is shown, which suggests the relative motions become more coherent as the substrate rigidity parameter ν 4s is increased. The slip-stick frictions can also be seen from Fig. 5 through the calculations of the forces, which shows the stability of the adhesion between the sensing sheet and the substrate, as illustrated in Fig. 3c and Fig. 3d.
In Figs. 6-8, the motions of the sensing sheet and the substrate are  Table 1. κ 1s ranges between − 50 M atom ps − 2 . and 50 M atom ps − 2 . shown as a function of time and κ 1s , which characterizes the sheetsubstrate interaction. In Fig. 6, for κ 1s = − 50M atom ps − 2 , − 10M atom ps − 2 , 10M atom ps − 2 , 50M atom ps − 2 , as the interaction between the sensing sheet and the substrate increases, the motions between the sensing sheet and the substrate become slightly more coherent, which suggests the improvement of the adhesion of the multi-layer films. This can also be seen in the corresponding calculations for the difference between the motions of the sensing sheet and the substrate in Fig. 7, as well as the time evolution of the force in Fig. 8. The calculation results for the motions have shown the behavior with a gradual stabilization along the time, except when ν 4s is small, suggesting the importance of the potential (or the stiffness) of the substrate [16]. The modeling presented here is partially in agreement with the previous experimental work, which shows the friction decreases when the substrates are thicker [21][22][23][24][25] based on the calculations of forces in Fig. 5 and Fig. 8, in addition to the stick-slip frictions [11,12]. However, the trend of the friction forces as a function of the number of atomic layers is not clear, which seems to be dependent on the substrate [24,25]. This can be studied in detail further by generalizing our modeling methodology to more than one atomic layer (i.e., more than one sensing sheet). The simulations here can also provide some information about the stiffness or the thickness of the sensing sheet and the substrate for the choice of the material as the substrate, to achieve mechanically robust wearable devices. The key is to figure out the relationship between the parameter ν 4s and the mechanical properties of the substrate, which can be explored in the future by experimental characterizations or materials modeling. In addition, the oscillations of the motions of the films and the forces between the AFM tip and the atom layer can be seen clearly, which suggests a slip-stick friction pattern, consistent with atomic-scale frictions, as shown in Figs. 3-8. In Figs. 9-11, the effect of V 1s (characterizing the interaction between the sensing sheet and the substrate) on the motions of the atomic layers and forces has been illustrated. When the magnitude of V 1s (with the negative sign, attraction forces) increases, the motions of the sensing layer and the substrate layer become more stabilized. When V 1s is between 0 M atom nm 2 ps − 2 and 40 M atom nm 2 ps − 2 , the difference between the motions of the sensing sheet and the substrate is still significant as shown in Fig. 10, which suggests V 1s is another important parameter similar to the ν 4s , and more significant than κ 1s . A negative V 1s can stabilize the motion as expected for an attraction energy. The stabilizing process between the sensing sheet and the substrate can also be seen from the force calculations in Fig. 11, in which the motions of the layers become stable when V 1s is 50 M atom nm 2 ps − 2 .
In Figs. 12-14, the behavior of the relative motion Δq as a function of the parameters V 1s , κ 1s , ν 4 and ν 4s has been studied and illustrated, at the fixed time t = 4000ps. As shown in Fig. 12(a), when ν 4s is fixed to a large value (corresponding to a more rigid substrate), V 1s will be a dominant parameter to determine the relative motion as compared with κ 1s . In addition, the motions of the two layers will be stable when V 1s has a relatively large negative value (strong attraction forces) or a relatively small positive value (weak repulsive forces). However, when ν 4s is small (Fig. 12b), the situation is more complicated, where V 1s need to be positive to stabilize the motions. Fig. 13 clearly shows that both V 1s and ν 4s need to be large to stabilize the motions between layers.
However, Fig. 14 suggests a more complicated situation when varying   Table 1. V 1s ranges between − 40 M atom nm 2 ps − 2 and 40 M atom nm 2 ps − 2 .   both ν 4 and ν 4s , the rigidities of the sensing sheet and the substrate. The main conclusion in Fig. 14 is that the motions can be stabilized when ν 4 is larger than ν 4s . Another stabilizing condition could be ν 4 is smaller than ν 4s when ν 4 is smaller than ~ 50 M atom nm − 2 ps − 2 .

Conclusions
In summary, a series of FE simulations for the two-layer graphene structure (one for the sensing sheet and the other for the substrate) has been performed according to the Newton's laws of motion to understand the friction between the layers. Another independent q-parameter for the substrate has been included to describe its motion explicitly, which can be used to distinguish the motions of the substrate from the sensing sheet. More importantly, by using this two-q model, the relative motions between the layers can be described properly, which is the key novelty of this work. The dependences of the motions and forces on V 1s , κ 1s , ν 4 and ν 4s have been studied. Depending on the parameters, the motions can be coherent, such as large values of V 1s , ν 4 and ν 4s , and both negative and positive V 1s can stabilize the motions. This suggests the importance of the rigidity of the layers and the interaction energy between the layers. However, under certain condition, such as the substrate with small υ 4s (mechanically fragile), the motions between them will be incoherent, thus causing problems for the proper adhesion between films. Moreover, the slip-stick frictions in the atomic layers can be found in the simulations, as seen in the AFM experiments in our two-q model. V 1s , ν 4 and ν 4s are important parameters for the interaction energy between sensing sheets and substrates and their potential energies to determine the stability of the motions, thus the adhesion performance of the multi-layer structures. The relationship between the film mechanical properties such as stiffness and the related parameter such as V 1s , ν 4 and ν 4s will therefore be very important for a better understanding of this type of dynamics, which will be investigated in detail in the future. In addition, the simulation methodology proposed here can be readily generalized further to simulate the friction among any number of atomic layers, thus providing a solid basis for the theoretical description of the motions in the real-world multi-layer devices.

Statement of originality
1) The paper has not been published previously, that it is not under consideration for publication elsewhere, and that if accepted it will not be published elsewhere in the same form, in English or in any other language, without the written consent of the publisher. 2) The paper does not contain material which has been published previously, by the current authors or by others, of which the source is not explicitly cited in the paper

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
All the code and data presented in the paper are available upon reasonable request.