Transition between the stick and slip states in a simplified model of magnetic friction

We introduce a simplified model of magnetic friction, and investigate its behavior using both numerical and analytical methods. When resistance coefficient $\gamma$ is large, the movement of the system obeys the thermally activated process. In contrast, when $\gamma$ is sufficiently small, the slip and stick states behave as separate metastable states, and the lattice velocity depends on the probability that the slip state appears. We evaluate the velocities in both cases using several approximations and compare the results with those of numerical simulations.


I. INTRODUCTION
The microscopic friction mechanism is an important subject of condensed matter physics and engineering [1][2][3][4], and various factors in this phenomenon, such as the lattice vibration and motion of electrons, have been studied [5][6][7][8][9].In particular, magnetic friction, the frictional force generated from the magnetic interaction between spin variables, has been extensively studied in recent years [10][11][12][13].To understand the mechanism of magnetic friction, many types of theoretical models have been proposed and investigated [14][15][16][17][18][19][20][21][22][23][24][25][26].These models differ from each other in several points such as the definition of dynamics, types of spin variables, and shape of the contact area.Accordingly, important features such as the relation between the frictional force and velocity vary with the choice of model.
In contrast, there is a well-known empirical law called the Dieterich-Ruina law for normal solid surfaces [27][28][29][30].This law generally has a complicated form that depends on the hysteresis.However, in the steady state, it can be expressed as the following simple relation: where A and B are constants.Note that most of the studies on the magnetic friction ignored the elastic deformation to consider only the influence of the magnetic interaction.Hence, it is difficult to investigate the change of the true contact area generated by the normal force.This is why the dependence on the normal force of the frictional force is not considered in Eq. (1).In our previous studies, we proposed models of magnetic friction that seemed to obey Eq. ( 1), at least in the steady state [31,32].In these models, magnetic structures behaved as a kind of potential barrier that prevented the lattice motion.Hence, we considered that the thermally activated process let the system obey the Dieterich-Ruina law, as in the case of normal solid surfaces [4].
However, even if the potential barrier exists, the frictional force, F , does not always obey the Dieterich-Ruina law.A well-known example is the Prandtl-Tomlinson model, which is one of the earliest theoretical models of friction.This model is composed of one particle moving on a sinusoidal potential pulled by a spring, which represents the contact point between the solid surfaces [33][34][35][36].If the temperature, T , is sufficiently low, this particle is trapped by the potential and cannot move smoothly.This point resembles the phenomenological explanation of the Dieterich-Ruina law.However, as previous studies have explained, frictional force F and velocity v obey the following rule: when F is smaller than a critical value, F c , and the system is trapped by the potential.Here, v c is the value of v when F = F c .This means that the prevention of the lattice motion by the potential barrier does not always result in Eq. (1).Note that the systems similar to the Prandtl-Tomlinson model always obey Eq. ( 2) when the Taylor expansion of the potential barrier exists, even if the potential is not the sinusoidal one.An example trying to explain the discrepancy between Eqs. (1) and (2) was proposed by Persson et al. [36].They considered a lubricated system and assumed that the first-order transition of the lubricant played a significant role.Generally, it is known that the friction process consists of two states, the stick state where the lattice is trapped by the potential, and the slip state where it moves smoothly.Persson et al. pointed out that the internal structure of the lubricant differs depending on whether the system is in the stick or slip state, and inferred that the system chooses the state with the lowest energy at any given time.Under this assumption, the point where the state with the lowest energy interchanges becomes the singular point of the effective potential, as shown in Fig. 1.The existence of this point prevents the system from obeying Eq. ( 2).Other studies also related the stick-slip motion of a lubricated system to the first-order transition of the lubricant [37][38][39].Magnetic bodies differ from these systems because they do not contain lubricants.However, the phase transition of the magnetization may affect the lattice motion.To understand the relation between the lattice motion and effective potential produced by the internal degrees of freedom of a solid, including a magnetic structure, a comprehensive investigation using both analytic and numerical approaches is required.
In our previous study, we succeeded in describing the behavior of a kind of infinite-range model in the thermodynamic limit using a mean-field analysis [32].However, we found that the relaxation time of this model diverged with increasing system size, and the actual behavior of the finite-size system in the steady state was different from that of the thermodynamic limit.Hence, it was important to understand the finite-size effect of this infinite-range model.However, an analytical discussion of this effect was difficult because this model had a complicated form.
Here, we introduce a simplified model to investigate a finite-size system.Specifically, we consider an infiniterange Ising model whose coupling constant J is the function of variable x, which represents the shift of the lattice.In our previous models, two magnetic bodies with different order parameters existed [31,32].However, for simplicity, our present model contains only one magnetic body with one order parameter, m.Although this model is largely simplified, it retains the essence of the magnetic friction, where the magnetic structure prevents the lattice motion.As we will explain later, we calculated the histogram of the state of this system, and found that the distinctions between the stick and slip states could be classified into two types.
In the reminder of this study, we first introduce the model in Sec.II.We then calculate its behavior in Sec.III, and finally summarize the study in Sec.IV.Supplementary discussions are also included as appendices A, B, and C. stick slip effective potential FIG. 1.The effective potential proposed by Ref. [36].They proposed that the transition between the stick and slip states occurs when the state with the lowest energy interchanges.

