Friction on layered media: How deep do phonons reach?

We theoretically study the frictional damping of a small probe object on a coated planar surface, analyzing the resulting phonon modes via a theory of viscoelasticity. Three different types of excitations are found to contribute to friction in distinct ways: traveling (3D) spherical waves, traveling (2D) surface waves, and evanescent waves. While traveling waves transport energy away from the probe, determined by long range elastic properties (wavelength), evanescent waves transform energy into heat in a near-field range, characterized by the size of the probe. Thus, fundamentally different behaviors are predicted, depending on coating thickness and material properties.

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][25][26][27][28][29][30][31][32][33][34][35][36], with a variety of interesting behaviors. For example, while adding graphene layers between sliding bodies can substantially reduce friction [33], both increases [24][25][26] and decreases [27,[30][31][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 ≤ 0 (see Fig. 1). The probe oscillates parallel to the surface, so that its x coordinate is X(t) = Re X 0 e iωt . The coupling with the surface causes a force acting on the probe, whose x-coordinate is F (t) = Re F 0 e iωt . The damping coefficient or friction coefficient of the probe Γ is defined as The specific form of coupling is not important for the conclusions of this manuscript, as detailed in the Supplemental Material (SM). We thus use a simple linear one, which allows us to find analytic expressions, where r ∥ = (x, y, 0) marks a position on the surface, and u x (r ∥ , t) is the x-component of the phonon field at z = 0 (introduced below). n 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 u, treated via a theory of viscoelasticity, i.e., using a Kelvin-Voigt model [39][40][41][42], 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][40][41], and related models have been influential in understanding friction [16,18]. It contains two fundamental modes with longitudinal, c L (ω) = 3K+4µ−iω(3ξ+4η) 3ρ , and transverse, c T (ω) = µ−iωη ρ , speeds of sound. K and µ (ξ and η) are the bulk and shear elastic (viscous) moduli, respectively, and ρ the mass density. c T and c L are complex due to finite viscous moduli, i.e., Eq. (2) contains phonon attenuation. Microscopically, such attenuation may be caused by phononphonon, phonon-electron, or phonon-defect scattering. For simplicity, we assume c L (ω)/c T (ω) = √ 3, independent of ω, a good approximation for solids [16,[43][44][45]. c 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 ∼ e iωr ∥ /cT0 /r ∥ (SM), and the damping coefficient reads with ζ ela ≈ 1.29, ζ vis ≈ 1.72, and c ′ T0 the real part of c T0 . Eq. (3) is valid for l ≪ λ, with phonon wavelength λ = 2πc ′ T0 /ω. For typical speeds of sound and l in the nanometer range, the corrections to Eq. (3) are small for ω ≲ 10 9 Hz. This number, 10 9 Hz, is large compared to typical frequencies excited by the probe, since they are expected to be in the range ω ∼ 10 2 . . . 10 5 Hz [11,14,15,19,20,[30][31][32]. Additionally, Eq. (3) assumes the quality factor Q(ω) = µ0 ωη0(ω) , 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 ∼ 10 4 [47,48] (see Table I for experimental parameters). Despite Q being large, phonon attenuation is essential for friction since, as shown below, it dominates the behavior in certain regimes.
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 η 0 . We discuss its physical origin below. Their relative weight is the dimensionless viscosity,η which depends on material properties and probe size l. The numbers of Table I 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(ω) = c T0 (ω)/c T1 (ω) and ρ 0 /ρ 1 [50]. It is insightful to start with the limit n = 0 [51], which yields a Dirichlet boundary condition (DBC) of u(r ∥ , −d, ω) = 0 [40,50]. 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 η ≫ 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 ≈ l. Indeed, we find that, for small d, the leading order of Γ is linear in d, (6) Γ vanishes as d → 0, because the coating, when placed on a stiff substrate, supports fewer phonons as d → 0. In the second step in Eq. (6), we used Eq. (5) to replace Γ ∞ , making apparent the mentioned saturation to Γ ∞ at d ≈ l. The dependence of Γ on d is short range, set by probe size l of Eq. (1). Refs. [19,27,32]. b Refs. [19,20,[30][31][32]. c Ref. [49].  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 η 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),ê z ·σ(r ∥ , −d, ω) = 0 with stress tensor σ and surface normalê z [50]. 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 ≈ Ql. Expanding this case for small d yields In this case, the probe excites traveling surface waves, ∼ e iωr ∥ /cT0 / √ r ∥ , 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 η 0 . As the sheet gets more compliant with d → 0, amplitudes of excitations get larger, and formally diverge as d → 0. The second equality of Eq. (7) implies the mentioned saturation length of d ≈ 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 socalled 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(ω) real, and use ρ 0 = ρ 1 [52]. 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 η = 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 → 0, Γ 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 [53] [54]. Γ thus varies between the bulk results of the substrate (d → 0) and of the coating (d → ∞).
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η/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 ≲ l and d ≳ Ql, and traveling surface waves What about the case ofη ≪ 1 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 ≪ λ the curves are similar to Fig. 2 (a). This is because, for n = 0, the coating does not support traveling waves for  5) and (6). For the chosen parameters, the first peak is at d/l = λ/4l ≈ 50.
d ≲ λ (SM). For d ≳ λ, traveling waves are excited, yielding a sharp transition between evanescent and traveling waves at d ≈ λ/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. This discussion allows identification of mechanism regimes, depicted in Fig. 5. The limit of d → ∞ is dominated by evanescent waves forη ≫ 1, and by traveling waves forη ≪ 1. This is indicated in Fig. 5 by using different colors in the inset graphs. The line n =η separates the same for d → 0, i.e., curves that begin with blue or red. This limit can be found from [55], Another separator is n = 1, dividing curves where Γ ∞ is larger than the limit of d → 0 from the opposite. The last curve, n = 3 √η , together with n =η, 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. 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. Fig. 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][14][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 [56]. We note however qualitative agreement with studies on graphene. Experiments [32] and simulations [27] reported friction on freely standing graphene flakes, i.e., n ≫ 1, finding a decrease with increasing number of layers. The same trend is observed for mounted graphene layers (or flakes) [27,[30][31][32]; graphene is much stiffer (µ ∼ 1 TPa [45,57,58]) than metal substrates (µ ∼ 10 GPa [49,59]), so that the case of n ≫ 1 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,60] 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,56] 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) -217133147 via SFB 1073 (Project A01). We thank Richard L.C. Vink for stimulating discussions. In mem-

