Roughness Evolution Induced by Third-Body Wear

Surface roughness is a key factor when it comes to friction and wear, as well as to other physical properties. These phenomena are controlled by mechanisms acting at small scales, in which the topography of apparently flat surfaces is revealed. Roughness in natural surfaces has been reported to conform to self-affine statistics in a wide variety of settings (ranging from earthquake physics to micro-electro-mechanical devices), meaning that the height profile can be described using a spectrum where the amplitude is proportional to its wavelength raised to a constant power, which is related to a statistical parameter named Hurst exponent. We analyze the roughness evolution in atomistic surfaces during molecular dynamics simulations of wear. Both pairs of initially flat and initially rough surfaces in contact are worn by a third body formed by particles trapped between them during relative sliding. During the first sliding stages, the particles trapped between the first bodies scratch the surfaces. Once the former becomes coated with atoms from the latter, the wear process slows down and becomes “adhesive like.” The initial particle sizes are consistent with the minimum size to be expected for the debris, but tend to grow by material removal from the surfaces and to agglomerate. We show that, for the particular configurations under consideration, the surface roughness seems to converge to a steady state characterized by Hurst exponent close to 0.8, independently of the initial conditions. Supplementary Information The online version contains supplementary material available at 10.1007/s11249-024-01833-9.


Introduction
Understanding the geometry and evolution of rough surfaces is an active research endeavor in tribology (Persson et al., 2004;Bonamy and Bouchaud, 2011;Renard et al., 2013;Candela and Brodsky, 2016;Aghababaei et al., 2023) as roughness mediates friction, wear and lubrication (Godet, 1984).Since the pioneering work of Bowden and Tabor, the response of sliding surfaces is known to depend on the real contact area between the surfaces (Bowden and Tabor, 2001), which, due to their roughness, is but a small percentage of the apparent contact area (Dieterich and Kilgore, 1994).In natural surfaces (Renard et al., 2013;Candela and Brodsky, 2016), across a wide range of scales, roughness appears to follow a fractal distribution (we refer in this case to "self-affine surfaces") which can be characterized in terms of a power law of the wavelengths whose exponent relates to the so-called "Hurst exponent", see Mandelbrot and Van Ness (1968).A number of factors have been put forward to explain how these particular statistics arise in nature: material heterogeneity (Sundaram et al., 2012), plastic mechanisms (Irani and Nicola, 2019;Hinkle et al., 2020), fracture (Bouchaud and Bouchaud, 1994), and corrosion (ASME, 2013), among others.Hinkle et al. (2020) proved that fractal roughness can arise from inelastic deformation induced by simple compression, combined with material discreteness and heterogeneity.
Recently, molecular dynamic (MD) models with simplified potentials have enabled in silico experiments in which the transition from wear regimes (from asperity plastic smoothing to fracture-induced debris creation) can be observed (Aghababaei et al., 2016(Aghababaei et al., , 2017;;Garcia-Suarez et al., 2023).In this context, roughness evolution of 1D surfaces was studied (Milanese et al., 2019(Milanese et al., , 2020)), reporting the generation of self-affine surfaces starting from contacting asperities that are subsequently sheared during surface relative sliding, which gives rise to a third body and whose rolling ends up wearing the surfaces through a process of tearing of shallow clumps of atoms (Milanese and Molinari, 2020), leading eventually to the attainment of a steady-state roughness with fractal characteristics.2022) performed 3D simulations, starting from a configuration with two surfaces in relative sliding motion with pre-formed third bodies in between (Figure 1).This work considered both surfaces that were initially rough and flat as well as third-body particles simultaneously, hence a setting that departs starkly from idealized configurations (Sorensen et al., 1996;Mo et al., 2009;Stoyanov et al., 2013;Eder et al., 2015;Yang et al., 2016;Sharp et al., 2017;Aghababaei, 2019).Among the many insights provided by these simulations, the appearance of cylindrical rolling particles as those observed in experiments (Zanoria and Danyluk, 1993;Zanoria et al., 1995a,b) and the attainment of an apparent steady-state roughness regime stood out.The steady-state is reached after substantial material transfer from the surface to the debris particles.This also translates into decimation of the larger topographical features (those associated to the wavelength of the order of the diameter of the debris); this also implies an apparent "flattening" of the surface that, unlike previous results (Sorensen et al., 1996;Spijker et al., 2011;Stoyanov et al., 2013), is not only associated to atomistic mechanisms but also to debris creation and its coating by surface atoms upon subsequent sliding.
Even though the original paper (Brink et al., 2022) was focused on quantifying frictional forces and surface wear, the surface topography state was extracted at regular increments of the simulations (see Methods section).This manuscript presents the post-processing of the roughness and characterizes the steady-state regimes quantitatively.Section 2 briefly reviews the computations' specificities and presents the topography spectral analysis techniques to be utilized.The surface changes as well as the corresponding Hurst exponent evolution are reported in Section 3 and discussed in Section 4. Conclusions and future work directions are presented in Section 5.

