# Breakdown of Reye’s theory in nanoscale wear

Joaquin Garcia-Suarez a, Tobias Brink b and Jean-François Molinari a
aCivil Engineering Institute, Materials Science and Engineering Institute,
École Polytechnique Fédérale de Lausanne (EPFL), CH 1015 Lausanne, Switzerland
b Max-Planck-Institut für Eisenforschung GmbH,
Max-Planck-Straße 1, 40237 Düsseldorf, Germany

###### Abstract

Building on an analogy to ductile fracture mechanics, we investigate the
energetic cost of debris particle creation during adhesive wear .
Macroscopically, Reye proposed in 1860 that there is a linear relation between
frictional work and wear volume at the macroscopic scale. Earlier work
suggested a linear relation between tangential work and wear debris volume
also exists at the scale of a single asperity, assuming that the debris size
is proportional to the micro contact size multiplied by the junction shear
strength. However, the present study reveals deviations from linearity at the
microscopic scale. These deviations can be rationalized with fracture
mechanics and imply that less work is necessary to generate debris than what
was assumed. Here, we postulate that the work needed to detach a wear particle
is made of the surface energy expended to create new fracture surfaces, and
also of plastic work within a fracture process zone of a given width around
the cracks. Our theoretical model, validated by molecular dynamics
simulations, reveals a super-linear scaling relation between debris volume
($V_{d}$) and tangential work ($W_{t}$): $V_{d}\sim W_{t}^{3/2}$ in 3D and
$V_{d}\sim W_{t}^{2}$ in 2D. This study provides a theoretical foundation to
estimate the statistical distribution of sizes of fine particles emitted due
to adhesive wear processes.

_K_ eywords Adhesive wear $\cdot$ Ductile fracture $\cdot$ Plasticity $\cdot$
Debris volume $\cdot$ Frictional work

##  1 Introduction

Adhesive wear is an unavoidable phenomenon at contacting surfaces subjected to
strong adhesive bonds (Rabinowicz, 1995, Burwell and Strang, 1952). It occurs
due to microscopic (adhesive) contacts that form wear particles during sliding
(Burwell and Strang, 1952, Archard, 1953, Rabinowicz, 1958), which are a
result of the surface roughness at small scales (Dieterich and Kilgore, 1994,
Renard et al., 2013, Bowden et al., 1939). In general, the severity of wear is
a function of the applied normal force and is therefore connected to friction
(Burwell and Strang, 1952, Archard, 1953, Rabinowicz et al., 1951), but a
parameter-free generally applicable model was not yet found (Meng and Ludema,
1995, Collins, 1993). In light of increasing environmental (Baensch-
Baltruschat et al., 2020, Grigoratos and Martini, 2015) and health (Kole et
al., 2017) concerns related to fine particle emissions, a mechanistic
understanding at the level of single wear particle formation is needed (Vakis
et al., 2018, Renouf et al., 2011).

Different mechanisms have been put forward to explain adhesive wear. These
include wear debris formation with experiments dating back to Archard
(Archard, 1953) and observed in many cases (Bhushan and Sundararajan, 1998,
Chung and Kim, 2003, Liu et al., 2010, Greenwood and Tabor, 1955, Brockley and
Fleming, 1965), plastic deformation of contacting asperities, which we may
call Holm’s mechanism (Holm, 2013), and more recently atom-by-atom attrition
(Gotsmann and Lantz, 2008, Bhaskaran et al., 2010, Sato et al., 2012, Jacobs
and Carpick, 2013, Stoyanov et al., 2014). Atomistic simulations (Stoyanov et
al., 2014, Sorensen et al., 1996, Zhong et al., 2013) generally show plastic
deformation, but not Archard’s debris formation mechanism. However, recent
work (Aghababaei et al., 2016) could reconcile those observations and revealed
that a transition between plastic deformation and debris formation is governed
by a critical length scale $d^{*}$. The parameter $d^{*}$ describes a minimum
junction between contacting asperities for the creation of a wear particle and
is related to material properties:

$\displaystyle d^{*}=\lambda\frac{\Delta\mathrm{w}}{\tau^{2}/2\mu}\,,$ (1)

where $\lambda$ is an order-1 prefactor that encapsulates the influence of the
geometry, $\tau$ is the shear strength, $\mu$ is the material’s shear modulus
and $\Delta\mathrm{w}$ is the per-crack decohesion work, which is in general
equal to the critical energy release rate $G_{Ic}$ for crack opening in mode
I, but can be replaced by two times the surface energy $2\gamma_{s}$ in the
case of brittle materials. The insight provided by this length scale for
adhesive wear has been used to elucidate nanoscale friction (Barras et al.,
2021, Brink et al., 2022), wear particle formation (Aghababaei et al., 2017,
Frérot et al., 2018, Brink et al., 2021), and surface morphology evolution
(Milanese et al., 2019, 2020).

The connection between friction and wear has always been intriguing. Archard
(Archard, 1953) hypothesized that the macroscopic wear volume created during
relative sliding of two surfaces is a linear function of the normal force
(similar to the tangential force in Coulomb friction) times the sliding
distance and inversely proportional to the material’s hardness. Reye (Reye,
1860) postulated a linear relation between the tangential work $W_{t}$ and the
total wear volume at the macroscopic level, which is also often found in
experiments (Rabinowicz, 1995, Burwell and Strang, 1952, Whittaker, 1947, Uetz
and Föhl, 1978, Fouvry et al., 2003, 2001). At the nanoscale, a similar linear
relation between a debris particle’s volume $V_{d}$ and the tangential work
$W_{t}$ required to form it was reported (Aghababaei et al., 2017). This
linearity was anchored on four assumptions:

  1. #1:

Bowder and Tabor’s frictional force argument (Bowden and Tabor, 2001): the
necessary peak force to overcome the “frictional force” arising from a
microcontact is approximately equal to the area of the contact times the shear
strength of the material, i.e., $F_{max}\approx\tau A_{contact}$. As, by
supposition, this contact leads to debris creation, the area of contact can
also be interpreted as a cross-section area of the wear particle
$A_{contact}\approx A_{debris}$, and thus it can be related to the junction
size (characteristic length of the contact patch), $d$. Combining everything:
$F_{max}\approx\tau A_{contact}\sim\tau d^{2}$.

  2. #2:

Effective sliding distance: in order to form a debris particle, the two
surfaces must slide relatively over a distance equal or close to the junction
size, $S_{eff}\approx d$.

  3. #3:

Tangential work: the work necessary to create the particle, $W_{t}$, is
approximately equal to the effective sliding distance times the peak force,
$W_{t}\approx F_{max}S_{eff}$.

  4. #4:

The volume of the debris particle is proportional to the junction size:
$V_{d}\sim d^{3}$. Moreover, assuming that the prefactor connecting the two
quantities is close to unity, one reaches $V_{d}\approx d^{3}$.

Thus, mathematically,

$\displaystyle\int F_{t}ds=W_{t}\approx F_{max}S_{eff}\approx\tau
d^{3}\approx\tau V_{d}\,,$ (2)

where $F_{t}$ represents the tangential force, leading to the proposed volume
estimator $\hat{V}_{d}=W_{t}/\tau$. The salient feature of this model is that
it predicts a linear scaling between debris volume and frictional work.
However, despite $\hat{V}_{d}$ comparing satisfactorily to simulations when
$d\approx d^{*}$ (Aghababaei et al., 2017), we will show in the following how
the performance of this estimate deteriorates as the size of the asperities
increases. In particular, we will discuss molecular dynamics simulations
results with contact junctions larger than $d^{*}$ that reveal super-linear
scaling, and will present theoretical arguments to rationalize our
observations. This work seeks to provide a more exhaustive description of the
wear process at the single asperity level from an energy-balance standpoint,
leveraging on notions of plasticity and ductile fracture mechanics.

The text is structured as follows. In Section 2, we detail the numerical
models that enabled the simulation results presented in Section 3. These
results show a super-linear scaling of debris size with frictional work.
Section 4 outlines a new framework to resolve the contradiction with the
previous theory that argued for a linear scaling. Further implications are
commented in Section 5, while Section 6 presents the final conclusions.

##  2 Methods

###  2.1 Numerical simulations

(a) (a)

(b) (b)

Figure 1: Evidencing the limitations of the debris volume estimate eq. 2
(Aghababaei et al., 2017) using both 3D simulations. Panel (a): scheme of the
3D simulations, including asperity diameter $D$, contact angle $\theta$ and
loading conditions (imposed velocity at top), model presented in Brink and
Molinari (2019). The two asperities overlap over a circle of diameter
$j=0.75d$ (junction size). Panel (b): comparison of debris volume measured in
silico (black) v. estimate (orange), $\theta=30^{\text{o}}$.

