Interplay of mechanics and chemistry governs wear of diamond-like carbon coatings interacting with ZDDP-additivated lubricants

Friction and wear reduction by diamond-like carbon (DLC) in automotive applications can be affected by zinc-dialkyldithiophosphate (ZDDP), which is widely used in engine oils. Our experiments show that DLC’s tribological behaviour in ZDDP-additivated oils can be optimised by tailoring its stiffness, surface nano-topography and hydrogen content. An optimal combination of ultralow friction and negligible wear is achieved using hydrogen-free tetrahedral amorphous carbon (ta-C) with moderate hardness. Softer coatings exhibit similarly low wear and thin ZDDP-derived patchy tribofilms but higher friction. Conversely, harder ta-Cs undergo severe wear and sub-surface sulphur contamination. Contact-mechanics and quantum-chemical simulations reveal that shear combined with the high local contact pressure caused by the contact stiffness and average surface slope of hard ta-Cs favour ZDDP fragmentation and sulphur release. In absence of hydrogen, this is followed by local surface cold welding and sub-surface mechanical mixing of sulphur resulting in a decrease of yield stress and wear.

D iamond-like carbon is an excellent solid lubricant that exhibits high wear resistance and superlubricity (i.e. friction coefficient µ below 0.01) under boundary 1,2 and mixed lubrication 3 with organic friction modifiers (FM), including oleic acid and glycerol. It has been increasingly applied as a protective coating in internal combustion engines, where operating conditions are extremely severe and thus the use of antiwear (AW) additives, e.g. zinc dialkyldithiophosphate (ZDDP) and FM additives, e.g. molybdenum dithiocarbamate (MoDTC), admixed with base oils is essential for reducing wear and friction. However, some additives (e.g. MoDTC 4 ) can cause massive wear of some DLCs. Particularly the role of ZDDP in the wear of DLC/ DLC contacts remains controversial [5][6][7][8][9] since research on these additives has mainly focused on conventional ferrous surfaces 10 . Consequently, a better understanding of interactions between such additives and DLC as well as an improvement of the AW performance of additives on DLC are technologically important.
ZDDP is one of the most indispensable and effective AW additives for automotive combustion engines 10 . However, since sulphur and phosphorous produce environmentally hazardous waste substances and degrade catalyst converters, alternative additives have been sought to replace ZDDP with other ash-and phosphorous-free AW agents for the sake of reducing the environmental load 11,12 . Nevertheless, owing to its extraordinary AW properties, ZDDP is still a principal AW additive. Substantial efforts have been made to identify chemical structures and formation mechanisms of ZDDP-derived tribofilms [13][14][15][16] , which is necessary not only to improve the AW performance of ZDDP but also to formulate alternative sulphur-and phosphorous-free lubricants.
Particular attention is paid to whether ZDDP is effective with DLC, which has been already investigated for various types of DLC coatings 5,6,9,17,18 , including hydrogen-free tetrahedral amorphous carbon (ta-C) and hydrogenated amorphous carbon (a-C:H). Vengudusamy et al. carried out comprehensive experimental studies on the effect of different types of DLC coatings, such as a-C:H, a-C, ta-C and doped DLCs [5][6][7] . They reported that ta-C lowers boundary friction, but exhibits higher wear compared with a-C:H and a-C. They concluded that DLC's wear resistance is mainly governed by the type of DLC, but the role of ZDDP as AW agent is unclear.
Interestingly, ZDDP-derived patchy films can form on DLC surfaces, but are less thick and durable than those observed on steel, probably due to weaker interactions of ZDDP with DLC surfaces 7,19 . The thicknesses and morphologies of these tribofilms are affected by intrinsic parameters of the DLC coatings, i.e. hydrogen content, sp 2 /sp 3 ratio, presence of dopants, surface roughness and mechanical properties 6,20 . The former three factors alter chemical reactivity of DLC surfaces and thus reaction rates of ZDDP, whereas the latter two factors affect the mechanical response of the system such as contact pressure and geometry, and elastic/plastic deformation of surface asperities.
Mosey et al. 21 simulated the cross-linking of bulk zinc phosphate precursors with ab-initio molecular dynamics (MD), showing that high hydrostatic pressures induce the formation of poly-phosphates. It would be tempting to ascribe the low quantity of ZDDP-derived polyphosphates on DLC to reduced contact pressures. However, DLCs are stiffer than steel and therefore ZDDP on DLC should experience higher contact pressures. This indicates that the chemical interactions between ZDDP and surfaces cannot be neglected in ab-initio modelling. Yet, only decomposition of other sulphur- 22 and phosphorous-containing additives 23 on iron surfaces have been simulated so far.
Here, we combine experiments and simulations to understand tribochemical reactions of ZDDP with DLC and the role of DLC's mechanical properties, surface roughness and chemistry. We study friction and wear of five different self-mated DLC coatings sliding in a ZDDP solution with a PAO base oil. First, a combination of sliding tests, surface chemical analyses, and contactmechanics calculations reveal the crucial role of asperity-scale contact mechanics on wear and friction. Surface asperities on soft DLC experience lower effective contact pressures (P l;eff ) than asperities on hard DLC. While 20-nm-thick tribo-patches form on the former and no visible wear is observed, the latter undergo severe wear and a drastic friction increase. Surface chemical analyses reveal significantly more S−C than Zn−S bonds on hard DLCs, indicating that ZDDP reacts preferentially with hard DLC coatings. Second, quantum-chemical simulations reveal that the kinetics of ZDDP's chemical-decomposition reactions is accelerated by P l;eff while it is independent of the DLC's bulk chemical structure. However, chemical differences between DLC surfaces alter the shear response of the frictional interface at high P l;eff (>5 GPa). Non-hydrogenated DLC surfaces are likely to cold-weld at high P l;eff owing to lower steric hindrance and strong attractive interactions. The cold-welded surfaces undergo shear-induced chemical mixing, resulting in sub-surface sulphur doping, which in turn weakens the DLC matrix (reduces its yield stress) causing severe wear. In contrast, local contact pressures on soft DLC are small enough (P l;eff ≲3:5 GPa) to prevent cold welding and sulphur doping. In this case, ZDDP chemisorbs and decomposes into atoms and small fragments that accumulate and cross-link to form anti-wear tribofilms. This study elucidates how the interplay between inter-asperity nano-mechanics and chemistry governs tribochemistry and macroscopic tribological behaviours of ZDDP/DLC systems. The effective local pressure at asperity contacts and the coatings' hydrogen content are the crucial mechanical and chemical factors, respectively. The former controls ZDDP's mechanochemistry, whereas the latter affects shear accommodation at the sliding interface.

