Thermoelastic instability on a frictional surface and its implication for size effect in friction experiments

Recent laboratory friction experiments on large rock samples revealed that dynamic weakening, a remarkable reduction in the friction coefficient at elevated slip rates, occurs at lower slip rates in larger samples. There is a large difference between the sizes of natural faults and those in laboratory experiments. Therefore, it is crucial to understand the effect of size on rock friction. In the field of tribology, the interaction between frictional heating and thermoelastic effect has long been investigated. It was shown that higher slip rates than the critical value Vcr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V}_{\mathrm{cr}}$$\end{document} causes growth of temperature and normal stress heterogeneity (thermoelastic instability), and Vcr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V}_{\mathrm{cr}}$$\end{document} is proportional to the wavenumber of the heterogeneity. Severely heterogeneous normal stress may cause concentration of frictional power, thus locally activating dynamic weakening and leading to macroscopic weakening. Because a larger sample hosts a perturbation of a smaller wavenumber, it is expected to weaken at a lower slip rate than a smaller sample. In this study, a new numerical method was developed for analysis of thermoelastic instability based on the definition of memory variables and numerical approximation to the integration kernel, for the 2-dimensional problem of a planar fault embedded in an infinite medium. This method was advantageous over the standard integral equation method in terms of numerical costs. Numerical solutions with the new method on sinusoidal perturbations in the normal stress were compared with previously derived steady-state solution and its stability for validation. The typical thermoelastic properties of gabbro yield Vcr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${V}_{\mathrm{cr}}$$\end{document} in a range of experimentally adopted slip rates, indicating that the thermoelastic effect may play an important role in high-velocity friction experiments. Because the temperature rise and the resulting normal stress change smear out after the friction experiments, measurement of the temperature distribution in a sample during a friction experiment is important for further understanding the dynamic weakening and scale effect of rock friction.


Introduction
One of the most important discoveries in fault mechanics over the past quarter-century is the dynamic weakening of a fault at coseismic slip rates (Tsutsumi and Shimamoto 1997a). Frictional heating and associated temperature rise along a fault activate various weakening mechanisms (Di Toro et al. 2011 and references therein). Although the thermal properties of host rocks have been shown to affect dynamic weakening (Yao et al. 2016), frictional power density (i.e., the product of slip rate V and shear traction τ ) is the source of the temperature rise and a good proxy for the on-fault temperature at steady state (Di Toro et al. 2011).
Experiments which provided the basis for the notion of dynamic weakening were typically conducted on centimeter-scale samples. However, whether these can be extrapolated to larger scales, such as a natural fault or kilometer-or meter-scale subfaults, remains unclear. Otsuki and Matsukawa (2013) argued that the macroscopic friction coefficient of a spring-slider system decreases with the length of the slider owing to heterogeneous fault motion consisting of small nucleation and dynamic rupture propagation. A similar effect was reported in the dynamic earthquake sequence simulations for a fault embedded in an elastic space incorporating dynamic weakening (e.g., Noda and Lapusta 2010;Noda et al. 2011). In their simulation, the shear stress supported by the fault during an interseismic period was controlled by the coseismic friction strength.
Recent friction experiments on meter-scale rock samples by Yamashita et al. (2015) revealed a sample size effect at the onset of dynamic weakening. A small rock sample with outer 4.0cm and inner 1.7cm diameters shows dynamic weakening at a frictional power density of approximately 10 −1 MJm −2 s −1 , whereas a large sample of 1.5m length and 0.1m width weakens at a smaller frictional power density of approximately 10 −2 MJm −2 s −1 . Yamashita et al. (2015) explained this discrepancy by the concentration of frictional power due to heterogeneous normal stress σ n , modeled on the basis of the measurements of shear stress heterogeneity and fault surface observations. In their model, fault gouge production increased the local σ n , and dynamic weakening reduced the gouge production rate. These effects were all local, and the development of the heterogeneity was due to prescribed parameters, such as the proportion of the area of high σ n to the nominal fault area. Thus, careful preparation of the fault surface might rule out size effects. In addition, the absence of such an effect in centimeter-scale samples remains unexplained.
In the field of tribology, interaction of frictional heating and thermal expansion has long been studied. Barber (1967;1969) proposed positive feedback in heterogeneity in temperature T and normal stress σ n on the sliding surface through thermal expansion. Frictional heating increases T and causes thermal expansion, resulting in a dynamic concentration of σ n if conductive heat transfer is not efficient enough. This process is called thermoelastic instability. Dow and Burton (1972) looked for a separable form of perturbation with sinusoidal dependency along a frictional surface in a thermoelastic thin blade sliding perpendicular to it on a substrate. They showed that there was a critical slip rate V cr such that the perturbation grows at higher slip rates V > V cr . In the limit of an adiabatic substrate, V cr is proportional to the wavenumber of the perturbation. Burton (1980) showed a steadystate amplitude of the normal stress perturbation, which diverges at V = V cr . These previous findings may be relevant to fault mechanics and the dynamic weakening at high slip rates, but rarely discussed in interpretation of rock friction experiments. Spontaneous concentration of σ n and thus frictional power due to the thermoelastic instability may lead to local activation of dynamic weakening in laboratory rock friction experiments.
The aim of the present study was to investigate if the thermoelastic instability occurs under conditions of typical rock friction experiments, and if it explained the sample size effect. For this purpose, a simple 2-dimensional thermoelasticity problem is analyzed and evolution of a Fourier-mode perturbation is calculated. In the analysis, an integral equation was solved numerically with a new efficient method developed using memory variables (Noda 2022). This new method of thermoelasticity can be implemented in future dynamic earthquake sequence simulations (e.g., Lapusta et al. 2000). Finally, we discuss the implications of the analyses on sample-size effect of friction experiments and on natural fault behavior.