We performed molecular dynamics simulations on different model asperity
geometries in LAMMPS (Plimpton, 1995, Thompson et al., 2022). In three
dimensions (3D), we modeled overlapping spherical asperities as described in
Brink and Molinari (2019). We used a modified Stillinger–Weber (Stillinger and
Weber, 1985) Si potential, with increased bond-angle stiffness (Holland and
Marder, 1998b), which better reproduces the fracture behavior at the cost of
the other material properties (Holland and Marder, 1998a, b). The integration
time step was $1\,\mathrm{fs}$. In order to obtain an isotropic sample, we
produced a glass by melt quenching using the procedure described in Refs.
(Fusco et al., 2010, Brink and Molinari, 2019). The material was found earlier
(Brink and Molinari, 2019) to have a critical length scale of
$d^{*}=18\,\mathrm{nm}$. Then, the geometry sketched in Fig. 1(a) was cut out.
Here, we used asperity diameters of $D=10,15,20,30,40,50,60\,\mathrm{nm}$. The
size of the bulk region in the top and bottom crystal was $103\times 82\times
15\,\mathrm{nm}^{3}$ each. In order to save computational time, we started
from asperities that were already in contact on a circular area with diameter
$d=0.75D$. Two sets of simulations were considered, one where the contact area
between the two asperities is aligned with the sliding direction
($\theta=0^{\circ}$) and one where it is inclined ($\theta=30^{\circ}$). After
equilibration for $100\,\mathrm{ps}$ at $T=300\,\mathrm{K}$, we applied a
normal pressure of $0.8\,\mathrm{GPa}$, imposed a tangential displacement
velocity of $20\,\mathrm{m/s}$ at the top boundary, and kept the bottom
boundary fixed. Langevin thermostats were applied over 4-Å-thick layers next
to the top and bottom boundaries. For the force calculation, the drag force
term of the thermostat was subtracted. The tangential force was computed as
the reaction force at the top boundary.

We estimate the debris volume by counting the number of atoms in the debris
particle and multiplying them by the average atomic volume in the bulk. We
also estimate the shear strength of the material via independent simulations,
finding
$\tau={\color[rgb]{0,0,0}7.9\,\mathrm{GPa}=}49.31\,\mathrm{eV/nm^{3}}$.
Numerical simulations were visualized using OVITO (Stukowski, 2009) and
Mathematica (Wolfram, 2000).

##  3 Numerical observations: deviation from linear trend

We conducted several sliding simulations of asperity–asperity contacts using
the 3D setup. Examples of differently-sized wear particles with
$\theta=30^{\circ}$ are shown in Fig. 2. The critical length scale for this
model material was $d^{*}\approx 18\,\mathrm{nm}$ (Brink and Molinari, 2019),
and asperities smaller than this did indeed plastify instead of emitting wear
particles (not shown here). We estimated the wear volume for simulations with
$d>d^{*}$ by counting the number of atoms per wear particle and extracted the
tangential force from the simulations. Figure 1(b) shows the measured wear
particle size as a function of the asperity diameter. It can be seen that the
volume estimator $\hat{V}_{d}=W_{t}/\tau$ from eq. 2 strongly underestimates
the resulting particle volume for large $d$. The quantitative agreement for
$d\approx d^{*}$ is quite good, however. The snapshots in Fig. 2 already
indicate that the relative amount of plastified material decreases and the
fracture process becomes more and more brittle, which might suggest that the
assumption #1 of the original model ($F_{max}=\tau A_{contact}$) could be
invalid.

It is clear that the numerical results are not linear, but proportional to a
power of the work with exponent greater than one. This means that the root
cause of this disagreement cannot be a faulty estimation of the shear strength
$\tau$, but must be associated to the breakdown of one of the four assumptions
discussed in Section 1.

Figure 2: Plastic strain for the 3D simulations with $\theta=30^{\circ}$. A
slice through the middle of the wear particle is shown to visualize the
plasticity in the bulk. (a) For the smallest asperity diameter close to
$d^{*}=18\,\mathrm{nm}$, significant plasticity occurs in the bulk of the
particle. (b)–(d) With increasing asperity size, the plasticity becomes more
localized and does not fill the whole particle.

##  4 New framework based on ductile fracture mechanics

The numerical results reveal the necessity to extend the current theory. We
note that the details of the crack propagation process are not explicitly
accounted for during the derivation of eq. 2. The theory of Linear-Elastic
Fracture Mechanics presupposes that the strain energy within the body goes
into breaking bonds between atoms, which in turns means that new surfaces are
formed and a crack propagates. Griffith’s original energy balance argument
(Griffith and Taylor, 1921) hinges on the assumption that the plastic
dissipation occurring due to the stress concentration at the crack tip is a
small percentage of the total energy being dissipated. To quantify the amount
of plasticity accompanying the fracture process, the plastic radius around the
crack tip, $r_{p}$ (Janssen et al., 2004), is used. This parameter is
traditionally presented as

$\displaystyle r_{p}$
$\displaystyle={1\over\pi}\left({K_{I}\over\sigma_{\mathrm{y}}}\right)^{2}\,,$
(3a) where $K_{I}$ is the mode I stress intensity factor, which reaches its
critical maximum value at the initiation of crack growth,
$K_{I}=K_{Ic}=\sqrt{G_{Ic}E^{\prime}}$ where $E^{\prime}=E$ (Young’s modulus)
in plane stress and $E^{\prime}=E/(1-\nu^{2})$ ($\nu$ being Poisson’s ratio)
in plane strain, and $G_{Ic}$ is a material parameter termed “critical energy
release rate” (in mode I). Thus $\displaystyle r_{p}|_{max}$
$\displaystyle={1\over\pi}{G_{Ic}E^{\prime}\over\sigma_{\mathrm{y}}^{2}}\,,$
(3b)

defines the characteristic size of this plastic region. The condition of
small-scale yielding (SSY) sets the range of validity of brittle fracture in
terms of the plastic radius and a characteristic geometric length of the
system. A generally agreed-upon test is checking if $L\geq 25r_{p}$
(Hutchinson, 1983) ($L$ being any of the characteristic geometrical in-plane
lengths involved, in this case we can take $L=d$). If so, the use of the
brittle approximation is warranted and crack plasticity can be ignored. If
that is not the case, plasticity can be taken into account by modifying the
brittle fracture criterion: in the brittle case, we had $G_{Ic}=2\gamma_{s}$,
while in the ductile one $G_{Ic}=2(\gamma_{s}+\gamma_{p})$, where the material
parameter $\gamma_{p}$ represents the plastic energy dissipation per unit of
new surface area.

Notice that $r_{p}$ and $d^{*}$ scale the same but differ in their prefactors
due to the difference in geometry and loading mode (tension versus shear). In
the case of single-asperity wear, $L$ is taken as the junction size $d$, and
SSY will not hold when $d/d^{*}\sim d/r_{p}|_{max}\ll 25$.

There is yet another parameter that appears in this context (Pineau and
Pardoen, 2007), the fracture process zone (FPZ) around a crack tip. It
displays the same scaling in terms of the mechanical properties, i.e.,
$d^{*}\sim\ell_{FPZ}$ too in the case of brittle materials :

$\displaystyle\ell_{FPZ}\sim{\Delta\mathrm{w}\over\tau^{2}/\mu}\,.$ (4)

This parameter represents the size of the damaged region that either
eventually nucleates a crack or along which the crack extends; it is
characterized by stress concentrations, micro-crack formation and/or other
degradation processes. The fact that the critical junction size scales in the
same fashion as the fracture process zone suggests a reinterpretation of the
former: stress concentrations around small asperities can only nucleate cracks
if the asperity itself is large enough to host a fracture process zone within
it, thus $d^{*}$ can be thought as characterizing the minimal geometry of the
system formed by interlocked asperities that can fit a FPZ that nucleates a
crack, which leads to third-body formation. Figure 3 represents schematically
three possible scenarios related to the prior discussion, each one of them
arising from changing the material properties while maintaining exactly the
same geometry and scale: panel (a) corresponds to a material with
$\ell_{FPZ}\gg d$ where stress concentrations can only lead to plasticity and
surface smoothing (Holm’s mechanism), in (b) $\ell_{FPZ}\approx d$ so
substantial plasticity gives rise to cracks and eventually leads to third-body
creation, and (c) represents a brittle material with $\ell_{FPZ}\ll d$ wherein
fractures nucleate before inelastic deformation takes place (but plasticity
still appears around the crack path).

Figure 3: Different behavior for the same system (in terms of scale and
geometry) of interlocking asperities, but changing the material. Each color
represents a different material with properties yielding different
$\ell_{FPZ}$, whose value can be interpreted as the diameter of the patterned
red area. (a) System with $d\ll\ell_{FPZ}$: all-encompassing plasticity. (b)
System with $d>\ell_{FPZ}$ while $d/\ell_{FPZ}\approx 1$: plasticity widely
present, but cracks can nucleate around spots of stress concentration (compare
to results from simulations using the phase-field method in Collet et al.
(2020)). (c) System with $d\gg\ell_{FPZ}$: localized plasticity around the
fracture path (compare to phase-field results from Carollo et al. (2019)).