Results
Reciprocating sliding tests. In order to investigate the impact of ZDDP on friction and wear of DLC coatings, reciprocating cylinder-on-disc sliding tests are performed for self-mated DLC/ DLC tribopairs in both pure poly-alpha-olefin (PAO-4) base oil and PAO-4 mixed with 1 wt% ZDDP. Five DLC coatings with varying properties are tested under the same conditions (Table 1). They are labelled by a-C:H, a-C, ta-C(51), ta-C(66) and ta-C(78). The values in parenthesis represent the hardness of ta-Cs in GPa.
The tests were repeated at least 4 times for each DLC to ensure the repeatability of our experiments. Significant differences in both friction and wear are observed between pure PAO and PAO + ZDDP lubrication. Figure 1a, b show the evolution of friction in pure PAO and PAO + ZDDP, respectively. For a-C and a-C:H, µ increases gradually and reaches a steady state with μ % 0:1 in PAO (Fig. 1a). Conversely, in PAO + ZDDP, the running-in periods are much shorter than those in PAO and µ decreases slightly to 0.08−0.09 (Fig. 1b). The three ta-Cs reach rapidly a steady, ultralow friction state (μ ¼ 0:01 À 0:03) in PAO. In PAO + ZDDP, the friction curve for ta-C(51) is almost the same as that in PAO. However, for harder ta-Cs, i.e. ta-C(66) and ta-C (78), the addition of ZDDP significantly changes their tribological performance. Friction increases gradually to μ ¼ 0:08 À 0:10 after 2000 and 12000 sliding cycles for ta-C(66) and ta-C(78), respectively. Figure 1c shows composite wear volumes W measured on discs and cylinders. In PAO, except for ta-C(78), the wear volumes are negligibly small. However, in PAO + ZDDP wear increases drastically for ta-C(66) and ta-C(78) in correlation with the increased µ (Fig. 1b). Optical images of wear tracks on both cylinders and discs after sliding show deep scratches for ta-C(66) and ta-C(78) (Fig. 1d). The rubbed surfaces are much rougher than the pristine surfaces (interferometer images in Supplementary Fig. 2). The average wear scar depths are about 0.5 and 1.5 μm for the disc and cylinder, respectively. Conversely, the wear scars of the softer DLCs are barely visible under the same observation conditions. Thus, ZDDP has a negligible effect on friction and wear of soft/moderately-hard DLCs but a catastrophic effect on those of harder ta-Cs.
Contact mechanics calculations. Our friction experiments are carried out with the same initial macroscopic cylindrical Hertzian contact pressure (P Hertz ¼ 210 MPa) for all tribo-pairs. However, the DLC coatings differ in their mechanical properties and surface topographies (Table 1) and thus local contact pressures P l at the asperity level can be different. Here, we perform contact mechanics calculations using the numerical boundary element method with experimental surface topographies to estimate these effects at the nanoscale (details in Methods). Surface topographies of all five DLC coatings are measured by atomic force microscopy (AFM) with a scan area of 20 × 20 μm 2 inside the wear tracks of cylinders and discs after running-in. The AFM topographies of cylinder and disc are brought into contact under an external pressure of 210 MPa. For the softer DLCs (a-C:H, a-C, and ta-C (51)), the P l probability density distributions have a high peak around P l ¼ 0 and decay rapidly, with most asperities experiencing a P l lower than a few GPa (Fig. 1e). In contrast, the P l distributions for ta-C(66) and ta-C(78) are broader and contact pressures P l >10 GPa occur frequently. Figure 1f shows the effective local contact pressure P l;eff (as the median of the P l distribution over all contact spots where P l >0) as a function of a material parameter E * h 0 rms (the product of the contact modulus E * and the RMS slope h 0 rms , discussed later on). For ta-C(66) and ta-C(78), P l;eff is 9.2 and 10.5 GPa, respectively. In contrast, for the three softer DLCs, P l;eff values are in the range 1.7-3.5 GPa. Interestingly, P l;eff for ta-C(51) is smaller than P l;eff for a-C (which has a larger h 0 rms ). The monotonic increase in P l;eff with E * h 0 rms indicates that P l;eff is not only a function of the Young's modulus E but also depends on surface topography.
Surface chemical analyses. Surface chemical analyses are performed using X-ray photoelectron spectroscopy (XPS) inside and outside the wear tracks of the discs to study the mechanisms underlying the tribological behaviours of DLC surfaces under PAO + ZDDP lubrication. The binding energies (BEs) of the C1s, O1s, P2p, S2p and Zn2p 3/2 photo-peaks are analysed carefully (see Methods for details of the fitting procedure and XPS spectrometer). Figure 2a shows the relative chemical compositions estimated from XPS spectra in the topmost 10 nm of the DLC discs. The ZDDP-derived elements, namely P, S, and Zn, are detected inside the wear tracks of all DLC surfaces (Fig. 2a), whereas very small traces of these elements are found outside the wear tracks. The atomic concentrations of P, S, and Zn inside the wear tracks decrease in the order: a-C:H>a-C>ta-Cs.
While the chemical states of C, O, P, and Zn on the five DLC coatings ( Supplementary Fig. 3) are hardly distinguishable, this is not the case for sulphur. Figure 2b shows S2p XPS spectra inside and outside wear tracks on flat specimens. Outside the wear tracks, similar S2p XPS spectra are observed for all DLC surfaces with a peak at BE = 162.4 eV (Fig. 2b, right column), which is very close to that of sulphur in neat ZDDP (BE = 162.5 eV). The other peak at BE = 168.5 eV indicates the presence of sulphates (presumably due to thermal oxidation of DLC surfaces and subsequent chemisorption of ZDDP). In contrast, inside the wear track of a-C:H (Fig. 2b, left column), the S2p spectrum has one single component (in blue) at BE = 162.3 eV corresponding to S−Zn 7,24 , or S−P and S=P in thiophosphates 24 . This peak is present also inside the wear tracks of the other DLCs, but the spectra contain another component (in red) at BE = 163.3 eV due to S−C bonds 25 . Note that C−S−H, C−S −C, or C−S−S−C bonding states are practically indistinguishable as their differences in the BEs are only about 0.1 eV 25 . The other component for ta-Cs is attributed to a peak at BE = 167.9 eV due to sulphur oxides 26,27 . Interestingly, the contribution from the S−C component at BE = 163.3 eV is much smaller than the one at BE = 162.3 eV for the a-C and ta-C(51), and thus sulphur in these coatings is preferentially bound to zinc and phosphorous (in the form of S−P and S=P). In contrast, for ta-C(66) and ta-C(78), the BE = 163.3 eV contribution is dominant, indicating the breaking of S−Zn bonds in ZDDP and the formation of S−C bonds. These results indicate that ZDDP is likely to stay intact on the three softer  , a-C (green), ta-C(51) (orange), ta-C(66) (red) and ta-C(78) (brown). c Composite wear volumes measured on the discs and cylinders by optical interferometry and d optical images of the wear tracks of the DLC-coated cylinders and discs after sliding. e Local contact pressure distributions calculated from contact mechanics calculations using AFM topographies measured inside the wear track after running-in (2000 cycles) with a scan area of 20 × 20 µm 2 . f Effective local contact pressure P l;eff as a function of E Ã h 0 rms , where E * is the contact modulus and h 0 rms is the combined root-mean-square slope between two DLC surfaces. P l;eff is defined as the median over all nonzero P l values and is obtained from a numerical boundary element solution of two experimental DLC topographies under a normal pressure of 210 MPa. Filled circles (with the same colour coding as in a and b mark the material parameter E Ã h 0 rms of each DLC coating. The boxes represent the range of the P l distribution (limited by the 25th and 75th percentiles).
DLCs, whereas it undergoes decomposition on the hard ta-Cs. Moreover, for ta-C(66), the XPS spectra after 5-minute sputtering (corresponding to a sputter depth of about 0.5 nm 28 ) remain similar to those before sputtering inside the wear track ( Fig. 2b) and show the presence of significant amounts of S−C bonds (see Fig. 2c). In contrast, for a-C:H, no peaks attributed to S−C bonds are observed. This indicates that sulphur is present in the subsurface region of hard ta-Cs but not of a-C:H. Figure 2d shows a clear correlation between wear volume and S−C/S−Zn ratio.
Scanning and transmission electron microscopy characterisations. It is still not clear whether ZDDP-derived tribofilms form on DLC surfaces. Thicknesses, morphologies and chemical compositions of ZDDP-derived tribofilms are scrutinised using Scanning Electron Microscopy (SEM), Energy Dispersive X-Ray Spectroscopy (EDX) and Transmission Electron Microscopy (TEM). Figure 3 shows representative SEM images of the five DLC surfaces after sliding in PAO + ZDDP. ZDDP-derived tribofilms form on a-C:H, a-C, and ta-C(51) (Fig. 3a-c), as evidenced by superimposing the chemical map of zinc. Quantification results of EDX spectra are available in Supplementary Fig. 4 and Supplementary Table 2. On a-C:H, ZDDPderived tribopatches are heterogeneous in size, ranging from 0.1 to 1 μm (Fig. 3a). On a-C, tribopatches are more homogenous in size and distribution (Fig. 3b). On ta-C(51), tribopatches are found in the least amount ( Fig. 3c), in agreement with the XPS results. In all these cases the tribo-patches form predominantly at the edges of the stroke ( Supplementary Fig. 5). For ta-C(66) and ta-C(78), scratches are observed and no ZDDP tribopatches are visible (Fig. 3d, e). Instead, small ZDDP-derived particles (marked by red circles) are detected using punctual EDX spectra. The chemical structure and thickness of the ZDDP-derived tribofilms on a-C:H are identified by inspecting a cross-section of the a-C:H coating that exhibits the largest tribopatches after sliding using Focused Ion Beam preparation and TEM analysis. A bright-field (BF) TEM image reveals that the ZDDP-derived tribolayer (marked by a dashed green ellipse) is amorphous and discontinuous (Fig. 4a). Its thickness is about 20 nm, i.e. 5−10 times thinner than typical pad-like tribofilms on steel surfaces under similar sliding conditions 14 , and in agreement with a previous ball-on-disc rolling/sliding experimental result 7 .  Table 3). Note that Combined BF-TEM images of the cross-section of a ta-C(66) sample after sliding show the formation of large grooves due to wear and filling of a groove with smeared wear particles (indicated by arrows in Fig. 4d). A clear boundary between the smeared layer of wear debris and ta-C bulk is observed (Fig. 4e). Chemical information is obtained by EDX inside the smeared layer and ta-C bulk (for ROI 4−8 specified in Fig. 4f). While ROI 4−6 contain sulphur and less zinc, ROI 7 and 8 contain mostly carbon and traces of oxygen and the chemical compositions are typical of the native ta-C. This indicates that sulphur is transported into the ta-C subsurface region. Interestingly, a higher magnification image ( Fig. 4g) reveals that the sulphur-rich bright regions (ROI 4−6) undergo phase transformation and exhibit ordered graphitic structures, which is consistent with a previous study by Berman et al. showing sulphur diffusion into diamond nanoparticles and subsequent formation of an onion-like structure 29 .
In summary, our friction experiments combined with surface analyses and contact-mechanics calculations suggest that friction and wear of DLCs in PAO + ZDDP depend on asperity-scale contact pressures. For softer DLCs, the majority of asperities experience mild local contact pressures (less than a few GPa), which seems to be sufficient for the growth of a 20-nm-thick, patchy tribofilm but insufficient for the complete decomposition of ZDDP. In contrast, for harder DLCs, higher P l;eff (>9 GPa) might promote breaking of Zn−S bonds in ZDDP and preferential formation of S−C bonds, accompanied by penetration of sulphur into the sub-surface region and massive wear.
Quantum-chemical calculations. The experimental correlation between wear and the presence of S−C bonds (Figs. 2 and 4) can be explained by weakening of ta-C via incorporation of small amounts of sulphur (Fig. 5). We carry out MD homogeneous shearing simulations of sulphur-doped DLCs using a third-order density-functional tight-binding (DFTB3) method 30 . Lees-Edwards boundary conditions 31 are employed in this nonequilibrium MD simulations to impose a shear flow in a representative volume element of a bulk a-C, a-C:H and ta-C system (see "Methods" for samples preparation).
For all three systems, the average shear stress τ (the yield stress) decreases with increasing sulphur concentration C S (Fig. 5b). While τ is almost identical for a-C and a-C:H, it is larger for ta-C (as expected). For a-C:H and a-C (Fig. 5b, red and blue curves, respectively), the yield stresses at C S ¼ 0 at% are about 11-12 GPa, i.e. almost the same as in ta-C at C S ¼ 8 at.% (Fig. 5b, green curve). Thus, significant sulphur doping of the hardest ta-Cs can lower their yield stress making the resulting sulphur-carbon phase much weaker than pure a-C:H and a-C. But how does sulphur doping occur and how do reactions of ZDDP depend on coating properties such as stiffness, surface topography and chemical structure? Our experiments and contact mechanics calculations (Figs. 1-4) indicate that local contact pressures P l between surface asperities vary with the contact modulus E * and RMS slope h 0 rms of the coatings. Differences in both P l and surface reactivity of the DLC surfaces may trigger different mechanochemical reactions of ZDDP and thus influence the macroscopic tribological behaviour. Thus, additional simulations are performed in order to elucidate the dependence of ZDDP decomposition reactions on contact pressure and DLC chemical structure. DLC/DLC tribocontacts in presence of ZDDP are modelled using both quasi-static quantum-mechanical molecular statics (QMS) and QMD simulations (computational details in "Methods"). Two types of DLC surfaces are considered: (i) an a-C:H with bulk density ρ = 2.0 g cm −3 and hydrogen concentration C H = 20 at% and (ii) an a-C with the same density. The a-C sample represents the near-surface region of ta-C coatings, since ta-C is covered with a several-nm-thick a-C layer due to the nature of the deposition process 32 and a tribo-induced sp 3 -to-sp 2 rehybridisation 33 . Figure 6a displays an initial configuration of the QMS system, where a ZDDP molecule with ethyl groups is positioned between two a-C:H surfaces. Unpassivated surfaces are considered since asperity collisions can remove passivating species mechanically, leaving behind reactive carbon atoms 33 . A full hydrogen passivation would inhibit reactions of ZDDP even under extremely high contact pressures (P z >10 GPa) ( Supplementary  Fig. 6). As the two surfaces come gradually into contact, ZDDP reacts with the upper surface even at a low contact pressure (P z % 0:1 GPa in Fig. 6b). The reactions start with the breaking of a Zn−S bond in the molecule and the formation of a S−C bond between the molecule and the upper surface (Fig. 6c). The second chemical bond forms between the Zn atom in the molecule and a C atom of the lower surface, resulting in the formation of a crosslinked configuration that bridges the tribogap at a slightly larger P z % 0:3 GPa (Fig. 6d). At P z ¼ 5:1 GPa, breaking of another Zn −S bond and formation of another S−C bond are observed (Fig. 6e). As the contact is reopened, the cross-linked ZDDP molecule is stretched (corresponding to negative P z ; grey line in Fig. 6b). This mechanical pulling induces breaking of a S−C bond on the lower surface and reforming of a Zn−S bond in the molecule, and eventually fragmentation of the ZDDP molecule into two smaller products (Fig. 6f). It is important to note that reactions of ZDDP on surfaces differ from those observed during simulations in the gas phase, where Zn−S bonds repeatedly break Bond-breaking and -forming processes in ZDDP are strongly affected by the contact pressure but not by the chemical structure of the DLC surfaces. Figure 6g shows the contact-pressure dependence of chemical reactions of ZDDP squeezed quasistatically between two a-C and two a-C:H surfaces. The values   6 Quasi-static contact-closing/opening simulations for a-C:H with ρ = 2.0 g cm −3 and CH = 20 at% and a-C with ρ = 2.0 g cm −3 in contact with a ZDDP molecule using DFTB3 quantum chemical calculations. a Initial atomic configuration for a representative a-C:H system. b Evolution of the normal pressure P z as a function of the z system dimension L z during the closing (black circles) and opening process (grey circles). c-f Evolution of atomic configurations for a representative trajectory of a-C:H. For C and H, different colours are assigned between the surface atoms (grey for C and white for H) and atoms in the ZDDP molecule (green for C and cyan for H). Zn, S, P and O are represented by violet, yellow, orange, and red spheres, respectively. The L z values and corresponding normal pressures of the snapshots c-f are indicated by labels in b. g Averaged numbers of S−C bonds Δn SÀC formed between ZDDP sulphur and DLC carbon atoms as well as h sulphur atoms released from ZDDP n Sulphur to the DLC after reopening the contact as a function of the contact pressure P z for a-C:H with ρ = 2.0 g cm −3 and C H = 20 at.% (red) and a-C with ρ = 2.0 g cm −3 (blue). The error bars represent standard error of the means. The vertical lines mark the effective local contact pressures P l;eff for the five experimental coatings (Fig. 1f).
hn SÀC i on the vertical axis represent the average number of bonds between sulphur and DLC carbon atoms that remain intact after reopening the contacts. Our simulations reveal that the formation of S−C bonds is the most dominant chemical reaction, in agreement with our XPS results (Fig. 2). While hn SÀC i increases with P z (0:5 ≤ P z ≤ 10 GPa), the pressure dependence of hn SÀC i in a-C:H and a-C is almost identical. For both DLC surfaces, chemical reactions of ZDDP follow the scenario described in Fig. 6. At P z ≲2 GPa, most ZDDP molecules chemisorb without decomposition (Fig. 6c) or decompose into two fragments by breaking of Zn−S bonds. A further increase of the contact pressure induces decomposition of ZDDP into atoms and small fragments.
Only sulphur atoms that are effectively released onto the DLC surface are relevant to the weakening of the DLC presented in Fig. 5. Therefore, we next calculate the average number of sulphur atoms hn Sulphur i that are not bound to any atoms in ZDDP and remain chemisorbed on the DLC surface after reopening the contact (Fig. 6h). As already observed for hn SÀC i; hn Sulphur i increases significantly with P z and no difference between a-C:H and a-C is detectable (Fig. 6h). Even though the a-C systems used here are representative of the surface region of ta-C, which typically has a low sp 3 content of about 10% 33 , we performed a comparison between ta-C with high sp 3 percentages p sp 3 (ranging from 48−67%) and a-C with p sp 3 ¼ 2 À 9% and observed that hn SÀC i and hn Sulphur i are not significantly affected by the sp 3 content of the hydrogen-free DLC surface (Supplementary Fig. 7).
Next, we perform QMD sliding simulations to explore the role of shearing (Fig. 7). A DLC/ZDDP/DLC tribointerface is studied at a sliding velocity of 30 m s −1 and various P z values (see snapshots for P z = 10 GPa in Fig. 7b, c). After 0.2 ns sliding, the two DLC surfaces are separated quasi-statically to ensure the consistency with the QMS results. For both a-C and a-C:H, hn Sulphur i increases about 2−3 times at all P z compared with hn Sulphur i in QMS simulations and the effect of the DLC's chemical structure on the effective number of sulphur atoms left on the surface is not significant (except for P z = 2 GPa) as shown in Fig. 7a.
Thus, the chemical difference between a-C:H and ta-C, i.e. presence or absence of hydrogen, hardly affects the kinetics of ZDDP's decomposition. The local contact pressures alone determine sulphur's chemical state on the different DLC surfaces. However, the DLC's chemical structure affects shear accommodation at the DLC/ZDDP/DLC interface at high contact pressures. Similar surfaces functionalized with ZDDPderived fragments form at small contact pressures (P z ≲5 GPa) on a-C:H and a-C ( Supplementary Fig. 8), whereas a clear structural difference is observed at P z = 10 GPa (Fig. 7b and c). While a layer of ZDDP-derived fragments still accommodates shear on a-C:H, the a-C surfaces are cold-welded ( Fig. 7c and Supplementary Movies 1 and 2). In this case, ZDDP-derived atoms/fragments diffuse underneath the surface as a result of mechanical mixing. Figure 8 provides further statistical evidence that the chemical difference between the DLC surfaces (i.e. the hydrogen concentration) has a significant effect on the shear response of the frictional interface. First, the weakest interface in each final configuration after MD sliding is determined and then the number of cold-welding C−C bonds n CW CÀC within this interface is calculated. Averaging over different MD runs yields the mean density of cold-welding C−C bonds hn CW CÀC i for a-C:H and a-C (Fig. 8a). hn CW CÀC i is almost zero at P z ≲2 GPa for both a-C:H and a-C, indicating the presence of a sliding interface with no significant cold welding. However, hn CW CÀC i increases drastically for a-C at P z ≳5 GPa, which is a typical contact pressure range reached only by ta-C(66) and ta-C(78) in our experiments (vertical lines in Fig. 8a). In this case, cold welding and shearinduced plastic flow increase the atoms' mobility and induce chemical mixing of the contaminated a-C surfaces. Figure 8b and c clearly shows that the depth distribution (z-axis profile) of the sulphur density at P z = 10 GPa is significantly broadened for a-C, while sulphur is still localised at the surfaces for a-C:H. is imposed along the x direction using Lees-Edwards boundary conditions. The system temperature is kept constant at T = 300 K. a Averaged numbers of sulphur atoms n Sulphur released from ZDDP to DLC after reopening the contact as a function of the contact pressure P z : The error bars represent standard error of the means. The vertical lines mark the effective local contact pressures P l;eff for the five experimental coatings (Fig. 1f). As references, n Sulphur for quasi-static calculations with no shear (Fig. 6h)

