Friction on layered media: How deep do phonons reach?

Sliding friction is a complex phenomenon, involving surface mechanics [1–8] and production of heat in the surrounding media. In the latter process, the relevant degrees of freedom have been found to include electronic [9–13] and phononic ones [8,11–23]. The role of phonon modes can be isolated by tuning them while keeping the surface properties and contact mechanics unchanged. This has been achieved by inducing a phase transition in the solid [13,19], or by changing an external electric field [20].

Studying friction on layers of different thickness is another way of achieving this goal; it is especially insightful as it not only allows material properties to be tailored while limiting changes in the contact surface, it also reveals how deep friction feels into the material. Many numerical and nanoscale experimental studies have been done on the effect of layers [24–36], with a variety of interesting behaviors. For example, while adding graphene layers between sliding bodies can substantially reduce friction [33], both increases [24–26] and decreases [27,30–32] in friction have been observed as the number of graphene layers is increased. Friction has also been observed to decrease with thickness for other layered materials, such as molybdenum disulfide and niobium diselenide, when they are weakly anchored to the substrate [32]. A numerical study, in contrast, found that friction increases with thickness, when the bottom layer absorbs incoming phonons [34]. Additionally, it has been reported that friction is smaller for strongly anchored samples compared to weakly anchored ones, suggesting that the changes in friction could be results of changes in local deformations [28,29,36]. The complexity of the observed behavior calls for a theoretical analysis that systematically addresses the dependence of phononic damping on layer thickness, boundary conditions, and material properties.

In this manuscript, we analyze friction of a nanoscopic object on a 3D planar coated substrate, treating phonon modes via a field theory of viscoelasticity that includes phonon attenuation. We find that friction arises due to traveling spherical waves, cylindrical surface waves, or evanescent waves, each dominating in different regimes of material properties and coating thickness. A finite viscosity thus not only results in phonon attenuation, but also, here more importantly, in losses from evanescent waves. Consequently, friction shows drastically different dependencies on coating thickness ranging from short range to long range behavior, and can increase or decrease with coating thickness. These regimes are determined by the phonon attenuation coefficients and the refractive index.

Consider a probe coupled to an isotropic solid filling the space $z\leq 0$ (see fig. 1). The probe oscillates parallel to the surface, so that its x coordinate is $X(t)=Re\{X_0 e^{i\omega t}\}$. The coupling with the surface causes a force acting on the probe, whose x-coordinate is $F(t)=Re\{F_0 e^{i\omega t}\}$. The damping coefficient or friction coefficient of the probe Γ is defined as $\Gamma=\omega^{-1}Im\{F_0/X_0\}$ [37]

Zoom In Zoom Out Reset image size

The specific form of coupling is not important for the conclusions of this manuscript, as detailed in the Supplementary Material Supplementarymaterial.pdf (SM). We thus use a simple linear one, which allows us to find analytic expressions,

$$\begin{equation} F(t)=n_\mathrm{A}\kappa\iint_{-\infty}^{\infty}\mathrm{d}^{2}{\boldsymbol{r_\parallel}} (-X(t)+u_x(\boldsymbol{r_\parallel},t))e^{-\frac{\boldsymbol{r_\parallel}^2}{l^2}}\!, \end{equation} \tag{ 1 }$$

where $\boldsymbol{r_\parallel}=(x,y,0)$ marks a position on the surface, and ${u}_x(\boldsymbol{r_\parallel},t)$ is the x-component of the phonon field at z = 0 (introduced below). $n_\mathrm{A}$ is the particle number per unit area [38], and κ is the coupling strength. Equation (1) contains a Gaussian envelope of width l, introducing a length scale of the interaction range or probe size.

The first term in eq. (1) is the force in absence of phonon excitations; it is in phase with X(t) and does not contribute to the damping coefficient Γ. The second term is the force due to excitations of the phonon field $\boldsymbol{u}$, treated via a theory of viscoelasticity, i.e., using a Kelvin-Voigt model [39–42],