###  4.1 Discussion of previous work in the context of fracture mechanics

One of the premises that led to both eq. 1 and eq. 2 is that the volume that
is plastified during sliding prior to debris creation, $V_{p}$, scales
proportionally to the volume of the region surrounding the contact patch, and
thus is similar to the final particle volume, $V_{p}\sim V_{d}\sim d^{3}$.
This was verified in simulations in which $d\approx d^{*}$ (Aghababaei, 2019),
where substantial inelasticity can be observed via post-processing or just by
looking at the permanent shape changes in the asperities prior to detachment.
This is not surprising: the transition from a state that is plasticity-
dominated to one where fracture also appears does not mean that plasticity is
excluded in the latter. Actually, substantial inelastic deformation can
accumulate before the cracks grow (Aghababaei, 2019) (see intermediate panel
in Figure 3). On the other hand, in the limit of $d\gg d^{*}$, fractures
develop before large deformations and inelasticity can occur.

Resorting again to the 3D simulations, let us visualize the extent of plastic
deformation as the asperity size increases. Figure 2(a) shows the smallest
size ($d\approx d^{*}$), where plasticity penetrates the bulk of the system
formed by the contacting asperities, while in Figure 2(d) traces of inelastic
activity are only found in a narrow region close to the crack path. The latter
hints at a ductile fracture, in which substantial plasticity accompanies the
crack tip trajectory, the volume of the plastified region at the end of the
crack propagation process being $V_{p}\sim d^{2}\ell_{FPZ}$, as the full crack
path area is equivalent to the asperity base ($\sim d^{2}$). Note that this is
a simplification, as the actual shape of the plastified region is more
complex, but we assume that the scaling of the plastic zone size is still
valid. Therefore, we find that $V_{p}\sim V_{d}$ only when
$d\sim\ell_{FPZ}\sim d^{*}$.

In order to understand that, we decompose the tangential work into
$W_{plastic}$ (work gone into permanent deformation anywhere in the asperity)
and $W_{debonding}$ (work invested in bond breaking). Energy dissipated as
heat in other processes can be neglected assuming close to quasi-static
loading conditions. The step-by-step reasoning goes as follows: start from

$\displaystyle W_{t}$ $\displaystyle\approx W_{plastic}+W_{debonding}\,,$ (5a)
so dividing by $\tau V_{d}$ gives $\displaystyle{W_{t}\over\tau V_{d}}$
$\displaystyle={\hat{V}_{d}\over V_{d}}\approx{W_{plastic}\over\tau
V_{d}}+{W_{debonding}\over\tau V_{d}}\,.$ (5b)

A volume around the crack tip is plastified while the crack propagates. We
assume a characteristic local plastic strain $\varepsilon_{p}$ after which the
crack propagates one step further, plastifying a new volume up to a strain
$\varepsilon_{p}$, and so on (cf. Fig. 3(c)). Recall that along the full crack
path a volume of $V_{p}\sim d^{2}\ell_{FPZ}$ will be plastified, so we can
express the total plastic work as $W_{\text{plastic}}\approx\tau
V_{p}\varepsilon_{p}$. Under the further assumption that
$\varepsilon_{p}=\text{const.}$ (which is likely material dependent), we
obtain the scaling

$\displaystyle{W_{plastic}\over\tau V_{d}}\sim{\tau
V_{p}\varepsilon_{p}\over\tau V_{d}}\sim{V_{p}\over
V_{d}}\sim{d^{2}\ell_{FPZ}\over d^{3}}\sim{d^{*}\over d}\,.$ (6)

The first term thus decays as the size increases, consistent with what we have
already seen in Figure 2. With the remaining term we find

$\displaystyle{W_{debonding}\over\tau V_{d}}$
$\displaystyle={2\gamma_{s}A_{created}\over\tau
V_{d}}\sim{2\gamma_{s}\over\tau d}\sim\left({\tau\over\mu}\right){d^{*}\over
d}\,,$ (7a) where $A_{created}$ is the total new area created by the
fractures. The factor inside the parenthesis does not change if the material
remains the same, so if only the size changes
$\displaystyle{W_{debonding}\over\tau V_{d}}$ $\displaystyle\sim{d^{*}\over
d}\,.$ (7b)

Thus, both parts of the tangential work have the same scaling and we reach

$\displaystyle{W_{t}\over\tau V_{d}}$ $\displaystyle={\hat{V}_{d}\over
V_{d}}=\mathcal{O}\left({d^{*}\over d}\right)\,.$ (8)

This result means that the estimate can perform well for the smallest
asperities ($d\approx d^{*}\Rightarrow\hat{V}_{d}\approx V_{d}$), but it may
underpredict the debris volume as the size increases ($d\gg
d^{*}\Rightarrow\hat{V}_{d}\ll V_{d}$).

Figure 4: Verifying the assumptions behind eq.(2). (a) Assumption #1:
Quadratic scaling of peak force with junction size in 3D: results for two
geometries (two contact angles ($\theta$)). (b) Assumption #2: Decreasing
effective sliding as a portion of junction size $d$. (c) Quadratic scaling of
tangential work with junction size.

###  4.2 Insights from a continuum beam model

Figure 5: Analytical model scheme. (a) Abstraction from contacting asperities
into Timoshenko beam system that can be halved resorting to anti-symmetry
conditions (Aghababaei and Budzik, 2020). (b) Idealized loading evolution:
linear up to point defined by $F_{max}$ and $S_{eff}$, followed by almost-
instantaneous drop (in reality, it would have a large slope proportional to
the Rayleigh wave velocity).

We supplement the previous discussion with an analytical continuum model
presented in Aghababaei and Budzik (2020), which helps us ascertain the system
behavior as its size increases and the continuum limit is approached. The
idealized model considers two interlocking asperities as a rectangular cross-
section Timoshenko beam (Timoshenko, 1922) having out-of-plane thickness
$\mathrm{t}$, while the other cross-section length corresponds to the junction
length $d$. The height of each asperity is $h/2$, so the interlocked system
has a total height of $h$, see Figure 5(a).

Assuming low-velocity, displacement-controlled loading conditions, we grant
that most of the elastic energy is released much faster than the time that it
takes for the remote loading to change once the failure conditions are
attained, i.e., we assume a sliding-force relation as in Figure 5(b). This
choice will be discussed at the end of the section.

Each asperity can be considered independently exploiting antisymmetry
conditions (Aghababaei and Budzik, 2020). We consider an imposed displacement
$u|_{y=h/2}$ instead of an imposed force. An estimate of the stiffness $K$ of
the system as a function of crack length $a$ was also provided in Aghababaei
and Budzik (2020) (the cracks of length $a$ are assumed to grow at the base of
each beam):

$\displaystyle\mathrm{K}(a)=\left[{4(h/2)^{3}\over
E(d-a)^{3}\mathrm{t}}+{h/2\over\kappa\mu(d-a)\mathrm{t}}\right]^{-1}\,,$ (9)

where $\kappa$ is the shear coefficient (Timoshenko, 1922), this dimensionless
parameter appears in Timoshenko theory as a correction factor to properly
account for the real distribution of shear stresses in the cross-section. Note
the different scaling of each addend in terms of $h/d$, the aspect ratio of
the asperity. This approximation assumes that the beams are clamped to a rigid
half-space, which is obviously not the case (the connection of the asperities
to the surfaces provides extra compliance); however, since we are primarily
interested in a scaling analysis, this simplification does not represent an
important drawback. Using the simple beam model, the total energy of the
system is twice the energy stored in each beam:

$\displaystyle E_{total}=\mathrm{K}(a)u^{2}|_{y=h/2}\,,$ (10)

thus the critical sliding $u=S_{eff}$ that triggers crack propagation can be
computed as the displacement at which the critical energy release rate
$G_{Ic}$ is attained:

$\displaystyle G={1\over\mathrm{t}}{\partial E_{total}\over\partial
a}=-{S_{eff}^{2}\over\mathrm{t}}{\partial\mathrm{K}\over\partial a}=G_{Ic}\,.$
(11)

If we further assume, on top of low aspect ratio, that the asperities are
initially defect-free (so $a/d\to 0$) and that unstable crack growth happens
immediately after $G_{Ic}$ is attained, see Figure 5 panel (b), this yields

$\displaystyle
S_{eff}=\left({G_{Ic}\mathrm{t}\over{\partial\mathrm{K}/\partial a}|_{a\to
0}}\right)^{1/2}=2\sqrt{G_{Ic}h/\kappa\mu}+\mathcal{O}\left[\left({h\over
d}\right)^{2}\right]\,,$ (12)

where the second addend represents terms that will be meaningful only if the
assumption of low aspect ratio was removed.

The total tangential work done over the system (two asperities) to the point
when sudden failure by unstable crack growth happens is

$\displaystyle
W_{t}=\mathrm{K}(0)S_{eff}^{2}=2\mathrm{t}dG_{Ic}=4\mathrm{t}d(\gamma_{s}+\gamma_{p})\,,$
(13)