II. MODEL
We here consider N Ising spin variables {σ i } interacting with each other by the following Hamiltonian: Here, x is a real value representing the shift of the lattice, m is the magnetization per spin: m ≡ i σ i /N , and J(x) is the periodic function of x.For simplicity, we impose a periodic boundary condition with period 2π on x.Note that we consider a model composed of one magnetic body with one order parameter, m, whereas our previous models had two magnetic bodies with different order parameters.If we try to consider the situation that two magnetic bodies interact with each other, we should construct a model containing at least their magnetizations m a and m b .The model containing only one order parameter like Eq. ( 3) corresponds to the case that the other parameter is fixed as m b = const.by the strong internal interaction, for example.Furthermore, the ferromagnetic order appears in the case of this model, whereas the antiferromagnetic one appeared in our previous studies.However, even in this simplified model, the magnetization, m, behaves as a potential barrier that prevents lattice motion.That is, when m gets larger, the Hamiltonian also has the larger value and trap the lattice motion more easily.Hence, we consider that the essence of the magnetic friction is not lost through this simplification.As a concrete form of J(x), two types, A and B, are considered.Type-A is a piecewise linear function given by the following: and type-B is a sinusoidal function: This model coincides with the infinite-range Ising model with the exception that coupling constant J depends on the shift in the system, x.In this study, we let J 0 = J 1 = 1.In both type-A and B models, the maximum of J(x) is J(0) = J 0 + J 1 = 2.If the lattice was fixed at x = 0, the critical temperature of the spin variables would be given as T c = 2J(0) = 4.In the following, we consider the cases that T < 4 to investigate the relation between the magnetization and lattice motion.We introduce the time development of spin variables by the Glauber dynamics, and define one Monte Carlo step (MCS) as the unit time.Namely, the acceptance ratio, w, of each updating of the spin variable is defined by the following equation: where δE is the change in energy.In addition, the time development of x under this Hamiltonian and external force F ex are introduced.Specifically, we let x obey the following overdamped Langevin equation: where R(t) is the white Gaussian noise, R(t)R(t ′ ) = δ(t−t ′ ).Here, we introduce coefficient N , considering the situation that N spins move simultaneously.Substituting Eq. ( 3) and letting f ex ≡ F ex /N , this equation can be transformed as follows: Considering Eq. ( 9), magnetization m prevents the lattice motion when J ′ (x) < 0. In the steady state, the frictional force balances the external force, f ex .Hence, to investigate whether the system obeys Eq. ( 1), (2), or other rules, the v-f ex relation should be considered.
In the thermodynamic limit, N → ∞, the randomforce term of Eq. ( 9) disappears, and the dynamics of x become deterministic: As we already mentioned, the Hamiltonian given by Eq.
(3) has almost the same form as the normal infiniterange Ising model.Therefore, the time development of m in N → ∞ is given by the following [40]: Note that m coincides with its ensemble average, m , in this limit, because the fluctuation of m is O(1/ √ N ).According to the above discussions, the behavior in the thermodynamic limit, N → ∞, can be described by two deterministic equations (10) and (11).