$$\begin{equation} \left(c_\mathrm{L}^2-c_\mathrm{T}^2\right)\nabla \nabla\cdot\boldsymbol{u}(\boldsymbol{r},\omega)+ c_\mathrm{T}^2 \nabla^2\boldsymbol{u}(\boldsymbol{r},\omega)=-\omega^2\boldsymbol{u}(\boldsymbol{r},\omega). \end{equation} \tag{ 2 }$$

The coupling to the probe enters eq. (2) via a time-dependent boundary condition. The resulting solution for u_x contains a part that is phase shifted with respect to X(t), yielding the damping coefficient Γ [16,17]. The oscillating probe excites phonons, whereby its motion is damped. The solution of this problem proceeds via the Green's function of eq. (2) (SM).

Equation (2) is a widely applicable model of phonon dynamics [39–41], and related models have been influential in understanding friction [16,18]. It contains two fundamental modes with longitudinal, $c_\mathrm{L}(\omega)=\sqrt{\frac{3K+4\mu-i\omega(3\xi+4\eta)}{3\rho}}$, and transverse, $c_\mathrm{T}(\omega)=\sqrt{\frac{\mu-i\omega\eta}{\rho}}$, speeds of sound. K and μ (ξ and η) are the bulk and shear elastic (viscous) moduli, respectively, and ρ the mass density. $c_\mathrm{T}$ and $c_\mathrm{L}$ are complex due to finite viscous moduli, i.e., eq. (2) contains phonon attenuation. Microscopically, such attenuation may be caused by phonon-phonon, phonon-electron, or phonon-defect scattering. For simplicity, we assume ${c_\mathrm{L}(\omega)}/{c_\mathrm{T}(\omega)}=\sqrt{3}$, independent of ω, a good approximation for solids [16,43–45]. $c_\mathrm{L}$ thus drops out of the discussion.

The bulk situation, where the coating thickness d in fig. 1 is infinite, was extensively studied previously [16,17,42,46]. In this case, the probe excites spherical waves $\sim e^{i \omega r_\parallel/c_\mathrm{T0}}/r_\parallel$ (SM), and the damping coefficient reads

$$\begin{equation} \Gamma^{\infty}=\frac{\kappa^{2}n_\mathrm{A}^{2}l^{4}}{c_\mathrm{T0}^{\prime\,3}\rho_0}\left(\zeta^{\mathrm{ela}} + \frac{\zeta^{\mathrm{vis}}\eta_0}{c_\mathrm{T0}^\prime\,\rho_0l}+\mathcal{O}\left(\frac{1}{Q^{2}},\frac{l^2}{\lambda^2}\right)\right)\!, \end{equation} \tag{ 3 }$$

with $\zeta^{\mathrm{ela}}\approx1.29$, $\zeta^{\mathrm{vis}}\approx1.72$, and $c_\mathrm{T0}^{\prime}$ the real part of $c_\mathrm{T0}$. Equation (3) is valid for $l \ll \lambda$, with phonon wavelength $\lambda=2\pi\,c_\mathrm{T0}^{\prime}/\omega$. For typical speeds of sound and l in the nanometer range, the corrections to eq. (3) are small for $\omega \lesssim 10^{9}\ \text{Hz}$. This number, $10^{9}\ \text{Hz}$, is large compared to typical frequencies excited by the probe, since they are expected to be in the range $\omega\sim10^{2}\text{--}10^{5}\ \text{Hz}$ [11,14,15,19,20,30–32]. Additionally, eq. (3) assumes the quality factor $Q(\omega)=\frac{\mu_0}{\omega\eta_0(\omega)}$, the ratio between phonon decay length and phonon wave length, to be large compared to unity. This seems a justified assumption as well, as seen in estimates of the order of $Q\sim10^4$ [47,48] (see table 1 for experimental parameters). Despite Q being large, phonon attenuation is essential for friction since, as shown below, it dominates the behavior in certain regimes.

Table 1:. Parameters typical for AFM experiments on solids: radius of contact l, frequency ω, mass density ρ, and shear elastic modulus μ.