Model setting
Fault planes are rough, anisotropic, and typically show lineations which are referred to as slickenlines (e.g., Renard et al. 2006). The local height of the fault plane may correlate with the local σ n (e.g., Schmittbuhl et al. 2006), and thus the distribution of σ n may be idealized to be 1-dimensional. Estimation of the distribution of σ n on the experimentally slid surfaces using pressure-sensitive sheets ( Supplementary Fig. 5b in Yamashita et al. 2015) also supports the correlation between σ n and slickenlines. The geometry of the system is shown in Fig. 1. The fault lies normal to the y-axis at y = 0 , and the slip is in the z -direction. The compressional normal stress on the fault at the initial state of a uniform temperature is denoted as σ n0 (x) . Frictional heating was treated as a planar concentrated heat source at y = 0.
For the sliding surfaces used in laboratory experiments, a significant lineation may be absent in the initial state. However, linear 3-dimensional thermoelasticity under the assumption of small deformation or negligible Péclet number can be treated with a spectral method using Fourier transformation, and the idealization of the 2-dimensional problem (Fig. 1) corresponds to the investigation of Fourier modes whose wavenumber vector is in the xy-plane. The analyses presented in this paper are on the evolution and stability of these modes and are thus useful in discussing the experimental sliding surfaces even before the development of slickenlines. The distribution σ n0 may have a length scale in the z-direction, and if a much longer slip than the length scale is considered, convective heat transport may become significant. In addition, high-stress patches may migrate relative to the host rocks (Ettles 1986). These effects were not included in the present model, as they may be of minor significance after the development of slickenlines.
The system of governing equations for the 2-dimensional, plane-strain uncoupled thermoelasticity consists of the following equilibrium equations: the following stress-strain relations, and the Fourier's law, where σ is the component of the stress tensor positive in tension, u is the component of the displacement vector, t is the time, T is the temperature change, is the Lamé's constant, µ is the shear modulus, α is the change in pressure per change in temperature under fixed-displacement boundary condition, D is the diffusion coefficient of temperature, C is the specific heat capacity, and ω is the heat source of frictional heating. In plane-strain problems, α = β/K , where β and K are volumetric thermal expansion coefficient and bulk modulus, respectively. The side rocks between the fault are assumed to be identical, and thus the fault plane works as a rigid and adiabatic boundary. Therefore, the present problem is identical to the limit of an adiabatic substrate of the problem analyzed by Dow and Burton (1972) and Burton (1980). They considered a thin blade treated as a plane-stress problem, which is mathematically analogous to a plane-strain problem If the initial state is selected as a reference from which the displacement is measured, the compressional normal stress change on the fault is expressed as: and the frictional heating is expressed as: where q is the frictional power density which converts into heat, τ is the shear traction, V is the slip rate, and f is the friction coefficient. For typical rocks, f is in the range of Byerlee's law (Byerlee 1978) and does not vary much before the onset of dynamic weakening. V may change spatially and temporally, but the long-term average was controlled in laboratory experiments. The focus of the present study is on the evolution of σ n over the timescale of an entire experiment rather than individual stick-slip events so that f and V are assumed to be constant during fault slip.
Dow and Burton (1972) and Burton (1980) assumed a separable form of the temperature field having sinusoidal dependency on x and exponential dependency on y (7) �σ n (x, t) = −σ yy (x, 0, t), and t . In the present study, �σ n was expressed in terms of spatio-temporal convolution of the frictional power density and Green's function, and the full solution was calculated. This expression is useful in considering variable V and f and thus applicable to simulations of earthquake sequences. Xiao et al. (2021) derived Green's function for a concentrated line heat source: where δ is Dirac's delta function and H is the Heaviside step function. The relevant component of Green's function for the normal stress change on fault σ yy can be expressed as: where G is the response in the compressional normal stress on the fault for the concentrated heat source and γ is a function of Poisson's ratio ν, Using G , �σ n is expressed with a spatio-temporal convolution as follows: Fourier transformation with respect to x, leads to where k denotes the angular wavenumber along the fault. The Fourier transforms of G , G can be expressed analytically in a simple form as: Equations (8), (14), and (15), and the assumptions of constant f and V yield If it is assumed that the fault slip starts at t = 0 , the expression of the Fourier-transformed normal stress change becomes:

Steady-state solution
The integral equation Eq. (16) has a steady-state solution � σ nss (k) . Substitution of � σ n (k, t) = � σ nss and definite integral yields where v is a nondimensionalized slip rate The steady-state solution (Eq. (19)) was previously reported as Eq. (11) in Burton (1980), the critical slip rate (Eq. (20)) is the limit of an adiabatic substrate of Eq. (10) in Dow and Burton (1972). Equation (19) indicates that for fixed k there are two different velocity regimes: "low velocity" ( v < 1 ) and "high velocity" ( v > 1 ). � σ nss has the same sign as σ n0 in the low-velocity regime, and the opposite sign in the high-velocity regime.
The exponential growth rate c of the separable perturbation by Dow and Burton (1972) is given by ( iβ in Eq. (20) in Burton (1980)): Just note that there is a typological error in Eq. (20) in Burton (1980). c is positive if v > 1 , which is a sufficient condition for the thermoelastic instability: the steady state in the high-velocity regime is unstable. If v < 1 and if the separable perturbation is dominant, then the normal stress perturbation approaches to the steady state scaled by its initial value. Then the amplitude of heterogeneity in the normal stress can be reduced by carefully preparing the initial sliding surface to be very flat in the low-velocity regime. It is noteworthy that the full solution is different from the separable perturbation as is evident if one compared ∂T /∂t just after onset of slip. In a later section, c is used to verify the newly developed numerical method.