that is, the energy necessary to grow two cracks over an area $td$ (the cross-
section area) on each asperity. See that if we further assume that $t\sim d$
(the out-of-plane thickness being of the same order as the other cross-section
length) we reach $W_{t}\sim d^{2}$. This result indicates that the energy
invested in sliding leading to debris creation can be controlled by fracture
and be proportional to area, instead of being controlled by plasticity and
proportional to volume.

It is clear that the second assumption eq. 2 rests upon ($S_{eff}\approx d$)
will enter in conflict with the latter result: see the $\sqrt{h}$
proportionality in eq. 12, i.e., the sliding distance to failure is not a
linear function of the size. This is consistent with our numerical
observations, Figure 4 panel (b). Note that the trend does not apply for the
smallest sizes ($d=20\,\mathrm{nm}$), in which the particle does not fully
detach as it sticks to the surfaces (see Figure 2(a) and Figure B.5(a) in the
supplementary material). For the larger sizes, which do lead to neat third-
body formation, the trend is met. See that both the numerical and analytical
results do yield $S_{eff}\approx d$ for the smallest asperities.

Likewise, Figure 4 panel (a) displays the maximum force measured in the
simulations, and reveals a quadratic scaling that appears consistent with the
Bowden and Tabor model (Bowden and Tabor, 2001), even though the continuum
model suggests

$\displaystyle F_{max}=\mathrm{K}(0)S_{eff}\approx\sqrt{G_{Ic}h\mu\kappa}\,.$
(14)

This indicates that the maximum force still seems to be dictated by the
significant plasticity that can be observed at the simulation sizes studied
here (see for example Fig. 2) and depends linearly (2D) or quadratically (3D)
on the contact size. Nevertheless, the sliding distance $S_{eff}$ is in
agreement with the continuum model, which also means that force dropoff after
reaching $F_{max}$ must become steeper for larger the asperities.

In conclusion, given the appraisals obtained from eqs. 12 and 14, the data
shows that the assumption #2 ($S_{eff}\approx d$) used to derive the original
volume estimator in eq. 2 is violated. We expect from the continuum model that
even assumption #1 ($F_{max}\approx\tau d^{2}$) would not be followed for much
larger asperity–asperity contacts, but the data presented in this paper likely
covers a range that is too small to clearly see this transition to the
idealized case used to derive eq. 14. These findings also make intuitive sense
if we think of crack propagation regimes being a function of the asperity
size:

  * •

Small junctions, $d\sim d^{*}$, strength-controlled propagation: the stress
intensity factor is strongly influenced by the surroundings’ geometry, and its
stiffness by extension. The process is driven by the remote sliding condition,
and the strain energy stored during prior deformation (along with the
stiffness of the system) is lost gradually at a rate proportional to the
imposed velocity.

  * •

Large junctions, $d\gg d^{*}$, toughness-controlled: even though the crack
nucleation will depend on the local geometry of the system, most of its growth
will take place far from the model’s edges, being effectively independent of
the geometrical features. The stored energy is released at a fast rate
(related to rapid crack propagation), as assumed in Figure 5(b), while the
crack grows in toughness-controlled conditions.

As the size of the system is increased (greater $d$) from an initial size
$d^{*}$, the relative level of toughness-controlled propagation will in turn
also increase, in detriment of ever smaller portions of strength-controlled.
We further substantiate this point in the supplementary discussion around
Figure B.4.

Finally, let us mention that the same trends reported herein are also observed
in the 2D simulations, Figure A.2: effective sliding as a portion of asperity
size $d$ decreases as the size increases while the peak force does scale
linearly as predicted by Bowden and Tabor’s model in 2D ($F_{max}\propto d$).

###  4.3 A new scaling relation

The system approaches a continuum-like situation similar to the one studied in
the prior section as the size of the junction increases over the threshold
value, hence the energy released at the crack tip becomes the dominant source
of dissipation, i.e.,

$\displaystyle W_{t}\approx G_{Ic}A_{created}\,,$ (15)

where $A_{created}$ represents the total new area, the sum of the two new
surfaces, one in each asperity. This work must be approximately equal to the
strain energy stored during sliding, it does not account for the extra work
that goes into the next stage: “rotating out” the particle.

Extrapolation of the attested scalings $V_{d}\sim d^{3}$ (3D) and $A_{d}\sim
d^{2}$ (2D) (Aghababaei et al., 2017) combined with eq. 15 leads to

$\displaystyle G_{Ic}d^{2}\sim G_{Ic}V_{d}^{2/3}\Rightarrow
V_{d}\sim\left({W_{t}\over G_{Ic}}\right)^{3/2}\text{ in 3D}\,,$ (16a)
$\displaystyle G_{Ic}td\sim G_{Ic}\mathrm{t}A_{d}^{1/2}\Rightarrow
V_{d}=\mathrm{t}A_{d}\sim{1\over\mathrm{t}}\left({W_{t}\over
G_{Ic}}\right)^{2}\text{ in 2D}\,.$ (16b)

Here, $\mathrm{t}$ is the thickness of the 2D system. These scaling relations
are the main contribution of this article. The corresponding pre-factors
require the estimation of the total crack area, which in turn must depend on
the specific geometry of the junction and the asperities (its size and
presence of the stress concentration spots).

Assuming that both the prefactors (called $\mathrm{k}$ in both cases) and the
fracture toughness were known, we would write

$\displaystyle\hat{V}_{d}$
$\displaystyle=\mathrm{k}\left(\frac{W_{t}}{G_{Ic}}\right)^{1.5}\,.$ (17a) for
the 3D case, and the 2D one $\displaystyle\hat{V}_{d}$
$\displaystyle={\mathrm{k}\over\mathrm{t}}\left({W_{t}\over
G_{Ic}}\right)^{2}\,,$ (17b)

See how eq. 17a captures the super-linear scaling of volume with work that we
see in the data, Figures 6(d) and 6(b), and how the estimators can provide
quantitative predictions if the unknowns $\mathrm{k}$ and $G_{Ic}$ are chosen
accordingly, Figures 6(a) and 6(c).

In conclusion, a super-linear scaling of debris particle volume with
tangential work is observed at the nanoscale as the size of the asperities and
the junction patch increases over $d^{*}$.

(a) (a) prediction for $\theta=0^{\text{o}}$

(b) (b) scaling for $\theta=0^{\text{o}}$

(c) (c) prediction for $\theta=30^{\text{o}}$

(d) (d) scaling for $\theta=30^{\text{o}}$

Figure 6:  Verifying the new scalings using 3D data. Panel (a): debris volume
measured in silico (y-axis) v. estimator (x-axis), 3D simulations with
$\theta=0^{\text{o}}$. The value of $\mathrm{k}$ is estimated from the crack
pattern seen in the simulations, leaving $G_{Ic}$ as the only fit parameter;
we found $\mathrm{k}=0.285$ and
$G_{Ic}=G_{Ic}=29.0\,\mathrm{J/m^{2}}=181\,\mathrm{eV/nm^{2}}$. Panel (b):
data from 3D simulations with $\theta=0^{\text{o}}$ in logarithmic scale
alongside trend line for linear scaling ($\sim W_{t}$) and super-linear one
($\sim W_{t}^{1.5}$). Notice the intermediate transition zone from one to the
other. Panel (c): similar to (a) for 3D simulations with
$\theta=30^{\text{o}}$. In this case
$G_{Ic}=61.2\,\mathrm{J/m^{2}}=382\,\mathrm{eV/nm^{2}}$, the difference with
respect to the other value being probably due to the different geometry-
dependent crack pattern that is not accounted for. Panel (d): similar to (b)
but for $\theta=30^{\text{o}}$.

##  5 Discussion

Figure 6(d) strongly substantiates the scalings presented in eq. 16(a). On the
other hand, Figure 6(b) reveals a transition from a linear exponent to the
super-linear one. It must be highlighted how the geometry (e.g., the angle
$\theta$) affects the extent of the transition region. The super-linear
scaling is evident for all simulations where $\theta=30^{\text{o}}$, including
those closer to $d^{*}\approx 18\,\mathrm{nm}$, whereas the transition does
not happen until $d\approx 2d^{*}$ in the case $\theta=0^{\text{o}}$. Hence,
these results also give credence to the idea that the transition may depend on
parameters as, e.g., the geometry (as this one controls the stress
distribution).

We have been chiefly concerned with the qualitative trends in the tangential
work – debris volume relation, but we have also proposed new estimators eqs.
17a and 17b. However, it remains to assess the fracture toughness $G_{Ic}$ via
independent MD simulations to fully gauge their ability to predict numerical
outcomes, instead of considering it a fit parameter. There is still no
consensus as to how to estimate this material parameter using molecular
dynamics but a number of options are currently being investigated (Stepanova
and Bronnikov, 2020, Tong and Li, 2020, Patil et al., 2016).