Data generation
Rough surfaces were generated by MD simulations of rigid particles rolling on surfaces made out of a silicon-like model material (Brink et al., 2022).We used LAMMPS (Thompson et al., 2022) with a modified (Holland and Marder, 1998a,b) Stillinger-Weber (Stillinger and Weber, 1985) potential using GPU acceleration (Brown et al., 2011;Brown and Yamada, 2013).While the modified potential does not reproduce the properties of silicon, it can model the brittle fracture (cleavage without dislocation activity) with low computational cost (for details see Brink and Molinari, 2019;Brink et al., 2022).
We used the setups shown in Figure 1.For two simulations (a-b), we prepared a nanocrystalline material with grain size of 3 nm, obtained by the Voronoi tesselation method (Derlet and Van Swygenhoven, 2003).For one simulation (c) we used a single-crystalline material with the (100) surface showing in z direction.Initially, the top and bottom bulk regions (first bodies) each had a size of around 54 × 54 × 20 nm 3 (total of around 8 million atoms in the final simulation cell).For one nanocrystal and for the single crystal, we started from synthetic rough surfaces (Wu, 2000) which were generated with a Hurst exponent of 0.8, a lower wavelength cutoff of 0.5 nm, an upper wavelength cutoff of 27 nm, no roll-off, and an RMS of heights of 2 nm using the software Tamaas (Frérot et al., 2020).The other nanocrystal had a flat surface initially.We then introduced four rigid particles each into the gaps between the surfaces (third bodies).These particles were polyhedra (rhombicuboctahedral shape), except for case (a), where we used round particles for comparison.Their diameter was chosen at ≈ 16 nm.Note that the particles get coated quickly by material picked up from the surfaces.Therefore (i) their initial shape does not matter and (ii) the simulation resembles more closely the adhesive wear case than the abrasive one (Brink et al., 2022).
The sliding simulations were performed with periodic boundary conditions along x and y.A layer of thickness 0.4 nm was fixed at the ends of the top and bottom surfaces, where a normal force of 7.69 µN was applied, corresponding to an average pressure of approximately 2.6 GPa or 8% of the hardness.Next to these boundary layers, another 0.4-nm-thick layer was used each to apply Langevin thermostats at room temperature with a damping constant of 0.01 ps.The center-of-mass velocity of the layer was subtracted from the thermostating calculation to avoid an artificial drag force.
Sliding was imposed on the top first body with a velocity of 20 m/s to a total sliding distance of 1 µm at an angle of 8.5 • off the x direction to avoid that the particles wear the same trench over and over again.
Every 0.1 µm sliding distance (5 ns), the simulation state was recorded.For post-processing, the third bodies are fixed in place while the surfaces are separated (Brink et al., 2022).The atoms on the surfaces of the first bodies were indentified using a surface mesh generation algorithm (Stukowski, 2014) implemented in Ovito (Stukowski, 2009).This algorithm is based on testing if a virtual probe sphere of radius 0.385 nm can penetrate the material or not (Stukowski, 2014).The resulting surfaces are then analyzed as described in the following.

2D surfaces
The height of a surface is defined by a function h Periodicity allows analyzing the surface using Fourier series and, since we work with discrete datasets, discrete Fourier transform.The first step is to interpolate among the original irregular mesh points to then evaluate the interpolant on a N × N regular mesh, N being the number of points in either direction.We use a piece-wise constant interpolant between atoms.See that by the end of this procedure one may obtain a set of points that do not reflect the periodicity of the original point cloud; this motivates the use of "windowing" discussed later.
Heights are known in discrete fashion.The position of the i-th surface atom comes given by a triple (x i , y i , z i ), z i = h(x i , y i ) being the height.The set of all points forms an unstructured mesh which must be re-sampled into a regular grid (sampling intervals ∆x = ∆y = 1 Å) before using discrete Fourier transform methods.Then, the height of the topography features can be expressed as: where the wavenumbers appear, horizontal q x = 2πn/L x and vertical q y = 2πn/L y , taking possible values indexed by n ∈ [0, . . ., N − 1].The amplitude corresponding to each combination of wavenumbers is computed as ĥqx,qy = x,y The heights are rescaled beforehand to guarantee , what amounts to ĥ0,0 = 0.The set of all Fourier coefficients will be referred henceforth as "the spectrum of the surface", and any individual coefficient as "a harmonic".