DAMPING COEFFICIENT
The damping coefficient can be calculated from a Green-Kubo relation [1][2][3][4], Note that in the second equality, we make use of the FDT [4][5][6][7][8][9], It is more convenient to obtain the Green's function in k ∥ = (k x , k y ) space, since the Kelvin-Voigt model ([Eq. (2)] in the main text) becomes an ordinary differential equation. Let us change the spatial integral to an integral over k ∥ space, The above expression for the damping coefficient can be further simplified if the homogeneity of the Green's function is assumed, i.e., G(r + r 0 , r ′ + r 0 , ω) = G(r, r ′ , ω) [10], The homogeneity thus entails k = −k ′ , i.e., Or, arXiv:2205.01151v2 [cond-mat.stat-mech] 23 May 2023 Plugging it into Eq. (4), one arrives at which we use to evaluate the damping coefficient both numerically and analytically.

BOUNDARY CONDITIONS
Regarding the geometry of the system, it has the surface at z = 0 in the xy plane, while extending infinitely in the x and y directions. For z > 0, it is a vacuum, i.e., no phonon can exist. The boundary conditions (BC) at z = 0 are given by [11][12][13][14] while other components are zero. Here, σ(r, ω) is the Cauchy stress tensor, which is defined as Note that the BC on the surface suggest that the displacement field u x (r, ω) is linearly proportional to the position of the probe X(ω), containing the phase information with respect to X(ω), which, in turn, yields a finite damping coefficient. For −d < z ≤ 0, we have the first layer (the coating) of the solid whose properties are characterized by c L0 (ω), c T0 (ω), and ρ 0 . The second layer (the substrate) expands for −∞ < z < −d with c L1 (ω), c T1 (ω), and ρ 1 . The BC at z = −d require that the displacement field and traction are continuous [15], whereê z is the unit vector in the z direction. The Dirichlet and Neumann BC in the main text are the special cases of the above BC.