Discussion
In this study, a combination of experiments and multiscale modelling reveals the crucial impact of both contact mechanics of surface asperities and DLC chemical structure on the tribochemistry of ZDDP decomposition and, in turn, on the different macroscopic tribological behaviours of DLC coatings in boundary lubrication with a ZDDP-additivated base oil (see a schematic representation of our findings in Fig. 9). While the macroscopic cylindrical Hertzian contact pressure is identical in all our experiments (P macro ¼ 210 MPa), the nanoscale mechanical response of surface asperities to the normal load varies with contact modulus E * and RMS slope h 0 rms ¼ ∇h j j 2 of the  coatings. According to Persson's contact mechanics theory 35,36 , the average local contact pressure for small P macro is given by P l ¼ h 0 rms E * % 0:59h 0 rms E * (irrespective of P macro ). A linear fitting of P l;eff obtained from our contact mechanics calculations gives a slope of 0.4 that is lower than the theoretical one. While this can be probably attributed to the lower degree of self-affinity of our DLC surfaces with respect to the theory and to overestimation of h 0 rms due to the presence of holes ( Supplementary Fig. 1), Persson's theory provides the rationale behind the results of our contact mechanics calculations. Interestingly, the ordering in P l;eff between a-C and ta-C(51) shows that P l;eff depends on both contact stiffness and surface topography. While ta-C(51) is much stiffer than a-C, h 0 rms for ta-C(51) is almost half of that for a-C. This smaller h 0 rms for ta-C(51) can decrease P l;eff compared with a-C, which is in excellent agreement with the experimental S−C/S −Zn ratios (Fig. 2d).
Our results suggest that the effective local contact pressure P l;eff is a crucial parameter to determine the kinetics of ZDDP decomposition as well as the wear and friction regime. For soft DLCs (such as a-C:H, a-C and ta-C(51)), at P l;eff ≲3 GPa, ZDDP chemisorbs on the surfaces by breaking of Zn−S and formation of Zn−C bonds and undergoes decomposition into zinc-and dithiophosphates (Fig. 9a). An anti-wear tribofilm, consisting of zinc polyphosphates, can grow via accumulation and subsequent cross-linking of such precursors 21 . During the tribofilm growth, the formation of a narrow shear zone in ZDDP-derived layers inhibits DLC's weakening by mechanical mixing. Eventually, nanoscale patchy anti-wear tribofilms form on the surfaces. Conversely, P l;eff exceeds 9 GPa on hard ta-Cs and such high local contact pressures aggressively promote ZDDP fragmentation (Fig. 9b). Mechanochemical sulphur doping of the sub-surface region becomes prominent and weakens the carbon matrix via shear-induced mixing, resulting in massive wear of ta-Cs. The difference in wear volume between a-C:H and hard ta-Cs is consistent with previous experiments 5 . Interestingly, P l;eff values ( Fig. 1f) are consistent with contact pressures for single-asperity DLC/steel sliding contacts where tribofilm growth was observed 16 . Since P l;eff is estimated as the median of all the local contact pressures in the real contact area from contact mechanics calculations, and thus characteristic of asperity-asperity contacts, this consistency suggests that the use of P l;eff can fill the gap between atomic and macroscopic experimental scales. P l;eff is determined by the contact parameter E * h 0 rms . Therefore, a direct link between coatings' properties and ZDDP's tribochemistry, which in turn governs the DLCs' macroscopic tribological behaviour, can be established. When E * h 0 rms exceeds a critical value (~25 GPa), tribofilm formation is replaced by sulphur doping (Fig. 9).
In addition to contact pressure, shear stress plays a decisive role in lowering the activation barriers for ZDDP chemical reactions and increasing the accessibility of reactive sites. Moreover, shear action at high contact pressures (P l >5 GPa) makes the effect of the DLC's chemical structure (i.e. presence or absence of hydrogen) apparent. At low contact pressures (P l <5 GPa), irrespective of presence or absence of hydrogen, shear can be accommodated in a narrow region consisting of ZDDP-derived species. Thus, since P l;eff ≲3:5 GPa for a-C:H, a-C and ta-C(51), these soft DLC surfaces are likely to remain intact, and the accumulation of ZDDP-derived species can lead to the growth of anti-wear tribofilms irrespective of the chemical differences between the DLC coatings.
However, at P l >5 GPa, for a-C and ta-Cs, the separation between two surfaces is insufficient to localise shear in the ZDDPderived tribolayer, and cold welding via C−C bonds (bridging the two surfaces) takes place. Under shear, mechanical mixing at the cold-welded tribointerface induces sulphur doping into the subsurface regions resulting in significant DLC's weakening and wear. This process is mechanically driven and leads to a nonequilibrium state with low shear resistance that accommodates the shear deformation. Like in other mechanically driven tribological processes 37,38 , this metastable material has a higher energy than the original one and sulphur penetration cannot be described thermodynamically. In contrast, on a-C:H, the ZDDP chemisorbed layer would be able to support P l >5 GPa owing to steric hindrance between surface hydrogen atoms. Although such high local pressures are not reached with our a-C:H, they would be possible for a hard ta-C:H.
In conclusion, this study sheds light on hitherto unknown multiscale physical processes and reveals how the interplay of mechanical, topographical and chemical factors determine the tribological behaviour of DLC coatings under lubrication with a ZDDP-additivated base oil. We find a knowledge-based design recipe for DLC coatings (hydrogen-free and of moderate hardness) that are wear-resistant in ZDDP-additivated oils while still allowing ultralow friction. Furthermore, the insights into ZDDP's fundamental mechanochemistry and the interplay of mechanics and chemistry could be transferable to ferrous or other nonferrous/metallic lubricated systems. Hopefully, this study will pave the way for further studies aiming at finding FMs and lubrication configurations/conditions that improve the wear and friction performance of DLC/ZDDP systems and, more in general, at controlling tribochemistry of additives by tuning microscopic contact mechanics and surface chemical structure. Such studies could facilitate the formulation of alternative, less harmful anti-wear additives that protect iron-based and DLC surfaces simultaneously.