The h 2
rms is an important parameter when it comes to test if a surface is "self-affine".This parameter represents an average squared height (Jacobs et al., 2017), and thus it conveys the magnitude of the topography oscillations.The 2D power spectral density (PSD) is a function of the wavenumbers defined using the harmonics' amplitudes which is equivalent to the magnitude of the Fourier transform of the height-to-height autocorrelation function (Jacobs et al., 2017), a consequence of Parseval's theorem (Evans, 2010).
The lack of periodicity associated to discreteness and interpolation can introduce spurious high-frequency oscillations in the spectrum of the surfaces (Jacobs et al., 2017).To avoid this issue, windowing is used.In this text, we use radial Hahn window, implicitly assuming that the roughness we are dealing with is isotropic.The radially-symmetric Hahn window is defined as (Jacobs et al., 2017) defined like this for x 2 + y 2 < min(L x , L y )/2 and equal to zero everywhere else.The modified "windowed" heights are given by h Hahn (x i , y i ) = w(x i , y i )h(x i , y i ) for all (x i , y i ) such that x 2 i + y 2 i < min(L x , L y )/2, and h Hahn (x i , y i ) = 0 otherwise.
The Fourier transform comes given in terms of the horizontal wavenumber q x and the vertical one q y .If the surface is isotropic, then the coefficients of the Fourier series depend on the wavenumbers through q r = q 2 x + q 2 y , meaning that the spectral amplitude C 2D must possess axial symmetry with respect to the origin of the q x − q y plane.Thus, for any fixed q r , we can define the radial average of 2D PSDs (implicitly assuming isotropy) as where N θ is a number of angular probes.We probe at θ ∈ [0, 2π/100, . . ., 198π/100], i.e., we average over N θ = 100 points for every fixed q r .
The PSDs that we obtain seem to reasonably satisfy this assumption.Thus, the 2D spectrum indexed by q x and q y is converted into a 1D one that depends on q r .If the surface is isotropic and self-affine, then the spectral amplitudes must scale as ∼ q −2(1+H) , where H is the Hurst exponent.
The value of H is obtained through the slope of the line fitted using logarithmic scales, discarding the roll-off phase (Jacobs et al., 2017).

1D line scans
We have also performed 1D scans on the surfaces.Their theory is briefly introduced next, for further details see Milanese et al. (2019) and the appendix A of Jacobs et al. (2017).
In the 1D case, given a height 1D scan h 1D (x) along the x-direction, the PSD (per unit length, at a wavelength q n ) of self-affine surfaces comes given as which has been discretized using a regular step ∆x = L x /N , thus x k ∈ [0, ∆x, . . ., (N − 1)∆x].Therefore, the discrete spectrum of the PSD contains discrete wavenumbers q n = 2πn/L x for n = 0, 1, 2, . . .Thus which, after averaging over many scans (Jacobs et al., 2017), must satisfy ∼ q −(1+2H) if the roughness is self-affine.
Results of 1D analyses must be averaged across many scans to render the results consistent with the 2D results (Jacobs et al., 2017).We use ten scans along x-direction and ten more along y.

Visualization of roughness evolution
For the geometrical setting and silicon-like material described in the Section 2.1, we show schemes of the roughness evolution in three different cases.The first, Figure 2, corresponds to the surfaces (top and bottom) that are initially flat, and whose bulk contains grain boundaries.The second one, Figure 3, features an initially-rough isotropic surface (bottom one) created with Tamaas (Frérot et al., 2020) with an initial Hurst exponent of 0.8 whose bulk material also contains grain boundaries.Finally, supplementary material Figure A.1.,is similar to Figure 3, but the bulk material is crystalline.
Every plot is accompanied by a scale to measure the amplitude of oscillations with respect to the mean.See that the color code remains the same, but the range of the scales changes between surfaces, since the magnitude of the topographical features evolves.The absolute position on the mean plane is marked in the vertical axis to better appreciate how the surface level descends as surface atoms are transferred to the coating of the third body.Note that the middle point between mean planes of the surfaces corresponds to the height equal to 0. Lighter colors highlight features that "stick out" of the surface, while darker ones penetrate into the bulk.
"Trenches" associated to the scratching by the third body are observed in both top and bottom surfaces after sliding by 0.1 µm (see second row in Figure 2).Bear in mind that the particles' trajectories wind over the whole surface, owing to the sliding direction being not aligned with either axis of the surface.As sliding progresses, the topography "homogenizes": the initially-flat surfaces become isotropically rough (third and fourth row in Figure 2) and the initially-rough ones evolve to a new state, similarly isotropic but characterized by lower height amplitudes (e.g., Figure 3).