Numerical methodology
In general, f varies during fault motion, as well as V and σ n . If a characteristic value f * is selected, the corresponding V * cr can be defined from Eq. (20). Equations (8), (14), (15), and nondimensionalization using V * cr , σ n0 , and Dk 2 leads to where η = � σ n / σ n0 is the nondimensional Fourier-transformed normal stress change caused by the thermoelastic effect, s = Dk 2 t is the nondimensional time, and ω the nondimensional Fourier-transformed frictional power density: This could be, in general, calculated every timestep from spatial distribution of the variables with the Fast Fourier Transform technique. Here the argument k was omitted because a focus was put on a single Fourier-mode perturbation in the present study. The assumption of constant V and f and the natural selection of f * = f yield Note that the steady-state solution (Eq. (19)) is nondimensionalized to The calculation of the convolution by a standard numerical integration requires the storage of history, which increases with time. To avoid this, Eq. (22) was approximately transformed into a system of ordinary differential equations by approximation of the integration kernel and definition of memory variables, which was proposed by Noda (2022) in the implementation of poroelastic rebound into a dynamic earthquake sequence simulation. If an approximation is found, then Eq. (22) can be written as where φ i a memory variable, and its derivative is expressed as: In this study, this is integrated in time using a secondorder exponential time-differencing method, assuming a piecewise constant source term (e.g., Noda and Lapusta 2010). The detailed numerical algorithm is described in the Appendix.
A numerical approximation of the integration kernel (Eq. (27)) was obtained by solving the least-squares problem in the present study. Alternatively, it may be possible to construct an approximation by discretizing the inte- where p is the integration variable of the Laplace transformation, and s c = 1/p corresponds to the characteristic decay time. The discretization of the integral in Eq. (30) yields an approximation of erfc √ s . However, in doing so, a definition of the collocation points and a quadrature rule is required, and the best choice is not clear. Equation (30) was not used in this study. The least-squares fitting of Eq. (27) eventually optimizes the collocation points s c = s i and the coefficient a i for each decaying mode for the desired range of nondimensional time. Figure 2 shows the numerical approximation (Eq. (27)) obtained by a numerical solver of least-squares problems, scipy.optimize.leastsq, using up to 20 memory variables. 1000 data points were sampled from s = 10 −6 to 30 at grid points regularly spaced on a logarithmic scale. The solution search was performed with a constraint of n i=1 a i = 1 which guarantees equality in Eq. (27) at s = 0 . With n = 20 , the difference between the two sides of Eq. (27) is less than 10 −7 (Fig. 2a). More analyses of numerical costs and errors by comparison with a standard integral equation method are described in the Appendix. Figure 3 shows simulation results with v ≤ 1 . The lowvelocity regime v < 1 corresponds to low V , large k , and small samples because a small sample cannot accommodate heterogeneity of long wavelength. The change of normal stress heterogeneity due to the thermoelastic effect η increases while d η/ds decreases (Fig. 3a), and it asymptotically approaches to the steady-state value (Fig. 3b). The solution with v = 1 shows an almost linear increase with time. The characteristic nondimensional time of the decay to the steady state is apparently longer for larger v . At very low slip rate v ≪ 1 , η is small and almost follows a limit of η ≪ 1 . Equations (25) and (26) lead to This is shown with a thick black line in Fig. 3b. The simulation results demonstrate that the steady-state solution was stable in the low-velocity regime. In this regime, the heterogeneity in the normal stress during the friction experiments can be mitigated by reducing the initial heterogeneity σ n0 by, for example, carefully preparing the sliding surface.

Simulation results
In the high-velocity regime v > 1 , however, η increases unlimitedly, d η/ds increases (Fig. 4a), and thus, the steady-state solution η = v/(1 − v) < 0 is not realized during the friction experiment as shown by Dow and Burton (1972). A semilogarithmic plot of the full amplitude of the nondimensional normal stress heterogeneity σ n0 + � σ n / σ n0 = η + 1 (Fig. 4b) indicates that it increases almost exponentially with time. The growth of heterogeneity causes the concentration of frictional power and dynamic weakening or partial opening of the fault plane. These nonlinear effects are not included here and require further research to better compare the simulation with experimental data and to constrain the behavior of natural faults.
The growth rate relative to the steady-state solution (Eq. (26)) can be defined as: This may be different from the theoretical prediction for the separable perturbation by Dow and Burton (1972) (Eq. (21)). Figure 5 represents the nondimensionalized difference in the growth rates (c − c DB )/k 2 D . It decays towards zero, indicating that the separable perturbation becomes dominant after then initial run-in phase affected by the abrupt onset of sliding. Also, this asymptotic behavior support validity of the numerical solution.