TRAVELING VS. EVANESCENT WAVE
Here, we show that traveling waves yield the elastic contribution, and evanescent waves the viscous contribution. Since any wave can be expressed by a superposition of plane waves, an ansatz of the transverse motion of the solid ([Eq. (2)] in the main text) can be written as ) is an amplitude of outgoing (incoming) wave. The longitudinal counterpart can be similarly written with q L = (k 2 ∥ − ω 2 /c 2 L ) 1/2 , which then gives us the full solution G = G T +G L . From the ansatz, it follows that if k ∥ ≤ ω/c T , the solution represents a traveling wave in the z direction. Contrarily, if k ∥ ≥ ω/c T , then it decays exponentially in the z coordinate, i.e., an evanescent wave, resulting from an adiabatic deformation (see below).
Notice also that the integral in Eq. (8) runs in k ∥ space, along which there may exist branches or singularities. Consequently, the integral can be divided into many parts according to the branch or/and singular points, where k s,n n ∈ {1, 2, . . . , N } are the branch and singular points of the Green's function. In what follows, we show that the terms running up to k s,N integrate traveling waves resulting in the elastic contribution, and the last term evanescent waves resulting in the viscous contribution. Let us exemplarily consider the transverse mode. The Green's function at z = z ′ (the argument is henceforth omitted) at the limit of d ≫ λ reads, The branch point is located at k s = ω/c T0 . Note that this branch point coincides with the singularity. The integration is thus divided into two parts. The first part integrates traveling waves, where erf(x) is the error function. Expanding it at a small ω and taking the imaginary part, one arrives at Note that one should perform the integral first and then the expansion, since the singularity renders these operations non-commute. Note also that the leading order of this contribution is independent of viscosity η 0 . The second part, on the contrary, takes care of the evanescent wave solution Here, we expanded the integrand for small ω first, since the integral over k ∥ and the limit of ω → 0 commute [4]. Performing the integration leads us to Taking the imaginary part, one finds the friction contribution from evansecent waves Unlike the elastic contribution from traveling waves, the leading order is linearly proportional to the viscosity η 0 .
Putting them together, we find It is noteworthy that the existence of a branch point or a singularity at a finite value of k ∥ directly indicates the existence of a traveling wave. If the Green's function exhibits no such point, then the corresponding wave is purely evanescent.

EVANESCENT WAVES AND LOCAL DEFORMATION
In this section, we show that evanescent waves represent the adiabatic local deformation due to the probe-sample interaction F , i.e., the Hertzian contact. The adiabatic local deformation can be obtained when the inertial term in [Eq. (2)] in the main text vanishes [11,13,14]. Exemplarily, let us consider the semi-infinite solid. The Green's function for the transverse mode reads This expression is precisely the Green's function in the rhs of Eq. (17). Note that since the probe-sample interaction acts parallel to the surface, this deformation represents the local shearing (e.g., puckering). Importantly, the local deformation contributes to the friction if and only if the solid is viscoelastic. This is because it is the imaginary part of the Green's function that gives rise to the damping coefficient, which is non-zero only for viscoelastic solids.

FUNDAMENTAL SOLUTIONS
For the Dirichlet BC, the Green's function at the surface (z = z ′ = 0), for d being the smallest length scale, is given by The leading order has no branch or singular point, meaning that it is evanescent. Since it is constant in k ∥ space, the above expression becomes a delta function in r ∥ space, For the Neumann BC, the Green's function reads One can easily tell it can be either a traveling or evanescent wave due to the existence of the singular point. By finding the real space counterpart, one can easily identify the dimensional nature of the solution. Assuming c T0 = √ 3c L0 , it is where K 0 is the Bessel function of the second kind of the zeroth order, a form of surface (2D) waves. At a large d ≫ λ on the other hand, the Green's function is given by This, again, can be either a traveling or evanescent wave, however, finding the exact expression of the real space counterpart seems difficult. A simpler way to see the dimensionality of the fundamental solution is to only consider the transverse mode, i.e., Eq. (14). In real space, it is given by which is a spherical (3D) wave solution. Upon arriving at these fundamental wave solutions, no assumption onη is made. In fact,η can only be defined after obtaining Eq. (26). This means, mathematically, that [Eqs. (6) and (7)] in the main text are valid for allη.