III. CALCULATIONS
In this section, we investigate the behavior of this system using both numerical simulations and analytic methods.In the simulation, x is updated by the stochastic Heun method with time interval δt = 1/(nN ), where n is an integer.This means that we update x every 1/n step of the Glauber dynamics.In other words, we repeat the updating of x by Eq. ( 9) n times after one updating of the spin variables.This process is introduced because the change in x during one step of the Glauber dynamics is large when γ is small.We let n = 20 for γ ≤ 0.01, and n = 5 for larger values of γ.
A. Behaviors in the thermodynamic limit First, we discuss this system in the thermodynamic limit, N → ∞.As we explained in Sec.II, the time development of the system is described by Eqs. ( 10) and (11).We calculate (the time average of) velocity v in the steady state by solving these equations using the fourthorder Runge-Kutta method with time interval δt = 10 −4 .Fig. 2 shows the relation between v and f ex for the type-A model at T = 3 and γ = 0.0025.Here, two types of initial conditions are considered, (x, m) = (0, 0.05) and (0, 1), and the velocity is averaged over 10 3 < t < 4×10 3 .As seen in Fig. 2, there is a range where the velocity in the steady state depends on the initial value.We also calculated the time development of the finitesize system using numerical simulation.In this simulation, we chose the type-A model, let T = 3, γ = 0.0025, and f ex = 0.1, and started from the perfectly ferromagnetic state with x = 0. To obtain the data with error bars, each quantity was averaged over 6400 independent trials.The result is shown in Fig. 3 (a).We also calculated the relaxation time τ in Fig. 3 (b).Here, τ is defined as the time when m becomes smaller than the threshold value, m = 0.5 (the gray dashed line of graph (a).)As seen in the fitting line of Fig. 3 (b), τ behaves as the exponential function of N , and diverges rapidly in the thermodynamic limit.Hence, it is natural that the system shows dependence on the initial condition in this limit.The divergence of the relaxation time was also observed in the infinite-range model of our previous study [32], and made the consideration of the finite-size system difficult.Hence, we should investigate our present model, which is much simpler, to understand the behavior of a finite-size system.Note that we also performed similar calculations for the type-B model, but omitted the graphs because no qualitative differences between Figs. 2 and 3 were found.Here, τ is defined as the time when m becomes smaller than the threshold value, m = 0.5 (the gray dashed line of graph (a).)Data are plotted as the purple square points, and the gray dash-dotted line expresses the fitting line, log τ = 0.0284N − 0.0141.