la	$\omega^{b}$	$\rho^{c}$	$\mu^{c}$
$1\ \text{nm}$	$10^{3}\ \mathrm{\frac{1}{s}}$	$10^{-23}\ \mathrm{\frac{kg}{{nm}^{3}}}$	$10\ \mathrm{\frac{kg}{{nm\,s}^{2}}}$
a References [19,27,32].
b References [19,20,30–32].
c Reference [53].

The two leading terms in eq. (3) have fundamentally different physical origins. The first, called the elastic contribution, corresponds to transport of energy by traveling waves, and it persists without phonon attenuation [16,42]. The second, called the viscous contribution, is proportional to viscosity $\eta_0$. We discuss its physical origin below. Their relative weight is the dimensionless viscosity,

$$\begin{equation} \tilde\eta=\frac{\eta_0}{c_\mathrm{T0}^{\prime}\rho_0 l}=\frac{\lambda}{2\pi Q l}, \end{equation} \tag{ 4 }$$

which depends on material properties and probe size l. The numbers of table 1 for typical AFM conditions imply $\tilde\eta\sim10^{5}$. We thus first discuss the case of $\tilde\eta\gg1$, where the bulk case of eq. (3) is dominated by the viscous contribution (using $c_\mathrm{T0}^{\prime} = \sqrt{\mu_0/\rho_0}+\mathcal{O}({Q^{-2}})$)

$$\begin{equation} \lim_{\tilde\eta\gg 1}\Gamma^{\infty}=\frac{\zeta^{\mathrm{vis}}\kappa^{2}n_\mathrm{A}^{2}l^{3}\eta_{0}}{\mu_{0}^{2}}. \end{equation} \tag{ 5 }$$

For finite values of d, we note that the boundary condition at the interface between coating (labeled with the subscript 0) and substrate (subscript 1) is determined by the refractive index $n(\omega)=c_\mathrm{T0}(\omega)/c_\mathrm{T1}(\omega)$ and $\rho_0/\rho_1$ [49]. It is insightful to start with the limit n = 0 (see footnote  ¹ ), which yields a Dirichlet boundary condition (DBC) of $\boldsymbol{u}(\boldsymbol{r_\parallel},-d,\omega)=\boldsymbol{0}$ [40,49]. In this case the substrate is much stiffer than the coating, and phonons are totally reflected at the interface, with phase shift π.

Figure 2(a) shows Γ as a function of d for n = 0 and $\tilde\eta\gg 1$, growing linearly in d, and saturating to the bulk value of eq. (5) for large d, with a cross-over length scale around $d\approx l$. Indeed, we find that, for small d, the leading order of Γ is linear in d,

$$\begin{equation} \Gamma=\frac{\pi\kappa^{2}n_\mathrm{A}^{2}l^{2}\eta_{0}\,d}{2\mu_{0}^{2}}+\mathcal{O}({d^3})\stackrel{\tilde\eta\gg 1}{=} \frac{\pi\Gamma^{\infty}}{2\zeta^{\mathrm{vis}}}\frac{d}{l}+\mathcal{O}({d^3}) \quad \text{(DBC)}. \end{equation} \tag{ 6 }$$

Γ vanishes as $d\to0$, because the coating, when placed on a stiff substrate, supports fewer phonons as $d\to0$. In the second step in eq. (6), we used eq. (5) to replace $\Gamma^\infty$, making the mentioned saturation to $\Gamma^\infty$ at $d\approx l$ apparent. The dependence of Γ on d is short range, set by probe size l of eq. (1).

Zoom In Zoom Out Reset image size

How is this short range decay possible over distances much smaller than the wavelength? The motion of the probe excites not only (attenuated) traveling waves, but also evanescent waves that decay within a range of l. In the presence of finite viscosity, they contribute to energy absorption, and thus to the damping coefficient Γ. In fig. 2, this mechanism outweighs the energy transported by traveling waves, which is why eqs. (5) and (6) carry $\eta_0$ as a factor.