We have shown in eq. 8 that a linear relation between tangential work and the
volume of worn material ($V_{d}\sim W_{t}$) can only hold when interlocking
asperities fully undergo inelastic deformation, which in turn is only the case
for the smallest asperities that lead to debris creation ($d>d^{*}$ while
$d\approx d^{*}$). Otherwise, we found $V_{d}\sim W_{t}^{2}$ in 2D (see
appendix) and $V_{d}\sim W_{t}^{1.5}$ in 3D. However, linearity between wear
volume and work on the macroscopic scale has been reported in many
experimental studies (Rabinowicz, 1995, Burwell and Strang, 1952, Whittaker,
1947, Uetz and Föhl, 1978, Fouvry et al., 2003). For example, Archard’s wear
model (Archard, 1953) states that

$\displaystyle V_{d}=\mathrm{k}{F_{n}\cdot S\over\mathrm{H}}\,,$ (18)

with $\mathrm{k}$ being the wear coefficient, $F_{n}$ the normal load, $S$ the
sliding distance, and $\mathrm{H}$ the material hardness. So if $F_{n}$ can be
related linearly to the tangential force $F_{t}$ (via any linear macroscopic
friction law), then it follows that $V_{d}\sim F_{t}S\sim W_{t}$.

A way to reconcile these facts is by accounting for the actual, statistical
process of wear particle formation during sliding on a rough surface
containing many asperity–asperity contacts. Contact simulations of rough
surfaces combined with the application of eq. 1 in a prior work (Brink et al.,
2021) have shown that most of the debris particles seem to arise from contacts
close to $d^{*}$. Because the microcontacts grow from an initially small size,
it stands to reason that their growth would be arrested by the emission of a
wear particle when reaching $d=d^{*}$. In such a case it would follow that
most of the worn mass stems from particles with sizes where the linearity is
recovered. Note that an initially narrow distribution of particle sizes can
later agglomerate into larger debris particles, giving rise to the plethora of
sizes that can be seen in experiments (see, for instance, Pham-Ba and Molinari
(2021) and Leriche et al. (2022)). It should be noted, though, that that model
assumes that wear particles are formed at isolated contact spots. It has been
proposed, however, that multiple, closely-spaced contacts can form a combined
wear particle in an even more efficient process, although this is likely only
the case under high normal load (Aghababaei et al., 2018, Pham-Ba et al.,
2020).

Even if the initial wear particle sizes were not restricted to the regime
where linearity holds, it could still be the case that this non-linear micro-
behavior gives rise to a linear relation as we upscale and more and more
contacts are considered simultaneously (from asperities to clusters, from
clusters to whole surfaces). Such a phenomenon, i.e., non-linear interactions
at the microscale resulting in a linear relation at the macroscale, is not
unheard of in the field of tribology. For instance it is well known since
Archard (Archard, 1957) and Greenwood and Williamson (Greenwood and
Williamson, 1966) that the sum of Hertzian contacts with a random distribution
of heights of contacting spheres—while strictly non-linear at the asperity
level (that is the circular contact area is a non linear function of the local
normal load)—gives rise to a linear dependence between real contact area and
macroscopic normal load. This suggests that future research should focus on
the collective behavior of asperity–asperity contacts on representative rough
surfaces.

An additional, important aspect is the effect of loading rate on the FPZ under
dynamic fracture conditions. It is well-known that the size of the FPZ depends
on the rate, so the effect of loading rate over wear should also be treated in
the framework presented here, for instance by adding correction terms. A
preliminary study on the influence of this parameter is carried out in the
appendix using idealized 2D simulations.

Finally, the possibility of other regimes where the volume–work scaling
relation changes is not ruled out. We have ascertained that the fracture
arrests before carving out the new particle completely, and that plastic
hinges form in the ligaments left between the arrested crack tip and the free
surface to finish the debris formation process (regard the strain distribution
in Figure 2 and the process depicted in Figure A.3). If the work that goes
into plastifying this region scales proportionally to the volume of the
asperities, then it may overcome the work associated to crack growth as the
main contribution to the overall tangential work. Note that this process is
geometry dependent, since the formation of the hinges occurs simultaneously
with the crack propagation in the 3D simulations, suggesting that the scaling
is not strongly affected in this specific case, at least. Techniques based on
the phase field method (Collet et al., 2020, Brach and Collet, 2021, Carollo
et al., 2019) seem ideal to further investigate such effects.

##  6 Conclusion

We have analyzed three-dimensional molecular dynamics simulations, at a single
asperity contact, resulting in the formation of adhesive wear particles. The
simulations revealed a super-linear scaling of debris size with tangential
work, contradicting a previous theoretical estimate with a linear scaling
(Aghababaei et al., 2017). A super-linear scaling could be observed because
the simulations cover a wide range of contact junction sizes, reaching clearly
above the critical junction size for debris formation (Aghababaei et al.,
2016). Our simulations indicate that for large contact junctions the process
of debris creation is more efficient than previously thought, that is, less
energy is required to form debris particles.

This has motivated the development of a new theoretical model, inspired by
ductile fracture mechanics, as the process of debris formation is always
accompanied by the propagation of cracks. The model accounts for two sources
of energy dissipation, the new surfaces formed along the crack path area and a
plastic volume around this crack path. The theory was shown able to explain
the trends observed in simulations.

It is important to highlight that in the limit of contact junction size
approaching the critical junction size, then the linearity between debris
volume and tangential work is fully recovered. The rationale is that the
critical junction size and the fracture process zone are essentially the same
length scale. This results in a total plastification of the contact junction.

Our findings provide a roadmap towards a quantitative framework to relate wear
debris volume and frictional work and should be further informed with
experimental observations.  The next task would be to devise independent MD
simulations to assess the fracture toughness. A particular point of interest
is that irrespectively of the asperity level mechanisms, a global linearity
between total wear volume and frictional work is generally observed. This
disconnection between asperity level mechanisms and global response requires
further studies and would certainly be enriched by considering the sliding
history.

## Acknowledgements

J. G.-S. and J.-F. M. gratefully acknowledge the support of the Swiss National
Science Foundation (grant 200021_197152, “Wear across scales”). J. G.-S. would
like to thank S. Z. Wattel for help with LAMMPS.

## Supplementary material

A Mathematica notebook (Wolfram, 2000) containing the computations leading to
results (including figures) shown in the text is provided as Supplementary
Material, and it can also be downloaded from the repository named ductile_wear
in the first author’s GitHub page github.com/jgarciasuarez. All other
materials necessary to reproduce results in this text can be obtained by
correspondence to the authors.