B. Behaviors of the finite-size system
We then investigated the behavior of the finite-size system in the steady state, mainly using numerical simulation.In the simulation, we started from the perfectly ferromagnetic state with x = 0 and a sufficiently weak external force, f ex , and then gradually increased f ex .At each value of f ex , the first 2 × 10 5 MCSs were used for the relaxation, and the next 8 × 10 5 MCSs were used for the measurement.We mainly investigated the cases where T = 1.5, N = 150 and T = 3, N = 300.The value of N at each T was chosen so that the relaxation time became smaller than the simulation time.To obtain the data with error bars, each quantity was averaged over 64 independent trials.
To investigate the relation between the magnetic structure and lattice motion, we first calculated the histograms of m and x.To make these histograms, we divided the interval of x into 100 parts, and did not color the points where the (relative) frequency was lower than a threshold value, 10 −7 .Fig. 4 shows several examples of these, in the case of the type-A model.As shown in Figs. 4 (a-3) and (b-3), when γ is large, the value of magnetization m changes continuously with x.This means that the lattice motion and weakening of m occur simultaneously.This process is thought to be caused by thermal activation if f ex is small.Conversely, in a case where γ is small, two peaks are observed in the histograms (see Figs. 4 (a-1), (a-2), (b-1), and (b-2).)One of these exists near the fixed point of Eqs. ( 10) and (11), and the other has the form of a band with relatively small m.These peaks are thought to indicate the stick and slip states, respectively.This fact means that the stick and slip states are clearly separated as different metastable states.Hence, the change between these states resembles the first-order phase transition, when γ is sufficiently small.We also calculated the histograms of the type-B model, but omitted the graph because the qualitative behavior was similar to Fig. 4. Namely, even in the type-B model, a first-order-like transition was observed in the small-γ system, and the thermal activation process was found in the large-γ system.The border between these two behaviors are discussed in the appendix A. Judging from Eq. ( 9), the velocity tends to have a large value when γ is small.Hence, it is thought that fast lattice motion prevents the slip state from changing into the stick state, and consequently these states are separated from each other.The point that the stick-slip motion is related to the first-order transition resembles the inference proposed by Ref. [36].However, the lattice motion of our model (at a small γ) occurs only in the slip state, while their theory considered that the lattice moved with the stick and slip states switching.Note that the stick and slip states of our model do not interchange smoothly because they are metastable states.