The opposite limit of boundary conditions, n → ∞, corresponds to a freestanding coating or a much more compliant substrate, and the waves obey a Neumann BC (NBC), $\boldsymbol{\hat{e}_z}\cdot\boldsymbol{\sigma}(\boldsymbol{r_\parallel},-d,\omega)=\boldsymbol{0}$ with stress tensor $\boldsymbol{\sigma}$ and surface normal $\boldsymbol{\hat{e}_z}$ [49]. Phonons are totally reflected without phase shift. This renders Γ fundamentally different from the case of n = 0 as shown in fig. 2(b); Γ diverges for small d and converges to the bulk value on a scale of $d\approx Ql$. Expanding this case for small d yields

$$\begin{equation} \Gamma=\frac{11\pi^{2}\kappa^{2}n_\mathrm{A}^{2}l^{4}}{64\mu_{0}|{\omega}|d}+\mathcal{O}{(d)}\stackrel{\tilde\eta\gg1}{=}\frac{11\pi^{2}\Gamma^{\infty}}{64\zeta^{\mathrm{vis}}}\frac{Ql}{d}+\mathcal{O}({d}) \quad \text{(NBC)}. \end{equation} \tag{ 7 }$$

In this case, the probe excites traveling surface waves, $\sim\,e^{i\omega r_\parallel/c_\mathrm{T0}}/\sqrt{r_\parallel}$, i.e., the coating oscillates like a freestanding 2D sheet. Energy is transported along the surface, rather than absorbed; eq. (7) is independent of viscosity $\eta_0$. As the sheet gets more compliant with $d\to 0$, amplitudes of excitations get larger, and formally diverge as $d\to0$. The second equality of eq. (7) implies the mentioned saturation length of $d\approx Ql$. With l on the scale of nanometers, this length is of the order of microns. It will be interesting to compare these surface modes to so-called puckering [32,36] or ploughing [28,35] identified in previous work.

The cases of n = 0 and n → ∞ provide a reference for the discussion of an arbitrary refractive index n. For simplicity, we assume that the quality factors of the coating and substrate are identical, making $n(\omega)$ real, and use $\rho_0=\rho_1$ (see footnote  ² ). For finite values of n, we note a numerical challenge in evaluating the longitudinal modes, so that we restrict the shown data to transverse waves (SM).

Figure 3 shows Γ as a function of d for various n, at $\tilde\eta=10^5$ and $Q=10^4$; Γ is monotonic in d, and stays within the bounds of the limiting cases of n = 0 and n → ∞. Importantly, for a thin coating, $d\to0$, Γ approaches the bulk value of the substrate, indicated as bars on the y-axis. These are obtained by eq. (3) with the material parameters of the substrate ^{3 , 4} . Γ thus varies between the bulk results of the substrate $(d\to0)$ and of the coating $(d\to\infty)$.

Zoom In Zoom Out Reset image size

The curves of fig. 3 up to $n=2^3$ can be understood by these two limiting cases, and a transition on the length scale l. These cases are thus dominated by evanescent waves. The dimensionless viscosity of the substrate equals $\tilde\eta/n$, and goes down for large values of n. The damping coefficient of the bulk substrate, i.e., the behavior at small d, is thus dominated by traveling waves at large n. Indeed, for n → ∞, the behavior of fig. 2(b) is approached. For intermediate values of n ($n=10^2$ in the graph), a two-step decay occurs, with evanescent waves for $d\lesssim l$ and $d\gtrsim Ql$, and traveling surface waves for $l\lesssim d\lesssim Ql$.

What about the case of $\tilde\eta\ll1$ where the limit d → ∞ is dominated by traveling waves? Despite less experimental relevance to AFM experiments, we include this insightful case in fig. 4 for completeness, focusing on n = 0 (see SM for n → ∞). As seen in the figure, for $d\ll \lambda$ the curves are similar to fig. 2(a). This is because, for n = 0, the coating does not support traveling waves for d ≲ λ (SM). For d ≳ λ, traveling waves are excited, yielding a sharp transition between evanescent and traveling waves at $d\approx\lambda/4$, followed by peaks. These are due to interference effects between outgoing and reflected traveling waves (SM). For d ≫ λ, the bulk value of eq. (3) is approached.