Scale effect in laboratory experiments
It was shown that V cr (Eq. (20)) defines the low-and high-velocity regimes in terms of thermoelastic effect (Dow and Burton 1972;Burton 1980). The interaction of frictional heating and the thermoelastic effect leads to instability in the high-velocity regime, which causes unlimited growth of the amplitude of normal stress heterogeneity, as long as the assumptions adopted here are valid. The resulting concentration of frictional power conceivably activates dynamic weakening and decreases spatially averaged frictional resistance. V cr is proportional to the angular wavenumber k , and a large sample can host perturbation of a small k . Therefore, it is expected that a larger sample has a smaller V cr and thus shows thermoelastic instability and associated weakening at a smaller V . If the width of the sample is W , the normal stress perturbation of the system has a minimum angular wavenumber of Then the minimum V cr of the sample is Based on the simulation results (Figs. 3 and 4), the characteristic time of the normal stress evolution t c may be defined as s = 1 when t = t c , t c = c −1 DB is another candidate, but it diverges at v = 1 and thus has difficulty in comparison with experimental conditions. The corresponding characteristic time for the perturbation of minimum angular wavenumber is then (36) t max c = W 2 4π 2 D .  (32)) from the theoretical prediction for the separable perturbation (Eq. (20')) found by Dow and Burton (1972) The thermoelastic properties of gabbro, which are often used in high-velocity friction experiments, are summarized in Table 1. The P-wave speed V P , S-wave speed V S , density ρ , thermal conductivity κ , and heat capacity per mass c p were obtained from Schön (2015) and the volumetric thermal expansion coefficient β was obtained from Robertson (1988). Using these parameters, V min cr and t max c are expressed as follows: where the Poisson's ratio ν is given by: The parameters listed in Table 1 and the friction coefficient in the range of Byerlee's law (Byerlee 1978) f = 0.7 , yield These formulae give, for example, V min cr = 1.32cm/s and t max c = 2.91s for a sample of W = 1cm , and V min cr = 1.32mm/s and t max c = 291s for a sample of W = 0.1m . These values drop within the frequently adopted range of laboratory friction experiments for subseismic to seismic slip rates, implying the importance of the thermoelastic effect in those experiments.
The size effect discovered by Yamashita et al. (2015) decreases the macroscopic frictional power at the dynamic weakening by about an order of magnitude when W is increased from 1.15cm to 0.1m . At constant normal stress and V , the thermoelastic instability predicts an inverse proportionality between sample size and weakening speed. Because the data by Yamashita et al. (2015) contain experiments at different σ , it is replotted against V for a more direct comparison ( Fig. 6a and b). The data for the large sample ( W = 0.1m ) and small samples ( W = 1.15cm ) are plotted in Fig. 6a and b, respectively, together with the estimation of V min cr . For f in Eq. (37), the average friction coefficients for V < 2mm/s and V < 20mm/s were used for the large and small samples, respectively. In addition, V min cr was compared with that in the experiments by Tsutsumi and Shimamoto (1997a), as shown in Fig. 6c. They used hollow-cylindrical samples of inner and outer diameters of 25 and 16 mm, respectively, W = 4.5mm . f in Eq. (37) is estimated from average of friction coefficients in V < 0.1m/s . These comparisons indicate that V min cr reproduces the inverse Volumetric thermal expansivity β 1.6 × 10 -5 K −1 *2 Thermal conductivity κ 2.63 Wm −1 K −1 *1 Heat capacity per mass c p 1.01 kJkg −1 K −1 *1 Fig. 6 Comparison of experimental data and predicted critical slip rate for the thermoelastic instability V min cr . "MG", "MD", and "G" indicate metagabbro, monzodiorite, and gabbro, respectively. The stress values represent the normal stress of the experiments. a Date for a large sample ( W = 0.1m ) from Yamashita et al. (2015). b Date for small samples ( W = 1.15cm ) from Yamashita et al. (2015). c Date for even smaller samples ( W = 4.5mm ) from Tsutsumi and Shimamoto (1997a) proportionality between the sample size and the slip rate of the onset of dynamic weakening but underestimates it by about an order of magnitude. This may be due to nonlinear effects, such as partial opening that keeps local frictional power non-negative, or too much idealization to the infinite medium and neglect of efficient thermal conduction in sample holders made of stainless steel. The importance of the finiteness of the system was pointed out and has been studied in tribology (e.g., Lee and Barber 1993;Du et al. 1997). It is crucial to conduct realistic thermoelastic analysis of the experimental system for rock friction for more detailed comparison. For example, the experiments by Yamashita et al. (2015) show a series of stick-slip events in which V changes dynamically for orders of magnitude. The earthquake sequence simulation with the thermoelastic effect is probably a powerful tool in modeling of friction experiments and in interpretation of experimental results.
The present analysis of the evolution of Fourier-mode perturbation in the normal stress is applicable before various weakening mechanisms activate and decrease f locally where the normal stress σ and frictional power τ V concentrate. The condition for the onset of dynamic weakening was not directly inferred in the present analysis, because not only the heterogeneity, but also the uniform overall temperature rise on the sliding surface must be considered. If a sample of a friction experiment were infinitely large and we could conduct an infinitely long experiment, then the temperature would unboundedly increase unless τ V rapidly decreased to zero. Typical studies on friction experiments report steady-state friction coefficients that depend on the thermal properties of the sample assembly, as demonstrated by Yao et al. (2006). A more detailed comparison accounting for thermal properties of the sample assembly (e.g., Tsutsumi and Shimamoto 1997b) and the experimental apparatus deserves future study.

Implication for natural fault behavior
Natural faults that host repeating earthquakes experience a vastly wide range of slip rates, from orders of magnitude below the plate convergence rate ( ∼ 1nm/s ) to the coseismic slip rate ( ∼ 1m/s ). Equation (20) indicates that the critical wavelength of the thermoelastic instability cr also varies inversely proportional to the slip rate, V : The corresponding characteristic time is With the thermoelastic properties in Table 1 and f = 0.7 , this formula yields cr = 1.32 × 10 5 m and t c = 4.59 × 10 15 s ∼ 145Ma for V = 1nm/s , and cr = 1.32 × 10 −4 m and t c = 4.59 × 10 −3 s for V = 1m/s . Note that thermoelastic instability takes place for normal stress perturbation of a longer wavelength than cr . These numbers indicate that thermoelastic instability is not important when considering long-term processes, such as steady-state plate convergence, if the fault always creeps. It may become important during a seismic event or in acceleration processes towards the seismic slip rate. Implementation of the thermoelastic effect to the analysis of fault motion may not be straightforward. The critical wavelength cr must be resolved in a numerical simulation based on continuum mechanics, unless parametrization or coarse-grained fault constitutive laws are developed. In the simulation of a seismic event, a submillimeter grid interval may be needed, which requires huge computational resources. In addition, the existence of pore water in the natural fault zone may cause significant differences from the laboratory room-dried friction experiments. First, a change in σ n causes a change in the pore pressure, scaled by Skempton's coefficient (e.g., Skempton 1954), which diffuses out with time. This poroelastic effect effectively decreases f in the undrained condition, and the effect of the subsequent fluid flow must be is considered. Second, frictional heating causes not only thermal expansion of rocks but also pore-pressure build-up, called thermal pressurization (e.g., Sibson 1973). The permeability of fault rocks varies by orders of magnitude depending on the type of host rock and degree of brittle deformation. Therefore, the dominance of the thermoelastic effect studied in this paper over the poroelastic effect is still unclear and requires further study.

Conclusion
Thermoelastic instability has long been investigated in the field of tribology, but rarely discussed in experimental studies in rock friction. Dow and Burton (1972) and Burton (1980) showed by looking for a separable form of perturbation that there is a critical slip rate V cr for unlimited growth of heterogeneity in the normal stress and that V cr is proportional to the wavenumber of the heterogeneity. In order to investigate possibility of the thermoelastic instability during laboratory rock friction experiments and if it explained the sample size effect, a full solution for the interaction between the thermoelastic effect and frictional heating on a planar sliding surface in an infinite medium was numerically calculated in a framework of (43) t c = ρc p cr 2-dimensional uncoupled thermoelasticity. The Fourier transformation along the sliding surface leads to a simple integral equation of the amplitude of the Fourier-transformed normal stress σ for each angular wavenumber k . The integral equation can be solved efficiently by approximating the integration kernel to the summation of exponential decays and by defining memory variables. The new algorithm is suitable for future implementation into earthquake sequence simulations. The numerical solution reveals that the steady-state solution is stable at low slip rates v < 1 , and σ grows in an unlimited manner at high slip rates v > 1 , consistently with the previous studies. A severely heterogeneous normal stress should cause a concentration of frictional power, local activation of dynamic weakening mechanisms, and thus macroscopic weakening of the sliding surface. A larger sample hosts heterogeneity with a smaller wavenumber, and thus, a smaller V cr . Therefore, it is expected that a larger sample will show dynamic weakening at a smaller slip rate, as observed in the laboratory friction experiments. Estimations of the critical slip rate and time scale of evolution of the heterogeneity based on typical thermoelastic properties of gabbro drops within a range of typical experimental conditions and reproduce the tendency of inverse proportionality between the sample size and the slip rate at the onset of dynamic weakening with underestimation of about an order of magnitude. This underestimation is possibly because of partial opening of the sliding surface or too much idealization to an infinite medium. It has been shown that thermoelastic effect probably plays an important role during laboratory friction experiments. Because the thermoelastic effect diffuses out after the friction experiment, measurement of the temperature distribution in samples during friction experiments is important for further understanding the dynamic weakening and scale effect of rock friction.

Appendix: Evaluation of numerical methodology
The numerical method adopted in the present study was evaluated. For reference, Eq. (25) is solved with a standard integral equation method (IE), and numerical solutions with the present method using memory variables (MV) are compared with those with IE in terms of the numerical error and computational time.
As the integration kernel in Eq. (25), erfc √ s has a square-root singularity at s = 0 , it should be regularized by integration by parts for numerical accuracy as: Discretization with a regular mesh of interval s , linear interpolation between the grid points, and approximation of the integral by the midpoint rule leads to where the numerical integration kernel K i can be expressed as: The numerical solution at s = m s is then This method is semi-analytic and second-order accurate.
On the other hand, numerical integration requires preparation of the integration kernel and the storage of history, which increases with the number of timesteps. For a long simulation such as an earthquake sequence simulation, the temporal convolution must be truncated with an affordable length of the time window.
In the present study, these issues are resolved by approximating the integration kernel by the summation of exponential decays (Eq. (27)) and the definition of the memory variables φ i (Eq. (28)). φ i is integrated with time using a second-order exponential time-differencing method, assuming a piecewise constant source term based on a predictor-corrector approach (e.g., Noda and Lapusta 2010). Suppose that the variables at s = m s are known; then, the integration of Eq. (29) with a constant source term yields: Superscript * indicates that this is a first-order estimation of the solution. η at s = (m + 1)�s can be estimated as: (51) η * ((m + 1)�s) = n i=1 a i φ * i ((m + 1)�s).
A second-order accurate estimation can be obtained by integrating Eq. (29) in the latter half of this time step by: and The variables with * * were adopted as the numerical solution.
Convergence analyses were performed to evaluate the accuracy of the proposed method and to select a numerical parameter s . The cases with v = 0.1 , 1 , and 10 are analyzed, and solutions at s = 1 are compared, which is in the middle of the decay to the steady state in stable cases. Values of s from 1 to 2 −22 were tested, and the solution with the integral equation method with s = 2 −23 was used as the reference solution from which the numerical error was defined. Figure 7a shows the relative numerical errors for s = 1 . The solid symbols represent the IE, showing secondorder convergence. The numerical error is larger for larger v . With v = 1 s = 2 −3 , the solution shows growing oscillation so that the numerical error is above the range of the vertical axis. The open symbols represent the solution with MV, which yields a numerical error comparable to IE for relatively large s . With decreasing s , the numerical error with MV decreases to s 2 as expected from the second-order accuracy of the integration scheme, but it stops decreasing at approximately 10 −8 . This is the numerical error caused by the approximation of the integration kernel (Eq. (27)). (52) The numerical error due to the approximation (Eq. (27)) in the long-term ( s ≫ 1 ) behavior can be evaluated for low-velocity cases by comparing the analytic and numerical steady-state solutions. Equation (28) at a steady state ( d/ds = 0 ) yields and Eq. (28) leads to The nondimensionalized steady-state solution (Eq. (26)) indicates that n i=1 2a i s i ≈ 1 , and the present approximation (Fig. 2) yields: The computational time is always shorter for MV than for IE as long as studied (Fig. 7b). MV yields a linear scaling between the computational time and the number of time steps m unless m is so small that the overhead cost becomes significant. On the other hand, IE shows a scaling of approximately m 2 for large m . Therefore, MV is more advantageous for a simulation with a larger number of time steps.
In the simulations presented in the body of this paper, I choose s = 10 −4 and m = 10 5 indicated by dashed lines in Fig. 7a and b. The expected numerical error is approximately 10 −8 for cases with low slip rates ( v < 1 ), and of the order of 10 −5 for the case with the highest slip rate ( v = 10 ). These error levels as well as errors in the steadystate solutions (Eq. (56)), are too small to affect the discussion in this study. MV saves computational time by approximately one order of magnitude relative to IE.
The accuracy of MV is comparable to that of the standard semi-analytic method (IE) if the numerical (54) φ iss = 2vs i 1 + η ss , n i=1 2a i s i = 1 + 6.42 × 10 −7 . Fig. 7 Evaluation of the new numerical method developed in the present study. "IE" and "MV" indicate the integral equation method and the method using memory variables, respectively. a Relative error at s = 1 for v = 0.1 (blue), 1 (black), and 10 (red) for different timestep s . b Computational time for different number of timesteps m approximation to the integration kernel is sufficient and requires a smaller amount of numerical costs than IE in terms of not only the computational time but also the memory requirement. In addition, the constraint of constant s can be easily removed for MV, whereas the preparation of the numerical integration kernel K i requires constant s and the introduction of an adaptive timestep is not straightforward for IE. The focus of the present study was on the evolution of a single Fourier-mode perturbation; however, the new method developed here is presumably useful for the full simulation of fault motion, including a sequence of stick-slips and earthquakes (e.g., Lapusta et al. 2000).

IE Integral equation method MV
Method using memory variables

Acknowledgements
The author thanks T. Shimamoto, F. Yamashita, and T. Yamaguchi for their discussions, and the editor and anonymous reviewers. Comments by reviewers were useful in revising the manuscript. Especially, the anonymous reviewer #1 gave a vital comment informing the author of a classical paper on thermoelastic instability. Reference to it and related papers improved the manuscript significantly.
Author contributions HN contributed to the conceptualization, development of methodology, coding of simulation, design and execution of the parameter study, analyses of simulation results, and composition of the manuscript. The author read and approved the final manuscript.

Funding
The present study was supported by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan under its Second Earthquake and Volcano Hazards Observation and Research Program (Earthquake and Volcano Hazard Reduction Research) and by JSPS KAKENHI Grant Number 21H05201.

Availability of data and materials
No observational data were used for this study. The experimental data in Fig. 6a and b are from the supplementary information for Yamashita et al. (2015), and the data in Fig. 6c are read from Fig. 3 in Tsutsumi and Shimamoto (1997a, b). The simulation code used in this study is available upon reasonable request to the author.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
There is no competing interest.