C. Evaluation of the effective potential under small γ
It is difficult to calculate the theoretical probabilities with which the stick and slip states appear because the lattice motion is a non-equilibrium phenomenon.However, when γ is extremely small, the change in x is much faster than that in m.In this case, it is expected that the effective free energy under a given m could be evaluated by taking a kind of time average over x.In this section, we evaluate the probability of the slip state using this effective free energy.
Under the assumption of this section, the time scales of m and x are completely separate.Hence, we should take time average of Eq. ( 11) to discuss the contribution of J(x) to the time development of m: Note that Eq. ( 11) itself describes the behavior of the thermodynamic limit, as explained in Sec.II.However, we use this relation for the following discussions assuming that the effect of the fluctuation of m on the effective free energy is small.In Eq. ( 12), I 1 is defined as follows: where p m (x) is the probability distribution of the Brownian motion given by Eq. ( 9) under a fixed m in the steady state.
We then evaluate the effective free energy per spin, a eff (m), using Eq. ( 12).To extract the contribution of the entropy, the effective "internal energy," u eff (m), is defined by the following equation: where s(m) is the entropy of the system per spin, under a given m: To evaluate u eff (m), we consider the time development of m in a case where the energy of the system per spin is really given as u eff (m), and compare it with Eq. ( 12).The probability that one up (down) spin is reversed at each step of the updating, P u→d (P d→u ), is given by the following: Here, the change in energy is evaluated based on the fact that reversing the up (down) spin indicates a decrease (increase) in m by 2/N .Thus, the change in the expectation value, m , is expressed as follows: Using the assumption that N is large and m can be approximated as m ≃ m , this equation can be transformed as follows: Comparing (19) with (12), we can evaluate u eff (m) as follows: To calculate I 1 in Eq. ( 20), we should first evaluate p m (x).In the steady state, the Fokker-Planck equation corresponding to Eq. ( 9) with a fixed m is given as follows: The solution of this equation is given as follows: where A, B, and x 0 are constants, and u loc is defined as follows: Using the value of p m (x 0 ), Eq. ( 22) can also be expressed as follows: In particular, considering that e N βu loc (x0,m) converges to zero as x 0 becomes sufficiently large, Eq. ( 24) can be transformed as follows: In the case of the type-A model, the integral of Eq. ( 22) can be calculated because u loc (ξ, m) is a piecewise linear function of ξ.We will discuss the concrete form of p m for this case in appendix B. Conversely, in the general cases, including the type-B model, the exact calculation of this integral is difficult.Hence, we evaluate the approximate value.
First, we consider a case where u loc has a local minimum, x m .In this case, the system is thought to be in the local equilibrium state near this point.Hence, this case is thought to correspond to the stick state.When u loc has a continuous derivative, x m is determined by the following equation: In the type-A model, on the other hand, the local minimum of u loc appears at the discontinuous point of J ′ (x), i.e., x = 0. Hence, x m is given by the following: and exists if and only if If N is sufficiently large, x is located on x m with high probability: Indeed, if x m exists, the integral of Eq. ( 25) has a nearly constant value in the neighborhood of x m : ∞ x e N βu loc (ξ,m) dξ ∼ e N βu loc (Xm,m) .
Here, X m is the smallest argument of the local maximum of u loc , which is larger than x m .Hence, p m is proportional to e −N βu loc (x,m) when x ≃ x m , and this function draws a sharp peak at x = x m .Substituting Eq. ( 29) into Eqs.( 13) and ( 20), we obtain the following: The second line of Eq. ( 31) can be derived by the following equation: Eq.( 32) is certified by Eq. ( 26) in the case of the type-B model, and by the fact that x m is constant in the case of the type-A model.
In the slip state where |m| < m th , the integrand of Eq. ( 25) decreases monotonically and rapidly.Hence, the integral can be evaluated as follows: Substituting Eq. ( 33) into Eq.( 25), we obtain the following: Here, constant B ′ ≡ B/(N β) is determined by the normalization condition.This equation indicates that the duration of stay at each x is inversely proportional to velocity v ∼ (f ex + J ′ (x)m 2 )/γ.In short, in the case of the type-B model, we use Eq.(31) in the stick state, and evaluate I 1 using Eq.(34) in the slip state.This evaluation can also be applied to the type-A model.However, as already explained, we do not adopt it for the type-A model because the integral of Eq. ( 22) is calculated in appendix B.
Using a eff (m), the probability density of m, P (m), is expressed as follows: We compare the result of Eq. ( 35) and the numerical simulation to investigate the validity of the above discussions.The results of the type-A model are shown in Fig. 5.Note that in the histograms of this figure, we did not distinguish the value of x.Hence, these histograms are the integration of those calculated in Fig. 4 over x, except for the fact that the parameters have different values.As seen in Fig. 5, the results of the simulation seem to approach the estimation of Eq. ( 35) with decreasing γ.  14), (15), and ( 20), this quantity is calculated as follows: According to Eq. ( 13), I 1 (0) is given as I 1 (0) = 0. Calculation of I ′ 1 (0) is explained in the appendix C. Using the result of this appendix, Eq. ( 36) can be transformed as follows: Hence, when T < 2J 0 (= 2), a ′′ eff (0) is negative and P (0) becomes the local minimum.In contrast, it becomes the local maximum when T > 2J 0 (= 2).We let m 0 and m 1 be the arguments of the local minima of a eff (m) in the areas where m > m th and m < m th , respectively, as shown in Fig. 6.Then, defining q 0 and FIG. 6.
Definitions of m0, m1, and m2.As explained in Sec.III C, m0 and m1 are defined as the arguments of the local minima of a eff (m) in the areas where m > m th and m < m th , respectively.In Sec.III D, we also introduce the value m2, which is the argument of the local maximum in the area where m th < m < m0.q 1 as follows: the probabilities with which the stick and slip states appear are nearly proportional to q 0 and q 1 , respectively.Hence, when both of these local minima exist, the slip state appears with probability Note that P slip = 0 when m 1 does not exist, and P slip = 1 when m 0 does not exist.The average velocity of the slip state is given by the following: where Using Eq. ( 41), we obtain the expectation value of the velocity of this system: Note that although the velocity of the slip state given by Eq. ( 41) itself has a large value under a small γ, the ensemble average of the velocity given by Eq. ( 43) is small when the slip state rarely appears.We compare the results of Eq. ( 43) with those of the numerical simulation.
The results of the type-A model are shown in Fig. 7.As seen in this figure, the estimation of Eq. ( 43) seems to become accurate with decreasing γ, as is the case with the histogram of Fig. 5. Calculations similar to those shown in Figs. 5 and 7 in the case of the type-B model are shown in Fig. 8.As seen in this figure, the results of the simulations seem to converge to those of Eqs. ( 35) and (43) with decreasing γ, as is the case with the type-A model.As already mentioned, in the case of the type-B model, we used Eqs.(31) and (34), the rough approximation of Eq. ( 22), to evaluate a eff , whereas we used Eq. ( 22) itself for the type-A model.Hence the accuracy of the theoretical evaluation is thought to be higher for the type-A model than for the type-B model.However, a comparison of Figs. 5, 7, and 8 does not show this tendency.It is difficult to infer the reason why such behavior occurs, but one possible factor is the point that the dependence of histograms on γ is too large to observe such tendency.Another point to note is the form of p m for the type-A model under Eq.( 34).In the case of the type-B model, both p m at |m| = m th + 0 evaluated by Eq. ( 29) and that at |m| = m th − 0 evaluated by Eq. ( 34) have a sharp peak near x = π/2.On the other hand, in the type-A model, p m is calculated as a piecewise constant function at |m| < m th if Eq. ( 34) is adopted.This evaluation has an apparently different form from Eq. ( 29).It is thought that such a difference makes the evaluation of p m and a eff near |m| = m th using Eqs.(31) and (34) inaccurate and the improved approximation is required for the type-A model.