2D PSDs
We also provide graphics of PSD evolution in the supplementary material.They clearly reveal "textured" surfaces during the first sliding stages, meaning that we see light spots in figures in supplementary material section C, indicating large magnitude oscillations that break the radial symmetry required by isotropy.Physically, this marked wavelengths are related to the deep grooves left by the wear particles during the first runs, which scratch the initial surface in an abrasive manner.As more sliding unfolds, the debris is coated with surface atoms, changing the wear mechanism to tearing of small flakes of material, induced by adhesion and the debris' rolling movement (Milanese and Molinari, 2020).The texture then fades away, yielding a PSD wavenumber distribution that seems reasonably angle independent, i.e., isotropic, see the final sliding PSD figures in supplementary material section C.This result is reassuring insofar it backs our subsequent analysis as to self-affinity.

Logarithmic slope of radial PSD (C iso ): fitting the Hurst exponent
As explained in the Methods section, isotropic surfaces depend solely on the modulus of the wavenumber vector (q x , q y ), not on the ratio q x /q y .This motivates the definition of a radial PSD, characterized by C iso , eq. ( 5).The Hurst exponent can be easily fitted when this function is expressed in log-log scales.This exercise has been pursued in supplemenrary material section B.

Hurst exponent evolution
Figures 4a to 4c show the roughness evolution corresponding to the cases presented in Figures 2, 3 and A.1.(supplementary material), respectively.The method to extract the exponent from the discrete data was outlined in Section 2.2.The initial values (before any sliding) are found to be, as expected, 0 in the case of initially-flat surface and approximately 0.8 in the two initially-rough cases.2D-analysis results represent top-andbottom average, while 1D ones correspond to averages over top and bottom surfaces, and twenty line scans each, ten along the x direction and ten more along y.

1D line scans
Since the range of data for the PSDs does not span many decades, it is necessary to verify the self-affinity using further measures in addition to C iso .Averaged 1D line scans of the final surfaces match reasonably well with the ones obtained directly from processing the complete 2D surfaces.See Table 1 for average Hurst exponents from 2D analysis (average between the values inferred from the top and bottom surfaces) and the 1D analysis (averaged over twenty scans, ten along the x direction and ten along y).
We provide plots of such results in supplementary material.Moreover, these scans yield an important piece of information: the Hurst exponents obtained from the 1D PSD predict well the ones of the height-to-height correlation function.This provides a consistency in self-affinity that palliates the lack of roughness information over many decades (Milanese et al., 2019), as would be desirable.Supplementary material figures, section B, show the fitting of the Hurst exponent from the 1D PSD on the left panel, while the right one shows the height correlation and the self-affine slope presumed from the Hurst exponent obtained from the PSD.We acknowledge that the exponent obtained from PSD analysis matches well the slope in the short correlation lengths, before the curves level.These plateaus already occur for height differences above the order of 1 nm, but this is simply a result of the low h rms values, which in turn are caused by the depletion of the bigger topography features during the wear process: correlations are only possible for features that actually exist, i.e., up to the order of h rms .