Zoom In Zoom Out Reset image size

This discussion allows identification of mechanism regimes, depicted in fig. 5. The limit of d → ∞ is dominated by evanescent waves for $\tilde\eta\gg1$, and by traveling waves for $\tilde\eta\ll1$. This is indicated in fig. 5 by using different colors in the inset graphs. The line $n=\tilde\eta$ separates the same for $d\to0$, i.e., curves that begin with blue or red. This limit can be found from ⁵

$$\begin{equation} \lim_{d\to 0}\Gamma= \Gamma^\infty n^3 \frac{\zeta^{\mathrm{ela}}+\zeta^{\mathrm{vis}}\frac{\tilde \eta}{n}}{\zeta^{\mathrm{ela}}+\zeta^{\mathrm{vis}}\tilde \eta}. \end{equation} \tag{ 8 }$$

Another separator is n = 1, dividing curves where $\Gamma^\infty$ is larger than the limit of $d\to 0$ from the opposite. The last curve, $n=\sqrt[3]{\tilde\eta}$, together with $n=\tilde\eta$, bounds two regimes with multiple transitions between traveling and evanescent waves. The latter occur because of the different d-dependence of traveling and evanescent waves.

Zoom In Zoom Out Reset image size

This framework identifies the main phononic mechanisms for damping between a probe and a coated planar surface, with vastly different behaviors. These regimes are expected to occur in experiments and simulations, whenever the thickness of coatings can be changed without changing contact mechanics. Figure 5 illustrates that pronounced dependence on layer thickness is expected, e.g., if the two materials show a refractive index very different from unity. The discussed isolated frequencies apply in non-contact measurements, where narrow frequency bands are excited [11,13–15]. In sliding experiments, one expects a spectrum of frequencies, so that the results reported here need to be averaged accordingly. It is important to note that the contributions by evanescent waves are frequency independent within the range given below eq. (3), making this average trivial. Friction from traveling waves has frequency-dependent features, see eq. (7), manifested as the peaks in fig. 4.

The presented model is simple and provides analytical results, which we expect to improve understanding of friction phenomena. The quantitative translation to the case of sliding motion needs to be investigated in future work [50]. We note however qualitative agreement with studies on graphene. Experiments [32] and simulations [27] reported friction on freely standing graphene flakes, i.e., $n\gg1$, finding a decrease with increasing number of layers. The same trend is observed for mounted graphene layers (or flakes) [27,30–32]; graphene is much stiffer ($\mu\sim1\ \text{TPa}$ [45,51,52]) than metal substrates ($\mu\sim10\ \text{GPa}$ [53,54]), so that the case of $n\gg1$ applies too. Previous work attributed this behavior to the enhanced rigidity for thicker samples [27,32], leading to the suppression of surface waves, in agreement with our findings.

The framework of eq. (2) naturally includes Hertzian contact theory [40,55] in the limit of small frequencies. It is interesting to remark that the identified evanescent waves can equally be found from the slowly moving distortion field of contact theory. Moving this distortion field dissipates energy due to viscosity (SM). This way, this approach can be linked to a variety of other approaches based on contact theory [4,7,18,50] and local deformations [28,29,36]. Natural extensions include nonlinear systems, where, e.g., deeper indentations [28,29,36] can be studied. Lateral confinement was also found to be of importance [8], and it provides an interesting additional possibility to identify involved phonon modes. This can be address in the presented framework in future work.

This work was funded by the Deutsche Forschungsgemeimschaft (DFG, German Research Foundation), grant No. 217133147 via SFB 1073 (Project A01). We thank Richard L. C. Vink for stimulating discussions. In memory of Philip Rauch.

Data availability statement: All data that support the findings of this study are included within the article (and any Supplementary files).