D. Comparison of cases with small and large values for γ
In a case where γ is large, as we saw in Fig. 4, the coexistence of the stick and slip states as metastable states breaks down.We evaluate the velocity of this case, assuming that the lattice motion is caused by the thermal activation process.Namely, we first define ∆a eff as the difference between the local maximum and minimum in the region where m th < |m| ≤ m 0 : and consider that the system is trapped in the stick state by the free energy barrier given as ∆a eff .Note that m 0 was already defined in Sec.III C, and m 2 is the argument of the local maximum in the area where m th < m ≤ m 0 , as shown in Fig. 6.Under this assumption, velocity v obeys the following relation: If external force f ex is stronger than a certain value, f c , the magnetic structure cannot trap the lattice motion, and Eq. ( 45) breaks down.The point f ex = f c is expressed as the critical value where the local maximum and minimum of a eff disappear and ∆a eff cannot be defined.We calculate the relation between − log(v/v c ) and f c −f ex using Eq. ( 45), and compare the results with those of the simulation.Here, v c is the value of v at f ex = f c .As the value of f c , we use the result of the theoretical evaluation even in the case of the simulation.The results are shown in Fig. 9.As seen in this figure, the v-f ex relation shows apparently different behaviors depending on γ, and is close to the results of Eq. ( 45) drawn as the black dashed lines (or curves,) when γ is large.The gray dash-dotted lines of Fig. 9 indicate the proportional relation of Eq. ( 2).In the case of the type-B model, these two lines are parallel to each other, which means that the evaluation of v using Eq.(45) obeys Eq. ( 2).This behavior results from the fact that Eq. (45) premises the continuous change between the stick and slip states.The value of log v given by Eq. ( 45) is not a linear function of f ex , because m 0 − m 2 decreases with decreasing f c − f ex .Instead, log v actually obeys Eq. (2) if J(x) does not have singular points.This mechanism resembles to that of the Prandtl-Tomlinson model.(For a more detailed explanation, see Ref. [36], for example.)In short, the lattice motion of our model in the large-γ case is essentially the same as that of the Prandtl-Tomlinson model, and the instantaneous phase transition proposed by Ref. [36] does not exist.Note that in the case of the type-A model, which has a singular point of J(x) at x = 0, the result of Eq. ( 45) does not obey Eq. ( 2).Finally, we calculate the quantity to investigate the difference compared to the Dieterich-Ruina law.This quantity becomes constant if the system obeys Eq. ( 1).When v is small, it is difficult to calculate R using the numerical differentiation of the simulation results, because the error bars of log v are large.Hence, we instead use the results of approximate evaluations based on Eqs. ( 43) and (45).The results of these two equations are compared in Fig. 10.Note that the values of R calculated by these equations do not depend on γ.As seen in Fig. 10, the results of Eq. ( 43), which correspond to the small-γ case, show a plateau in the small-f ex range.In our previous studies, we considered that the system obeyed the Dieterich-Ruina law in this range.However, even in this plateau, the value of R is not always constant.Hence, strictly speaking, the Dieterich-Ruina law does not always hold even in this range.It is possible that the models of our previous studies also de- viated from the Dieterich-Ruina law, but we could not find the difference using only the v-f ex graphs, like those shown in Figs.7 and 8.(b).We briefly discuss this deviation using Eq.(43).When f ex is small, the stick state appears with high probability.It means that q 0 ≫ q 1 , and consequently Eq. ( 43) gives the following relation: Here, we ignored slowly varying factors.Using Eq. (47), R is calculated as follows: Hence, in the small-f ex range, R stays constant and Eq. ( 1) holds, if a eff (m 0 ) and a eff (m 1 ) is the linear function of f ex .In other words, deviation from Eq. ( 1) in this range expresses the nonlinearities of the effective free energies, a eff (m 0 ) and a eff (m 1 ), as the function of f ex .Note that in the large-γ case, the discrepancy between Eq. ( 1) and the actual velocity is caused by the similar mechanism to that of the Prandtl-Tomlinson model, as we already explained.