Discussion
The main insight we can extract, c.f. Milanese et al. (2019), it is that, for this material and system size and configuration, the topography of the three surfaces seems to converge to a steady state, characterized by H ≈ 0.8, independently both of the initial conditions (rough or flat) and of the microstructure of the bulk material.The latter remark is in agreement with the findings of Hinkle et al. (2020): material heterogeneity can have a strong influence on wear at this scale (Wattel et al., 2022), but it cannot be the controlling factor of roughness evolution; rather, discrete deformation mechanisms bear responsibility (Irani and Nicola, 2019;Hinkle et al., 2020); in our case, they amount primarily to wear by tearing of flakes of material (Milanese et al., 2020;Milanese and Molinari, 2020).
It is also remarkable how the initially-flat surface, Figures 2 and 4a appears to converge faster to the steady-state roughness regime, with smaller oscillations around a mean slightly greater than 0.8.This seems to indicate that the "memory" of the previous roughness in the other two cases takes longer to be erased.This is reminiscent of the memory length-scale that appears in rate-and-state friction laws (Scholz, 2019).In this context, the frictional state of the interface changes dynamically, reaching the new one after a transient.The extent of this transient is considered a function of the existing microcontacts, i.e, of the roughness.The numerical results seem to reflect this, since the "microcontact population" of the flat surface is very different from the one of the intially-rough surface.Micromechanically, this could be related to the need of attaining an intermediate "indifferent" roughness state (corresponding to H ≈ 0.5, in which the probability of the vertical position of the next atoms is equally probable to be below or above the current one, i.e., the roughness follows a standard random walk).A similar observation is made by Hinkle et al. (2020).From the prior samples, it seems that the roughness may need to oscillate around these values before reaching the steady-state.In order to evolve into the intermediate regime, the features associated to the roughness H ≈ 0.8 have to be erased by wear, while in the flat case they are directly created by wear.
As mentioned in the introduction, simulations involving more ductile materials (either zinc, copper or aluminum) converged to surface welding in which both surfaces were joined together, plastically deforming to engulf the particles in their midst.For this brittle materials, the destiny may be the same in the long term: we acknowledge an ever-growing volume of the coated third bodies (Brink et al., 2022), which could ultimately agglomerate and bridge the gap between surfaces forming the aforementioned shear-band state.Hence, this steady-state could be conceived as an "intermediate asymptotics" state, that may seem locally stable but that can devolve into a shearband-like state over longer timescales.In practice, this final state may be avoided by other mechanisms acting on the said extended time spans.One of such is passivation: the coating atoms may react with the atmosphere and form in turn a composite layer that prevents the flake tearing process, thus deactivating the welding.

Final remarks
The roughness evolution induced by third-body wear has been studied using large-scale molecular dynamics simulations.Using conventional post-processing techniques (Jacobs et al., 2017), the self-affinity of the resulting surfaces has been verified, yielding in addition evidence as to the existence of convergence to roughness with H ≈ 0.8.This final state, induced under certain circumstances (e.g., brittle enough material and abrasive wear particles), appears to display similar roughness characteristics independently of both the initial topography (either flat or rough with a higher surface roughness H ≈ 0.8) and bulk lattice structure (either single crystal or presence of grain boundaries).
We stress that this independence from bulk structure is consistent with Hinkle et al. (2020), and that these authors also reported the existence of a transient "random-walk" state (H ≈ 0.5) during the evolution of roughness prior to attainment of a steady state.We also emphasize that Hinkle et al. (2020) run purely plasticity-driven simulations, wherein neither third bodies nor fracture were present.

Figure 1 :
Figure 1: Setup of the sliding contact simulations.In all cases, two blocks of bulk material (yellow: bottom block, red: top block) were put in contact with rigid wear particles (gray).Blue areas indicate the boundaries where force and displacement were imposed.We used different initial setups: (a) Nanocrystalline first bodies (grain boundaries indicated by black and gray atoms) starting from a flat surface and round wear particles, (b) nanocrystalline first bodies starting from artificial surface roughness and polyhedral wear particles, and (c) the same as (b) but with single crystalline first bodies.

Figure 2 :
Figure 2: Snapshots of surface evolution (units in Å unless otherwise stated): silicon surfaces (top and bottom) initially flat.Corresponding Hurst exponent evolution in Figure4a.Note changing scales.Vertical axes mark the (evolving) mean height of each surface.In either surfaces, white means topographical features "bulging out", while blue means penetrating into the surface bulk.

Figure 3 :
Figure 3: Snapshots of surface evolution (units in Å unless otherwise stated): silicon surfaces (bottom) initially rough, nanocrystalline bulk.Corresponding Hurst exponent evolution in Figure 4b.Note changing scales.Vertical axes mark the (evolving) mean height of each surface.In either surfaces, white means topographical features "bulging out", while blue means penetrating into the surface bulk.

Figure 4 :
Figure 4: Evolution of the Hurst exponent of Si surfaces worn by wear particles during sliding: three configurations.

Table 1 :
Comparison final Hurst exponents, derived from two different methods.