Methods
Materials. Five types of DLC coatings, labelled by a-C:H, a-C, ta-C(51), ta-C(66) and ta-C(78), were prepared by the company HEF IREIS. The numbers in the parentheses represent the hardness values of ta-C in GPa. Chemical compositions, mechanical and structural properties, and roughness parameters of the DLC coatings are tabulated in Table 1. The coatings had a thickness of 2 μm and were deposited onto mirror-polished steel discs and steel cylinders. The steel discs had a diameter of 25 mm and a thickness of 7 mm, and was made of M2 steel (0.85%C, 6%W, 5%Mo, 4%Cr, 2%V, hardness 64 HRC). The 100Cr6 steel cylinders had a radius of 5 mm and a length of 10 mm. For both geometries, an intermediate layer was applied onto the steel surface to improve adhesion of DLC films. In order to obtain a smooth surface and remove surface defects, the DLC-coated discs and cylinders were polished prior to tribological tests. The surface roughness was measured with a tactile profilometer using a filter cut-off of λ c ¼ 80 μm. The mechanical properties of the materials were determined by nanoindentation using a Nanoindenter XP ® . Indents were done with a diamond Berkovich indenter (tetrahedral geometry, 115.12°between edges). The maximum applied load was 450 mN assuming a Poisson's ratio of 0.20. Continuous stiffness measurement method was used. An oscillating displacement with a frequency of 32 Hz and an amplitude of 1 nm was superimposed at the continuous displacement. The tests were performed using a "constant strain rate" procedure. The ratio of the loading rate to the load was kept constant and was equal to 3 10 À2 s −1 . With a Berkovich indenter, the strain is constant during loading experiments 39 . Seven tests are made for each sample. Simultaneous measurement of the applied normal load, the contact stiffness versus displacement was recorded. Thus, hardness H (related to a contact pressure) and reduced elastic modulus E * were continuously calculated and analysed versus plastic depth. The mechanical properties of the DLC coatings on the cylinders are assumed to be the same as those on the discs. PAO-4 as a base oil and primary-C4 ZDDP (1 wt.%, sourced from Total, France) were used. We note that the concentration of 1 wt% is typical for commercial engine oils and that NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24766-6 ARTICLE NATURE COMMUNICATIONS | (2021) 12:4550 | https://doi.org/10.1038/s41467-021-24766-6 | www.nature.com/naturecommunications primary C4 ZDDP is more thermally stable 40 and its tribofilm growth rate is slower than secondary ZDDP 15 . The kinematic viscosity of PAO-4 is 3.8 cSt at T = 100°C.
Tribological tests. Cylinder-on-disc sliding tests were performed using a homemade linear reciprocating tribometer 41 . All tests were carried out using self-mated DLC/DLC tribopairs. The contacting surface between the cylinder and disc (flat coupon) was immersed in 3 ml of lubricant. A good contact between the disc and cylinder is ensured for every test thanks to an auto-alignment spring system of the cylinder holder. The friction coefficient µ is calculated as the average of the instantaneous friction coefficient over one complete cycle (1000 points per cycle). The testing conditions used were as follows: applied load F N ¼ 23 N (corresponding to an initial maximum cylindrical Hertzian contact pressure P Hertz of 210 MPa), temperature T = 110°C, average speed v ¼ 0:1 m s −1 , stroke length l ¼ 10 mm, and test duration t ¼ 90 min (27,000 cycles). P Hertz ¼ 210 MPa is the same for all DLC coatings. Since the DLC coatings are only 2 μm thick and much thinner than the steel substrate, the macroscopic contact pressure is dominated by mechanical properties of the steel substrate. Indeed, the diameter of the wear scars is equal to that predicted by the Hertz theory. In this test configuration, the lubrication regime corresponds to boundary lubrication (BL) at low speeds, and transitions into mixed lubrication (ML) at higher speeds, as the thickness of the lubricant film increases (see estimation of EHL film thickness and a detailed discussion about the lubrication regime in Supplementary Note 9). As shown in the following section, sliding induced local temperature rises in the BL sections of the stroke stay well below 5°C.
Calculation of flash temperatures. In our reciprocating sliding tests, flash temperatures at the sliding interface do not affect ZDDP's decomposition and macroscopic friction and wear. For rectangular contacts (i.e. cylinder on flat) 42 , the flash temperature rise (T f ) at a dimensionless position X on a homogeneous body subject to a uniform heat flow q ¼ μWjv disk Àv cylinder j 4bl is calculated by where l and b are the half width and length of the rectangular contact, respectively, K is the thermal conductivity, μ is the friction coefficient, W is the normal load, and v is the sliding speed of the homogenous body relative to the contact. The thermal diffusivity κ is calculated from κ ¼ K=ρC p , where ρ is the mass density and C p is the heat capacity. The dimensionless parameters X, B, and L are defined as vx 2κ , vb 2κ and vl 2κ .
x is the position on the body along the sliding direction. First, we consider the temperature rise of a homogeneous body consisting of steel and correct for the presence of a thin DLC coating in a following step. The material parameters for 100Cr6 steel ρ steel ¼ 7:704 g cm −3 , K steel ¼ 26:1 W m −1 K −1 , and C p;steel ¼ 0:446 J g −1 K −1 are taken from the literature 43 .
In the following, we choose a reference frame with a stationary cylinder in contact with a moving disk. The integral in the formula for T f X ð Þ is evaluated numerically for the stationary situation (v ! 0) providing the maximum temperature rise T fmax;cylinder on the cylinder as well as for velocities jvj>0 providing the maximum temperature rise T fmax;disk on the disk. Now, we correct both flash temperatures for the presence of a coating. For a layered body consisting of a thin film (DLC) and substrate (steel), the maximum flash temperature for DLC T fmax;DLC is calculated from the ratio of T fmax;DLC to that for the homogeneous steel T fmax;steel 44 : for 0:1 ≤ Pe DLC ≤ 10:0 for Pe DLC <0:1, According to articles by Robertson and his co-workers 45,46 , the correlations between the density ρ and Young's modulus E of a DLC coating are given by ρ ¼ 1:37 þ E 2=3 =44:65 for a À C and ta À Cs; ð7Þ ρ ¼ 0:257 þ 0:011 E þ 511 ð Þ =44:65 for a À C:H: The heat capacities C p are taken from the literature 47 . The thermal conductivities K for the DLCs are estimated by a linear fitting of the relation between experimental values of K and E 48 . A friction coefficient μ = 0.02 is used for ta-C(51) and μ = 0.1 for the other coatings. T fmax is defined as a maximum T f during the stroke.
At the end, the maximum flash temperature rise T fmax must be the same for both stationary cylinder and moving disk, and can be calculated from 49 .
Most relevant are the flash temperature rise in the BL regime (i.e. at low sliding speeds near the edges of stroke) where reactions of ZDDP happen and thus our simulation models are valid. The detailed results are given in Supplementary Note 10. An average maximum flash temperature in the BL regime (18% of the stroke) T fmax;BL is 2.3°C for a-C:H. For ta-C(51), the BL regime accounts for 32% of the stroke and T fmax;BL is 0.5°C. The temperature rise is negligibly small for all DLCs, and indeed no correlations of T fmax;BL with DLC's properties and observed phenomena were found. For example, for a-C:H, tribo-patches form predominantly near the ends of the stroke while T f is larger at the middle of stroke (i.e. at higher sliding speeds) and regardless of the local temperature the hard ta-Cs undergo much more wear.
Basic surface analysis. After the friction tests, the discs and cylinders were cleaned twice in an ultrasonic bath with n-heptane (Chimie Plus: >99%) for 10 min to remove residual oil. The sample was handled with metallic tweezers without touching the surface of the disc. The wear volumes of the disc and cylinder were measured using an optical white light interferometer (Contour GT-K1, Bruker). The surfaces were observed without any conductive coating by SEM using a FEI XL30-FEG equipped with an Everhardt-Thornley secondary electron detector and operating under high vacuum. The acceleration voltage was set between 2 and 5 kV. Chemical composition analyses were carried out by EDX using an Oxford Instruments X-max silicon drift detector (80 mm² ultra-thin window). Quantitative analysis of the EDX spectra was performed using the Aztec software.
X-ray photo electron spectroscopy. XPS analysis was performed on the discs using an ULVAC-PHI Versa Probe II spectrometer equipped with a monochromatic Al Kα X-ray source with a beam diameter of 200 μm. The binding energy scale was calibrated with respect to the C1s photo-peak at a binding energy of 284.8 eV. The error of binding energies is estimated as ±0.1 eV. First, a survey spectrum was performed using a pass energy of 187.85 eV to identify all the elements present on the surface. Afterwards, narrower scans with a range of 20 eV were acquired using a pass energy of 23.5 eV to accurately identify the chemical state of each element and to perform quantitative analysis using PHI Multipack software. The contribution of the background was approximated by the Shirley method and Wagner sensitivity factors (corrected for the transmission function of the apparatus) were used for the calculation of the atomic concentrations. An Argon ion beam was used to sputter the DLC surfaces at an angle of 45°. The ion beam was accelerated at 250 V with a current density of 120 nA mm -2 . We sputtered the surfaces for 5 min. The corresponding sputter depth is~0.5 nm 28 .
FIB-TEM analysis. Dual Beam Focused Ion preparation was performed on a thin cross-section of the wear track of the disc employing Manutech-USD. A platinum gaseous precursor was first used to deposit a platinum layer onto the DLC surface, to protect it during milling and avoid any re-deposition. This layer was deposited in two steps. First, a low-energy electron beam was used to avoid damaging the DLC surface. Afterwards, a Ga ion beam was used at high energy to increase the growth rate of the protective platinum film and accelerate the process. The a-C:H thin foils obtained were analysed by TEM using a Jeol 2010F-UHR TEM equipped with a Schottky field emission gun operating at 200 kV. Bright-field images were acquired using a Gatan Orius 100 CCD camera and HAADF images in the Scanning Transmission Electron Microscopy (STEM) mode. EDX (80 mm² silicon drift detector from Oxford Instruments) was used for local chemical analysis (spectra or elemental maps). The ta-C(66) thin foil was analysed using a FEI Titan ETEM G2 electron microscope operated at 300 keV and equipped with a Cs image aberration corrector. To avoid contamination prior to the analysis, the samples were plasma-cleaned with Argon for 30 s. Compositional analyses were conducted in TEM mode with an energy dispersive X-ray spectrometer (SSD X-max 55 mm² from Oxford Instruments).
Contact mechanics calculations. Elastic-plastic contact of two rough DLC surfaces were modelled by using experimental AFM topography maps with discrete height profiles h x;y over a scan area of 20 × 20 μm 2 (including 512 × 512 data points). Local contact pressures p x;y and elastic displacements u x;y were calculated under an external normal pressure of 210 MPa (equal to a maximum cylindrical Hertzian contact pressure used in our experiments) using the PyCO code developed by Pastewka and co-workers. Plasticity of surface asperities was taken into account if the local contact pressure p x;y went beyond the material hardness 50 . However, no significant differences in P l;eff were observed between elastic-plastic and full elastic calculations. The details of numerical techniques used in our contact mechanics calculations are described somewhere else [51][52][53][54][55] . Young's moduli and hardness values are tabulated in Table 1. A typical value of 0.2 was used for the Poisson's ratio of DLC. For each pair, 25 contact mechanics calculations were carried out by moving a topography with respect to the other along the x axis with a step of 20 pixels. We assumed that the load is entirely supported by asperities and no fluids are entrained in the contact zone. Thus, P l;eff is only valid for boundary lubrication, i.e. at the ends of the stroke. We employed AFM topographies measured inside the wear tracks after running-in since they are representative of DLC surfaces in a steady, low-friction state. A contact mechanics calculation of DLC topographies measured outside wear tracks showed that significant amounts of surface asperities undergo plastic deformation (corresponding to an initial high friction in the Fig. 1b). However, the high friction state lasts only for the first few hundred cycles and µ drops immediately. For example, for a-C:H, the RMS heights and slopes outside the wear track are 2−3 times larger than those inside the wear track after running-in, which results in a much larger P l;eff , and no correlation between P l;eff and E * h 0 rms was found.
Quantum chemical calculations. Bond-breaking and -forming of ZDDP in contact with DLC surfaces were studied using third-order density-functional tightbinding (DFTB3) calculations 30 as implemented in the Atomistica software suite 56 . Atomic forces are extracted from quantum chemical calculations in order to describe complicated bond-breaking and -formation processes under tribological conditions. The accuracy of Slater-Koster parameters for modelling the interactions between ZDDP and DLC surfaces is confirmed by first-principles densityfunctional theory (DFT) calculations ( Supplementary Fig. 10).
In this study, three types of DLC surfaces were considered: a-C:H with the density ρ = 2.0 g cm −3 and hydrogen content C H = 20 at%, a-C with ρ = 2.0 g cm −3 , and ta-C with ρ = 3.0 g cm −3 . The aims of the simulations in Figs. 5 and 6-8 are different. In Fig. 5, the bulk shearing of sulphur-doped DLCs aims to understand weakening of the bulk regions under shear. Thus, the three computational models were chosen to exactly correspond to experimental a-C:H, a-C, ta-C(78), respectively. In contrast, in the simulations of Figs. 6-8, we considered the former two DLC models due to the following reasons. Our focus is on understanding tribochemical reactions of ZDDP on DLC surfaces. In general, hydrostatic compression of a-C increases its density and sp 3 percentage in the bulk and as a result forms ta-C with a high sp 3 percentage p sp 3 38 . However, a superimposed shear motion drives the system away from its equilibrium phase. At the sliding interface, a ta-C surface layer undergoes sp 3 -to-sp 2 rehybridisation resulting in a significant decrease of p sp 3 33 . This tribolayer with p sp 3 % 10% should be universal for any hydrogen-free DLC coatings and therefore we modelled the reactive zone on top of the ta-Cs as a low-density a-C.
All DLC samples were generated by quenching a melt from 5000 to 0 K at a constant rate of 1 K fs −1 . For Fig. 5, sulphur-doped ta-C samples with initial densities ρ = 2.0 and 3.0 g cm −3 and a same initial cell size of 14 × 14 × 14 Å 3 were produced as starting configurations. For the DLCs in Figs. 6-8, three samples with a cell size of 15 × 15 × 10 Å 3 were generated. For each sample, three surfaces were then created by cutting them perpendicular to the z-axis at different z coordinates, separating the resulting slabs by 2 nm and introducing a ZDDP molecule (with ethyl groups) into the created vacuum region. Before inserting the ZDDP molecule, the DLC slabs were fully relaxed. Thus, in total, nine independent quasi-static molecular statics contact-closing/reopening trajectories were generated for each type of DLC. During QMS simulations, the cell size L z in the z-direction is decreased with a step 4L z ¼ À0:2 Å and once the system reaches a desired contact pressure P 0 , L z is increased to the initial one (L 0 ¼ 30:0 Å) in steps 4L z ¼ þ0:2 Å. For each step, the entire system is relaxed so that all force components are below 5:0 10 À3 eV Å −1 with the FIRE algorithm 57 . This quasi-static contact-closing/ opening simulation mimics an experimental situation in which two asperities come into contact and then detach from each other in the limit v→0 m s −1 . However, the effect of shear stress is not explicitly taken into account within this quasi-static simulation.
MD shearing simulations of two systems (ZDDP/DLC and sulphur-doped ta-C) were carried out with Lees-Edwards boundary conditions 31 . The system temperature T was kept constant at 300 K using a Peters thermostat 58 , and the equations of motion were integrated with a time step Δt = 0.5 fs using the velocity-Verlet algorithm 59 . The system pressure was controlled using a Berendsen barostat 60 . The systems were sheared at a constant speed along the x axis, and the averaged shear stress τ was calculated from the component τ zx of the stress tensor.
Although boundary conditions in our atomistic simulations were matched to local experimental conditions (e.g., temperature and local contact pressures), the sliding speed of v MD ¼ 30 m s −1 in MD simulations is two orders of magnitude larger than the experimental values (v exp ¼ 0:157 m s −1 ), as in previous MD studies [61][62][63] . The use of such high sliding speeds is necessary to simulate a long sliding distance, sufficiently sample phase space, and especially describe tribochemical reactions and subsequent phase transitions at sliding interfaces 2,61 .
However, as long as the sliding speed is well below the speed of sound in solids 64 , the heat generated at sliding interface can be rapidly dissipated to surrounding bodies (which is often modelled by coupling the thermostat to the system). This is a rationale for reliable MD modelling of friction of materials. In our study, v MD is two to three orders of magnitude smaller than the speed of sound in DLCs. Indeed, we did not observe any artefacts of the higher sliding speed on the motion of atoms.
In addition, when the temperature is low enough, thermal activation of chemical reactions is negligible. This means that the chemical reactions of ZDDP are dominated by the competition between the mechanical strengths of the different chemical bonds and independent of the sliding speed. Thus, the observations in our MD simulations should be transferable to our low-speed experiments. This argument is confirmed by the observation that our MD shearing simulations reveal similar results as our quasi-static calculations.

Data availability
Source data generated in this study are provided in the Source Data file. Data generated in this study have been deposited under the Zenodo repository at https://doi.org/10.5281/ zenodo.4898983. Source data are provided with this paper.

Code availability
The density-functional tight-binding calculations were performed using the open-source software Atomistica. The source code is available under GNU GPL v2 at http://www. atomistica.org. The density-functional theory calculations were carried out using the open-source software CP2K. The source code is available under GNU GPL v2 at https:// www.cp2k.org. The contact mechanics calculations were carried out using the opensource software PyCo. The source code is available under MIT license at https://github. com/ComputationalMechanics/ContactMechanics.