THE DAMPING COEFFICIENT FOR AN ARBITRARY INTERACTION AT THE LIMITING CASES
Let us consider a general way to formulate the probe-sample interaction F . Such an interaction arises from the underlying pairwise probe-atom interaction g(X, r ∥ , ω) (it can easily be generalized for vectors as well), Expanding the pairwise interaction in u x , one finds If u 2 x (r ∥ , ω) ≪ r ∥ − X is assumed, the first term can be seen as the pairwise interaction in phase with the motion of X, and the second term out of phase. Now one can calculate the damping coefficient similarly as Eq. (2).
Let us consider the case of small d. For the Dirichlet BC, we use Eq. (22) for the Green's function and perform the integral. Because the Green's function is independent of k ∥ for the leading order of d, the damping coefficient is given simply by the integration of the interaction, The linear dependence of the damping coefficient on the coating thickness d is thus universal to any type of interaction. For the Neumann BC, the Green's function is given by Eq. (24), which, for c T0 = √ 3c L0 , reduces to Note that g(X, k ∥ ) is azimuthally symmetric, so that it can be factored out of the integral over θ. The integral over k ∥ can be done by means of the residue theorem due to the singularities at yielding the elastic contribution from traveling waves. The remaining integrals that are associated with the residue theorem contributes to higher order terms in ω (i.e., O(d/λ), see Ref. [4] for the detailed calculation). The resulting damping coefficient for the leading order of ω and d is The above expression is universal to any types of interaction for the leading order of ω. In fact, the elastic contribution is agnostic to the form of interaction, the characteristic of which we named universality in Ref. [4]. To see this, let us consider the locations of the singularities at ω → 0; the singularities are approaching to k ∥ = 0 as ω → 0. This means the Green's function behaves like a delta function, and therefore the interaction has to be evaluated at k ∥ = 0. By Eq. (28), g(X, k ∥ = 0) = F (X) can be found; the damping coefficient does not depend on the form of the interaction. For detailed discussions and derivations of the universality of the elastic contribution, we refer the readers to Ref. [4].

TRANSVERSE VS. FULL SOLUTIONS
In the main text, the cases for arbitrary n are calculated with only the transverse waves considered. In Fig. 1, we present, for a few selected n values, a comparison between the full (transverse and longitudinal) and the transverse wave solutions. Note that in the plot, the damping coefficients calculated from the transverse solution are scaled by a factor which is nothing but the ratio of [Eq. (5)] in the main text to the viscous contribution in Eq. (20). From Fig. 1, it is clear that the transverse waves can capture the essential physics when it comes to the friction calculation.

FINDING THE PEAKS
Forη ≪ 1, the damping coefficient Γ exhibits peaks at d ≈ λ, resulting from resonances of the excited waves. We identify that there are two types of resonances: body waves and surface waves. The former is due to interference of outgoing and reflected waves. The surface wave resonance, on the other hand, is related to the singularity of the Green's function in Fourier space (k ∥ , ω) [16][17][18].
Let us begin with the peak of the surface wave. From Eq. (26), one finds that the position of singularity is at .
The peak position of surface wave is thus .
The peaks of body waves can be found by considering the plane waves, i.e., evaluating the Green's function at k ∥ = 0, from which the peak positions are obtained, Figures 2 (a) and (b) show the peak positions for n = 0 and n → ∞, respectively. The detected peak positions agree with our predictions as shown in the subsequent plots. Note that in Fig. 2 (b) we do not observe the plateau from the viscous contribution before the peaks. This is becauseη and Q are inversely related to each other, resulting in lQ ≈ λ; the peaks occur before the first plateau of the viscous contribution fully establishes.

DEFINING THE REGIMES IN [FIG. 5] IN THE MAIN TEXT
The scaling behavior of Γ can be understood by studying lim d→0 Γ and Γ ∞ , Keep in mind that ζ ela /ζ vis ≈ 1. The numerator is the damping coefficient of the substrate and the denominator that of the coating. To begin with, let us consider the numerator, from which a separating line n =η (40) can be defined. If n ≫η, the elastic contribution of the substrate is much larger than the viscous contribution of the substrate. Another line can be found when comparing the elastic contribution of the substrate with the viscous contribution of the coating, n = 3 η.
If 3 √η ≲ n ≲η, one finds ζ visη ≲ ζ ela n 3 ; the viscous contribution of the coating is smaller than the elastic contribution of the substrate. Because the latter decays at a slow rate d ≈ Ql, it is visible as shown in [ Fig. 3] for n = 100 in the main text. Another case is when 3 √η ≲ n ≲ 1. The elastic contribution of the substrate is larger than the viscous contribution of the coating. In this case, the damping coefficient remains constant until d ≈ λ.
Comparing the viscous contribution of the substrate to the elastic contribution of the coating defines yet another line, A case of interest would be 1 ≲η ≲ n −2 , where the viscous contribution of the substrate is much larger than the elastic contribution of the coating. This should result in two plateaus, one from the viscous contribution of the coating at d ≈ Ql, the other from the elastic contribution at d ≈ λ similar to those seen in Fig. 2 (a). However, this regime is inaccessible since Q becomes inevitably large for smallη so that Ql becomes comparable to or larger than λ as discussed above. The line defined by Eq. (42) is thus not shown in [Fig. 5] in the main text.