## References

  * Aghababaei (2019) R. Aghababaei.  On the origins of third-body particle formation during adhesive wear.  _Wear_ , 426-427:1076–1081, 2019.  doi: https://doi.org/10.1016/j.wear.2018.12.060.  22nd International Conference on Wear of Materials.
  * Aghababaei and Budzik (2020) R. Aghababaei and M. K. Budzik.  Fracture modes of brittle junctions under shear.  _Extreme Mechanics Letters_ , 35:100644, 2020.  doi: https://doi.org/10.1016/j.eml.2020.100644.
  * Aghababaei et al. (2016) R. Aghababaei, D. H. Warner, and J.-F. Molinari.  Critical length scale controls adhesive wear mechanisms.  _Nature Communications_ , 7(1):11816, Sept. 2016\.  doi: 10.1038/ncomms11816.
  * Aghababaei et al. (2017) R. Aghababaei, D. H. Warner, and J.-F. Molinari.  On the debris-level origins of adhesive wear.  _Proceedings of the National Academy of Sciences_ , 114(30):7935–7940, July 2017.  doi: 10.1073/pnas.1700904114.
  * Aghababaei et al. (2018) R. Aghababaei, T. Brink, and J.-F. Molinari.  Asperity-level origins of transition from mild to severe wear.  _Physical Review Letters_ , 120:186105, May 2018.  doi: 10.1103/PhysRevLett.120.186105.
  * Archard (1953) J. F. Archard.  Contact and rubbing of flat surfaces.  _Journal of Applied Physics_ , 24(8):981–988, 1953.  doi: 10.1063/1.1721448.
  * Archard (1957) J. F. Archard.  Elastic deformation and the laws of friction.  _Proceedings of the royal society of London. Series A. Mathematical and physical sciences_ , 243(1233):190–205, 1957\.
  * Baensch-Baltruschat et al. (2020) B. Baensch-Baltruschat, B. Kocher, F. Stock, and G. Reifferscheid.  Tyre and road wear particles (trwp)-a review of generation, properties, emissions, human health risk, ecotoxicity, and fate in the environment.  _Science of the Total Environment_ , 733:137823, 2020.
  * Barras et al. (2021) F. Barras, R. Aghababaei, and J.-F. Molinari.  Onset of sliding across scales: How the contact topography impacts frictional strength.  _Physical Review Materials_ , 5:023605, 2021.  doi: 10.1103/PhysRevMaterials.5.023605.
  * Bhaskaran et al. (2010) H. Bhaskaran, B. Gotsmann, A. Sebastian, U. Drechsler, M. A. Lantz, M. Despont, P. Jaroenapibal, R. W. Carpick, Y. Chen, and K. Sridharan.  Ultralow nanoscale wear through atom-by-atom attrition in silicon-containing diamond-like carbon.  _Nature Nanotechnology_ , 5(3):181–185, Mar. 2010\.  doi: 10.1038/nnano.2010.3.
  * Bhushan and Sundararajan (1998) B. Bhushan and S. Sundararajan.  Micro/nanoscale friction and wear mechanisms of thin films using atomic force and friction force microscopy.  _Acta Materialia_ , 46(11):3793–3804, July 1998\.  doi: 10.1016/S1359-6454(98)00062-7.
  * Bowden and Tabor (2001) F. P. Bowden and D. Tabor.  _The friction and lubrication of solids_ , volume 1.  Oxford university press, 2001.
  * Bowden et al. (1939) F. P. Bowden, D. Tabor, and G. I. Taylor.  The area of contact between stationary and moving surfaces.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 169(938):391–413, 1939\.  doi: 10.1098/rspa.1939.0005.
  * Brach and Collet (2021) S. Brach and S. Collet.  Criterion for critical junctions in elastic-plastic adhesive wear.  _Physical Review Letters_ , 127:185501, Oct 2021.  doi: 10.1103/PhysRevLett.127.185501.
  * Brink and Molinari (2019) T. Brink and J.-F. Molinari.  Adhesive wear mechanisms in the presence of weak interfaces: Insights from an amorphous model system.  _Physical Review Materials_ , 3:053604, May 2019.  doi: 10.1103/PhysRevMaterials.3.053604.
  * Brink et al. (2021) T. Brink, L. Frérot, and J.-F. Molinari.  A parameter-free mechanistic model of the adhesive wear process of rough surfaces in sliding contact.  _Journal of the Mechanics and Physics of Solids_ , 147:104238, 2021.  doi: https://doi.org/10.1016/j.jmps.2020.104238.
  * Brink et al. (2022) T. Brink, E. Milanese, and J.-F. Molinari.  Effect of wear particles and roughness on nanoscale friction.  _Physical Review Materials_ , 6(1):013606, Jan. 2022.  doi: 10.1103/PhysRevMaterials.6.013606.  Publisher: American Physical Society.
  * Brockley and Fleming (1965) C. A. Brockley and G. K. Fleming.  A model junction study of severe metallic wear.  _Wear_ , 8(5):374–380, Sept. 1965.  doi: 10.1016/0043-1648(65)90168-7.
  * Burwell and Strang (1952) J. T. Burwell and C. D. Strang.  On the empirical law of adhesive wear.  _Journal of Applied Physics_ , 23(1):18–28, 1952\.  doi: 10.1063/1.1701970.
  * Carollo et al. (2019) V. Carollo, M. Paggi, and J. Reinoso.  The steady-state archard adhesive wear problem revisited based on the phase field approach to fracture.  _International Journal of Fracture_ , 215(1):39–48, 2019.
  * Chung and Kim (2003) K.-H. Chung and D.-E. Kim.  Fundamental Investigation of Micro Wear Rate Using an Atomic Force Microscope.  _Tribology Letters_ , 15(2):135–144, Aug. 2003\.  doi: 10.1023/A:1024457132574.
  * Collet et al. (2020) S. Collet, J.-F. Molinari, and S. Brach.  Variational phase-field continuum model uncovers adhesive wear mechanisms in asperity junctions.  _Journal of the Mechanics and Physics of Solids_ , 145:104130, 2020.  doi: https://doi.org/10.1016/j.jmps.2020.104130.
  * Collins (1993) J. A. Collins.  _Failure of Materials in Mechanical Design: Analysis, Prediction, Prevention_.  Wiley-Interscience publication. Wiley, 1993.
  * Dieterich and Kilgore (1994) J. H. Dieterich and B. D. Kilgore.  Direct observation of frictional contacts: New insights for state-dependent properties.  _Pure and Applied Geophysics_ , 143(1):283–302, 1994.
  * Fouvry et al. (2001) S. Fouvry, P. Kapsa, and L. Vincent.  An elastic–plastic shakedown analysis of fretting wear.  _Wear_ , 247(1):41–54, 2001.  doi: https://doi.org/10.1016/S0043-1648(00)00508-1.
  * Fouvry et al. (2003) S. Fouvry, T. Liskiewicz, P. Kapsa, S. Hannel, and E. Sauger.  An energy description of wear mechanisms and its applications to oscillating sliding contacts.  _Wear_ , 255(1):287–298, 2003.  doi: https://doi.org/10.1016/S0043-1648(03)00117-0.
  * Frérot et al. (2018) L. Frérot, R. Aghababaei, and J.-F. Molinari.  A mechanistic understanding of the wear coefficient: From single to multiple asperities contact.  _J. Mech. Phys. Solids_ , 114:172–184, 2018.  doi: 10.1016/j.jmps.2018.02.015.
  * Fusco et al. (2010) C. Fusco, T. Albaret, and A. Tanguy.  Role of local order in the small-scale plasticity of model amorphous materials.  _Physical Review E_ , 82:066116, 2010.  doi: 10.1103/PhysRevE.82.066116.
  * Gotsmann and Lantz (2008) B. Gotsmann and M. A. Lantz.  Atomistic Wear in a Single Asperity Sliding Contact.  _Physical Review Letters_ , 101(12):125501, Sept. 2008.  doi: 10.1103/PhysRevLett.101.125501.  Publisher: American Physical Society.
  * Greenwood and Tabor (1955) J. A. Greenwood and D. Tabor.  Deformation Properties of Friction Junctions.  _Proceedings of the Physical Society. Section B_ , 68(9):609–619, Sept. 1955.  doi: 10.1088/0370-1301/68/9/305.  Publisher: IOP Publishing.
  * Greenwood and Williamson (1966) J. A. Greenwood and J. B. P. Williamson.  Contact of nominally flat surfaces.  _Proc. R. Soc. Lond. A_ , 295:300–319, 1966.  doi: 10.1098/rspa.1966.0242.
  * Griffith and Taylor (1921) A. A. Griffith and G. I. Taylor.  Vi. the phenomena of rupture and flow in solids.  _Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character_ , 221(582-593):163–198, 1921.  doi: 10.1098/rsta.1921.0006.
  * Grigoratos and Martini (2015) T. Grigoratos and G. Martini.  Brake wear particle emissions: a review.  _Environmental Science and Pollution Research_ , 22(4):2491–2504, 2015.
  * Holland and Marder (1998a) D. Holland and M. Marder.  Ideal brittle fracture of silicon studied with molecular dynamics.  _Physical Review Letters_ , 80:746–749, 1998a.  doi: 10.1103/PhysRevLett.80.746.
  * Holland and Marder (1998b) D. Holland and M. Marder.  Erratum: Ideal brittle fracture of silicon studied with molecular dynamics [Phys. Rev. Lett. 80, 746 (1998)].  _Physical Review Letters_ , 81:4029–4029, 1998b.  doi: 10.1103/PhysRevLett.81.4029.
  * Holm (2013) R. Holm.  _Electric contacts: theory and application_.  Springer Science & Business Media, 2013.
  * Hutchinson (1983) J. W. Hutchinson.  Fundamentals of the Phenomenological Theory of Nonlinear Fracture Mechanics.  _Journal of Applied Mechanics_ , 50(4b):1042–1051, 12 1983.  doi: 10.1115/1.3167187.
  * Jacobs and Carpick (2013) T. D. B. Jacobs and R. W. Carpick.  Nanoscale wear as a stress-assisted chemical reaction.  _Nature Nanotechnology_ , 8(2):108–112, Feb. 2013\.  doi: 10.1038/nnano.2012.255.
  * Janssen et al. (2004) M. Janssen, J. Zuidema, and R. Wanhill.  _Fracture Mechanics: Fundamentals and Applications_.  CRC Press, 2004.
  * Kole et al. (2017) P. J. Kole, A. J. Löhr, F. G. Van Belleghem, and A. M. J. Ragas.  Wear and tear of tyres: a stealthy source of microplastics in the environment.  _International journal of environmental research and public health_ , 14(10):1265, 2017.
  * Leriche et al. (2022) C. Leriche, S. Franklin, and B. Weber.  Measuring multi-asperity wear with nanoscale precision.  _Wear_ , 498-499:204284, 2022.  doi: https://doi.org/10.1016/j.wear.2022.204284.
  * Liu et al. (2010) J. Liu, J. K. Notbohm, R. W. Carpick, and K. T. Turner.  Method for Characterizing Nanoscale Wear of Atomic Force Microscope Tips.  _ACS Nano_ , 4(7):3763–3772, July 2010.  doi: 10.1021/nn100246g.  Publisher: American Chemical Society.
  * Meng and Ludema (1995) H. Meng and K. C. Ludema.  Wear models and predictive equations: their form and content.  _Wear_ , 181-183:443–457, 1995.  doi: https://doi.org/10.1016/0043-1648(95)90158-2.  10th International Conference on Wear of Materials.
  * Milanese et al. (2019) E. Milanese, T. Brink, R. Aghababaei, and J.-F. Molinari.  Emergence of self-affine surfaces during adhesive wear.  _Nature communications_ , 10(1):1–9, 2019.
  * Milanese et al. (2020) E. Milanese, T. Brink, R. Aghababaei, and J.-F. Molinari.  Role of interfacial adhesion on minimum wear particle size and roughness evolution.  _Physical Review E_ , 102:043001, 2020.  doi: 10.1103/PhysRevE.102.043001.
  * Morse (1929) P. M. Morse.  Diatomic molecules according to the wave mechanics. ii. vibrational levels.  _Physical Review_ , 34:57–64, Jul 1929.  doi: 10.1103/PhysRev.34.57.
  * Patil et al. (2016) S. P. Patil, Y. Heider, C. A. Hernandez Padilla, E. R. Cruz-Chú, and B. Markert.  A comparative molecular dynamics-phase-field modeling approach to brittle fracture.  _Computer Methods in Applied Mechanics and Engineering_ , 312:117–129, 2016.  doi: https://doi.org/10.1016/j.cma.2016.04.005.  Phase Field Approaches to Fracture.
  * Pham-Ba and Molinari (2021) S. Pham-Ba and J.-F. Molinari.  Creation and evolution of roughness on silica under unlubricated wear.  _Wear_ , 472-473:203648, 2021.  doi: https://doi.org/10.1016/j.wear.2021.203648.
  * Pham-Ba et al. (2020) S. Pham-Ba, T. Brink, and J.-F. Molinari.  Adhesive wear and interaction of tangentially loaded micro-contacts.  _International Journal of Solids and Structures_ , 188-189:261–268, 2020.  doi: https://doi.org/10.1016/j.ijsolstr.2019.10.023.
  * Pineau and Pardoen (2007) A. Pineau and T. Pardoen.  2.06 - failure of metals.  In I. Milne, R. O. Ritchie, and B. Karihaloo, editors, _Comprehensive Structural Integrity_ , pages 684–797. Pergamon, Oxford, 2007\.  doi: https://doi.org/10.1016/B0-08-043749-4/02109-1.
  * Plimpton (1995) S. Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  _Journal of Computational Physics_ , 117(1):1–19, 1995.  doi: https://doi.org/10.1006/jcph.1995.1039.
  * Rabinowicz (1958) E. Rabinowicz.  The effect of size on the looseness of wear fragments.  _Wear_ , 2(1):4–8, 1958.  doi: https://doi.org/10.1016/0043-1648(58)90335-1.
  * Rabinowicz (1995) E. Rabinowicz.  _Friction and wear of materials_.  Wiley, New York, 1995.
  * Rabinowicz et al. (1951) E. Rabinowicz, D. Tabor, and F. P. Bowden.  Metallic transfer between sliding metals: an autoradiographic study.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 208(1095):455–475, 1951\.  doi: 10.1098/rspa.1951.0174.
  * Renard et al. (2013) F. Renard, T. Candela, and E. Bouchaud.  Constant dimensionality of fault roughness from the scale of micro-fractures to the scale of continents.  _Geophysical Research Letters_ , 40(1):83–87, 2013.  doi: https://doi.org/10.1029/2012GL054143.
  * Renouf et al. (2011) M. Renouf, F. Massi, N. Fillot, and A. Saulot.  Numerical tribology of a dry contact.  _Tribology International_ , 44(7):834–844, 2011\.  doi: https://doi.org/10.1016/j.triboint.2011.02.008.
  * Reye (1860) T. Reye.  Zur Theorie der Zapfenreibung.  _Der Civilingenieur_ , 4(1860):235–255, 1860\.
  * Sato et al. (2012) T. Sato, T. Ishida, L. Jalabert, and H. Fujita.  Real-time transmission electron microscope observation of nanofriction at a single Ag asperity.  _Nanotechnology_ , 23(50):505701, Nov. 2012.  doi: 10.1088/0957-4484/23/50/505701.  Publisher: IOP Publishing.
  * Sorensen et al. (1996) M. R. Sorensen, K. W. Jacobsen, and P. Stoltze.  Simulations of atomic-scale sliding friction.  _Physical Review B_ , 53(4):2101–2113, Jan. 1996\.  doi: 10.1103/PhysRevB.53.2101.  Publisher: American Physical Society.
  * Stepanova and Bronnikov (2020) L. Stepanova and S. Bronnikov.  A computational study of the mixed–mode crack behavior by molecular dynamics method and the multi – parameter crack field description of classical fracture mechanics.  _Theoretical and Applied Fracture Mechanics_ , 109:102691, 2020.  doi: https://doi.org/10.1016/j.tafmec.2020.102691.
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber.  Computer simulation of local order in condensed phases of silicon.  _Physical Review B_ , 31:5262–5271, Apr 1985.  doi: 10.1103/PhysRevB.31.5262.
  * Stoyanov et al. (2014) P. Stoyanov, P. A. Romero, R. Merz, M. Kopnarski, M. Stricker, P. Stemmer, M. Dienwiebel, and M. Moseler.  Nanoscale sliding friction phenomena at the interface of diamond-like carbon and tungsten.  _Acta Materialia_ , 67:395–408, Apr. 2014.  doi: 10.1016/j.actamat.2013.12.029.
  * Stukowski (2009) A. Stukowski.  Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool.  _Modelling and Simulation in Materials Science and Engineering_ , 18(1):015012, dec 2009.  doi: 10.1088/0965-0393/18/1/015012.
  * Thompson et al. (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  _Comput. Phys. Commun._ , 271:108171, Feb. 2022.  ISSN 0010-4655.  doi: 10.1016/j.cpc.2021.108171.  https://lammps.org/.
  * Timoshenko (1922) S. Timoshenko.  On the transverse vibrations of bars of uniform cross-section.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 43(253):125–131, 1922.  doi: 10.1080/14786442208633855.
  * Tong and Li (2020) Q. Tong and S. Li.  A concurrent multiscale study of dynamic fracture.  _Computer Methods in Applied Mechanics and Engineering_ , 366:113075, 2020.  doi: https://doi.org/10.1016/j.cma.2020.113075.
  * Uetz and Föhl (1978) H. Uetz and J. Föhl.  Wear as an energy transformation process.  _Wear_ , 49(2):253–264, 1978.  doi: https://doi.org/10.1016/0043-1648(78)90091-1.
  * Vakis et al. (2018) A. Vakis, V. Yastrebov, J. Scheibert, L. Nicola, D. Dini, C. Minfray, A. Almqvist, M. Paggi, S. Lee, G. Limbert, J. Molinari, G. Anciaux, R. Aghababaei, S. Echeverri Restrepo, A. Papangelo, A. Cammarata, P. Nicolini, C. Putignano, G. Carbone, S. Stupkiewicz, J. Lengiewicz, G. Costagliola, F. Bosia, R. Guarino, N. Pugno, M. Müser, and M. Ciavarella.  Modeling and simulation in tribology across scales: An overview.  _Tribology International_ , 125:169–199, 2018.  doi: https://doi.org/10.1016/j.triboint.2018.02.005.
  * Whittaker (1947) E. J. W. Whittaker.  Friction and wear.  _Nature_ , 159(4042):541–541, 1947.
  * Wolfram (2000) S. Wolfram.  _The mathematica book_ , volume 4.  Cambridge University Press Cambridge, 2000.
  * Zhong et al. (2013) J. Zhong, R. Shakiba, and J. B. Adams.  Molecular dynamics simulation of severe adhesive wear on a rough aluminum substrate.  _Journal of Physics D: Applied Physics_ , 46(5):055307, Jan. 2013.  doi: 10.1088/0022-3727/46/5/055307.  Publisher: IOP Publishing.

##  Appendix A Supplementary material: 2D simulations

###  A.1 Numerical setting and results

(a) (a)

(b) (b)

Figure A.1: Evidencing the limitations of the debris volume estimate eq. 2
(Aghababaei et al., 2017) using both 2D simulations. Panel (a): scheme of the
2D simulations, including junction size and loading conditions. The geometry
of the system is parametrized in terms of $d$, in particular $R=3d/4$, thus
the gap between surfaces is $3/2d$.  Panel (b): comparison of debris volume
measured in silico (black) v. estimate (orange) in 2D simulations ,
$v_{top}=0.02r_{0}/t_{0}$. The shaded orange region corresponds to maximum and
minimum values of the estimates presented in (Aghababaei et al., 2016). The
shaded blue region correspond to maximum and minimum values of the eq. 17b.

Due to the large size of the simulations (on the order of 12 million atoms),
we were, however, not able to run many different geometries and could not run
statistically independent repetitions for a given asperity size. Thus, we also
used the computationally more affordable 2D simulations. Here, we could run a
number of realizations of each case by initializing the atoms with different
random velocity distributions consistent with the desired temperature.

These simulations hardly represent a practical situation (neither the
simplified potential nor the 2D configuration are realistic), but they are not
exempt of academic interest: in first place, they allow us to test the 2D
estimator, eq. 17b and, moreover, the geometry in this case, fig. A.1(a),
features no sharp edges, so the lack of stress concentration leads to more
variety in crack paths, what allow us to test the scaling in a disadvantageous
situation. Finally, in this suite of virtual experiments we were able to
repeat situations varying only the loading velocity $v_{top}$, so a first
assessment of the result sensitivity to rate effects can be carried out.

For each size, six realizations were run and the results are reported in
Figure A.1(b) in black with error bars representing the range of values. The
volume estimator (eq. 2) is shown in orange. It also covers a range of values
indicated by the shaded region due to fluctuations of the total tangential
work in the different realizations.

We used a modified Morse potential of the form (Morse, 1929, Aghababaei et
al., 2016)

$\displaystyle{V(r)\over\epsilon}=\begin{cases}(1-\exp{\\{-\alpha(r-r_{0})\\}})^{2}-1&\text{\,
if \,}0<r\leq 1.1r_{0}\\\ {c_{1}\over 6}r^{3}+{c_{2}\over 4}r^{2}+{c_{3}\over
6}r+c_{4}&\text{\, if \,}1.1r_{0}<r\leq r_{cut}\\\ 0&\text{\, if
\,}r>r_{cut}\end{cases}\,.$ (19)

We used $\alpha=3.93\,r_{0}^{-1}$ and expressed all data in reduced units of
$\varepsilon$ and $r_{0}$. The cutoff distance was chosen as
$r_{cut}=1.422\,r_{0}$ and the parameters $c_{1},\,c_{2},\,c_{3},\,c_{4}$ were
selected to ensure continuous energy and force at $r=1.1\,r_{0}$, as well as
zero energy and force at the cutoff. This potential was used in Aghababaei et
al. (2016) and labeled “P6”. The timestep is chosen to be $\Delta
t=0.005t_{0}=0.005\,r_{0}/\sqrt{\varepsilon/m}$ (where $m$ is the atomic
mass).

We then used atoms on a hexagonal lattice and cut out the geometry sketched in
Fig. 1(c). The height of the interlocking asperity is $2R=3d/2$. The box size
is chosen proportionally to $d$: its total height $H$ equals $8d$, while its
horizontal length is $L=5d$. These values are chosen to guarantee that the
model’s horizontal edges and the periodic boundary conditions (PBCs) on the
vertical ones do not affect the local stress state at the junction. Langevin
thermostats are located on top and bottom layers to enforce a constant
temperature $T=0.025\,\epsilon/k_{B}$, where $k_{B}$ represents the Boltzmann
constant. The bottom layer is constrained to remain fixed, while a normal
pressure of $0.03\,\epsilon/r_{0}^{3}$ and a sliding velocity are applied to
the top one. These simulations are run over a minimum sliding distance equal
to $3d$, to ensure enough sliding so as to trigger debris creation. The
sliding velocity is $0.02\,r_{0}/t_{0}$, additional simulations that were run
at velocity $0.05\,r_{0}/t_{0}$ are also reported. The pre-formed junction
sizes considered in these simulations range from $d=30\,r_{0}$ to
$d=100\,r_{0}$ in increments of $10\,r_{0}$. The minimal value of $d$ that
yields debris creation is $d^{*}\approx 15\,r_{0}$, meaning all these
simulations lead to third-body formation. The number of atoms in the largest
box is $\approx 40600$ which makes these 2D simulations relatively affordable,
thereby allowing six independent realizations of each size for statistics.

We estimate the debris volume (area in 2D), by counting the number of atoms in
the debris particle and multiplying them by the average atomic volume (area)
in the bulk.

We also estimate the shear strength via independent simulations, finding
$\tau=0.78\,\epsilon/r_{0}^{3}$.

Figure A.2: Verifying the assumptions behind eq. 2. (a) Assumption #1: maximum
force. Error bars represent the range of values obtained in numerical
simulations ($v_{top}=0.02r_{0}/t_{0}$), dashed line represents the values
from Continuum Mechanics, solid blue line represents the linear trend defined
by the last two points. $G_{Ic}$ in our case represents a fit parameter once
we fix an assumed total crack length ($l_{crack}=6d$, see Figure A.3), its
value $4.86\,\epsilon/r_{0}^{2}$ is inferred directly from simulations by
fitting the linear data presented in panel(c). (b) Assumption #2: effective
sliding distance. Error bars represent the range of values obtained in
numerical simulations ($v_{top}=0.02r_{0}/t_{0}$), the dashed line represents
the values from Continuum Mechanics. (c) Work scaling: linear increase in
tangential work necessary to generate debris, $v_{top}=0.02r_{0}/t_{0}$
(squares) and $v_{top}=0.02r_{0}/t_{0}$ (circles), notice apparent rate
effects (deviation from linear trend) for largest size and larger velocity.

Taking advantage of the fact that the height that appears in eq. 12 can be
written as $h=3/4d$, Figure 1(c), Figure A.2 panel (a) depicts the numerical
results for $S_{eff}/d$ from the 2D model in Figure 1(a) concurrently with a
scaling trend consistent with eq. 12. Note that the prefactor is adapted
slightly to account for the shape of the asperities in the simulation not
being prismatic.

Even though the 2D results we have obtained agree well with the predicted
trends, we remark that for the largest size we tested, $d=100\,r_{0}$, there
is a difference between the results for the two velocities considered,
Figure4(a). This hints to a growing sensitivity to loading rate as the
asperity size increases.

###  A.2 Example of fracture process

Figure A.3: Snapshots of the third-body formation process (for size
$d=90\,r_{0}$), images include cumulative dimensionless sliding $\Delta s/d$.
Process runs from left to right, from top to bottom.

##  Appendix B Further evidence of propagation regimes

Figure B.4: Influence of crack propagation during wear. (a) Scheme of the
displacement-controlled system with “handshake” pre-formed junction (as
presented in Pham-Ba et al. (2020)), all junctions $h=6r_{0}$. (b) Snapshots
of crack propagation leading to debris detachment. (c) Tangential force
evolution for different junction sizes, the blue triangles highlight the slope
of the unloading in the case $d/h=5.33$ (slow unloading proportional to
imposed velocity) and $d/h=10$ (fast unloading proportional to wave velocity).
(d) Tangential work evolution for different sizes.

This study uses quasi-2D simulations, with atoms arranged on an FCC lattice,
featuring rectangular pre-formed junctions (Pham-Ba et al., 2020) (the
trajectory of the asperities prior to contact was neglected). All units are
expressed in the corresponding potential reduced units, similarly to the 2D
simulations presented in the body of the article, see Section 2. These
correspond to displacement-controlled virtual experiments, where the top
surface sliding velocity is $0.01r_{0}/t_{0}$. Five different junction widths
$d$ were considered, $d/r_{0}=30,\,32,\,35,\,40,\,60$ ($r_{0}$ is the distance
between the atom center and the potential well), while fixing the junction
height $h/r_{0}=6$, thus $d/h=5,\,5.33,\,5.83,\,6.67,\,10$. All simulations
run until the total sliding distance reaches $50r_{0}$. The critical junction
size was $d^{*}\approx 35r_{0}$. All the other details can be checked in Pham-
Ba et al. (2020).

The evolution of the tangential force for each size is visualized in Figure
B.4(c). The qualitative dissimilarities are evident. In the lightest color,
the smallest size, $d/r_{0}=30$ and $d/h=5$, was chosen so as not to generate
debris (plastic asperity smoothing, since $d^{*}\approx 31r_{0}$). All the
others do lead to particle separation (notice the different convexity of the
work evolutions in Figure B.4(c)).

See that, as the junction width increases, gradients in the force profile
become acuter; for the two largest sizes, there is a sudden drop, which
entails drastic loss of stiffness and fast energy release, coherent with
unstable crack propagation. The following bump represents work invested not in
separation, but in rotating the already-formed particle out of the hollow left
behind by the atoms that now belong to the contour of the debris.

The qualitative change evidences the transition between the strength-
controlled regime and the toughness-controlled one.

##  Appendix C Snapshots from 3D simulations for $\theta=0^{\circ}$

Figure C.5: Plastic strain for the 3D simulations with $\theta=0^{\circ}$. A
slice through the middle of the wear particle is shown to visualize the
plasticity in the bulk. (a) For the smallest asperity diameter close to
$d^{*}=18\,\mathrm{nm}$, significant plasticity occurs in the bulk of the
particle. (b)–(d) With increasing asperity size, the plasticity becomes more
localized and does not fill the whole particle.