IV. SUMMARY
In this paper, we considered a simplified model of magnetic friction in which the lattice motion is prevented by the magnetization of the system.Numerical simulations showed that the mechanism of the lattice motion differed depending on the resistance coefficient, γ.When γ is small, the stick and slip states exist as metastable states separated from each other, and the lattice velocity is determined by the probability of the slip state.In this case, the change between these two states resembles the first-order phase transition.As explained in Sec.III B, this transition has a different form from the inference of Ref. [36], which did not consider metastability.On the other hand, when γ is large, such a coexistence of two states is not observed, and the lattice motion is caused by the thermal activation process.This difference affects the relation between the lattice velocity and frictional force.As explained in Sec.III D, although the behavior of the lattice velocity at a small γ value under a sufficiently weak external force seems to be close to the Dieterich-Ruina law, Eq. (1), it does not always exactly coincide with this law.
Magnetic structures such as the magnetization are a kind of internal structure of the lattice.Hence, it is expected that the results of this study could be applied not only to magnetic friction, but also to the friction of normal solid surfaces.Such a generalization of this study could be the subject of future work.In particular, whether the stick and slip states are separated as metastable states or interchange with each other smoothly should also be investigated in the case of the usual solid surfaces.To apply our study to the realistic friction, simplification of the model seems to become the problem.For example, the overdamped Langevin equation assumes that the inertial term is much smaller than the damping term.Hence, our investigation on this equation under small damping constant γ seems to be strange.However, it does not mean that our study failed to grasp the essential behavior of the friction.In our model, the velocity of the slip state is fast when γ is small.Considering that the system does not have sufficient time to change the magnetization under the fast lattice motion, high velocity of the slip state is thought to be the direct cause of the separation of the stick and slip states.Hence, it is possible that such separation can be observed even though γ itself does not have a small value.In addition, in the case of the system which does not have the parameter corresponding to γ, the velocity of the slip state is thought to be the alternative criterion which judges whether such separation exists.We should consider these points carefully in future work.
The system size N should be also remarked.In this study, we considered the case that N is several hundred.If N gets larger, the relaxation time of the system diverges rapidly and the stochastic transition between two metastable states cannot be observed.On the other hand, if N = 1, metastable states themselves do not appear because there are no internal structures such as m.Hence, the coexistence of the metastable stick and slip states discussed in our study can be observed only when N is not too large nor too small.This restriction seems to be related with the true contact area of the solid surfaces.In the case of the friction of the normal solid surfaces, it is known that the frictional force is generated by the asperities, the junctions of protrusions of the surfaces.Hence, the true contact area of the surfaces is much smaller than the apparent one, but it is never composed of only one atom or molecule [2].In short, both the system we investigated in this study and the true contact area of the solid surfaces are few-body systems, and our result is an example that friction can have the behavior peculiar to such systems.Considering this point, our result implies that the dependence on the size of the contact area has an important role on the friction.
From the standpoint of statistical mechanics, it is also interesting that the infinite-range Ising model, which shows the second-order phase transition in the equilibrium state, shows the first-order-like transition when it is combined with the lattice motion.In a case where the model has the first-order phase transition or does not have any phase transition in the equilibrium state, whether the system behaves the same as our present model remains unclear.This point should also be studied in future.
In the case of the type-A model, using constants a 1 , a 2 , b 1 , and b 2 , Eq. ( 22) is expressed as follows: Solving these equations, we obtain relations between constants a 1 , a 2 , and a.When N is large, e −N βα+,mπ is much smaller than e −N βα−,mπ .Hence, we ignore the first term on the right-hand side of Eq. (B4b).Under this approximation, Eq. (B1) can be rewritten as follows:

FIG. 2 .
FIG. 2.Relation between (the time average of) velocity v and external force fex of the type-A model at T = 3 and γ = 0.0025.The red solid and blue dashed curves indicate the cases where the initial condition is given as (x, m) = (0, 0.05), and (0, 1), respectively.
FIG. 3. (a)Relaxation of magnetization m of the type-A model at T = 3, γ = 0.0025, and fex = 0.1.The red open square, green circular, blue triangular, and black closed square points represent the data with N = 100,200,300, and 400, respectively.(b)Relaxation time τ under the same condition as graph (a).Here, τ is defined as the time when m becomes smaller than the threshold value, m = 0.5 (the gray dashed line of graph (a).)Data are plotted as the purple square points, and the gray dash-dotted line expresses the fitting line, log τ = 0.0284N − 0.0141.

Comparing Figs. 5 (
a) and (b), the point m = 0 changes from the argument of the local minimum of P (m) to that of the local maximum with increasing T .It results from the change of the sign of a ′′ eff (m) at m = 0. Using Eqs. (

FIG. 7 .
FIG. 7. v-fex relation of the type-A model with differentγ values at (a)T = 1.5 and N = 150, and (b)T = 3 and N = 300.The red square, green circular, and blue triangular points represent the results of simulations at γ =0.0025, 0.01, 0.04, respectively.The red solid, green dashed, and blue dashdotted curves show the results of the theoretical evaluation using Eq.(43) for each value of γ.

FIG. 8 .
FIG. 8. v-fex relation of the type-B model with different γ values at T = 1.5 and N = 150.The meanings of the points and curves are the same as those in Figs. 5 and 7.

FIG. 9 .FIG. 10 .
FIG. 9. Relation between − log(v/vc) and fc − fex in cases for (a-1)type-A model at T = 1.5 and N = 150, (a-2)type-A model at T = 3 and N = 300, (b-1)type-B model at T = 1.5 and N = 150, and (b-2)type-B model at T = 3 and N = 300.The red open square, green circular, blue triangular, and purple closed square points indicate the results of simulations at γ = 0.0025, 0.16, 2, 16, respectively, and the black dashed lines or curves were calculated using Eq.(45).We also plotted the proportional relation of Eq. (2), as shown by the gray dash-dotted lines.