# Modeling in tribology: Recent advances, applications, and open questions

Lars Pastewka, Antonis Vakis, Ramin Aghababaei, Andreas Almqvist, Giuseppe Carbone, Michael Chandross, Daniele Dini, Stefan Eder, Hendrik Ehrich, James Ewen, et al.

► To cite this version:

Lars Pastewka, Antonis Vakis, Ramin Aghababaei, Andreas Almqvist, Giuseppe Carbone, et al. Modeling in tribology: Recent advances, applications, and open questions. Tribology International, 2025, pp.111326. 10.1016/j.triboint.2025.111326. hal-05397999

HAL Id: hal-05397999

https://hal.science/hal-05397999v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Modeling in Tribology: Recent Advances, Applications, and Open
Questions

Lars Pastewka<sup>a</sup>, Antonis I. Vakis<sup>b</sup>, Ramin Aghababaei<sup>c</sup>, Andreas Almqvist<sup>d</sup>, Giuseppe Carbone<sup>e,n</sup>, Michael Chandross<sup>f</sup>, Daniele Dini<sup>g</sup>, Stefan J. Eder<sup>h,i</sup>, Hendrik J. Ehrich<sup>h,i</sup>, James P. Ewen<sup>g</sup>, Nicola Menga<sup>e</sup>, Jean-François Molinari<sup>j</sup>, Gianpietro Moras<sup>k</sup>, Lucia Nicola<sup>l</sup>, Marco Paggi<sup>m</sup>, Carmine Putignano<sup>e</sup>, Michele Scaraggi<sup>o,p</sup>, Vladislav A. Yastrebov<sup>q</sup>, Martin H. Müser<sup>r,*</sup>

**Abstract**

Recent advances in modeling have enhanced our ability to make quantitative predictions for tribologi-
cal phenomena, thereby unraveling relevant mechanisms. Algorithmic innovations, including those based
on multiscale methods and machine learning, have been especially impactful, for example in overcoming
long-standing bottlenecks that hinder simulations of systems with strong coupling across disparate scales.
However, traditional modeling approaches, such as boundary-element techniques, have also progressed and
continue to yield new insights. This article reviews developments from the past decade, examining how both
new and established methods have deepened our understanding of experimental results and have furthered
theoretical approaches in key tribological areas, including contact mechanics, lubrication, metal friction,
and tribo-chemistry. Selected applications, such as tunable interfaces and energy harvesting, illustrate the
broad influence of recent developments on the field.

Keywords: modeling; tribology

<footnote><sup>a</sup>Department of Microsystems Engineering, Faculty of Engineering, University of Freiburg, Georges-Köhler-Allee 103,
D-79110 Freiburg, Germany</footnote>

<footnote><sup>b</sup>Computational Mechanical and Materials Engineering, Faculty of Science and Engineering, University of Groningen, 9747
AG Groningen, The Netherlands</footnote>

<footnote><sup>c</sup>Department of Mechanical and Production Engineering, Aarhus University, Katrinebjergvej 89F, 8200 Aarhus N, Denmark</footnote>

<footnote><sup>d</sup>The Division of Machine Elements, Luleå University of Technology, Luleå SE97187, Sweden</footnote>

<footnote><sup>e</sup>Department of Mechanics, Mathematics and Management, Polytechnic University of Bari, Via Orabona 4, Bari 70125, Italy</footnote>

<footnote><sup>f</sup>Computational Materials Science, Sandia National Laboratories, Albuquerque, New Mexico 87123, United States</footnote>

<footnote><sup>g</sup>Department of Mechanical Engineering, Imperial College London, South Kensington Campus, Exhibition Road, London
SW7 2AZ, United Kingdom</footnote>

<footnote><sup>h</sup>AC2T research GmbH, Viktor-Kaplan-Straße 2/C, 2700 Wiener Neustadt, Austria</footnote>

<footnote><sup>i</sup>Institute for Engineering Design and Product Development, TU Wien, Lehárgasse 6 – Objekt 7, 1060 Vienna, Austria</footnote>

<footnote><sup>j</sup>Institute of Civil Engineering, Institute of Materials Science and Engineering, École Polytechnique Fédérale de Lausanne
(EPFL), Lausanne, Switzerland</footnote>

<footnote><sup>k</sup>Fraunhofer IWM, MikroTribologie Centrum µTC, Wöhlerstraße 11, 79108 Freiburg, Germany</footnote>

<footnote><sup>l</sup>Department of Industrial Engineering, University of Padova, Padua, Italy</footnote>

<footnote><sup>m</sup>IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy</footnote>

<footnote><sup>n</sup>CNR Institute for Photonics and Nanotechnologies U.O.S. Bari, Physics Department M. Merlin, Via Amendola 173,
I-70126, Bari, Italy</footnote>

<footnote><sup>o</sup>Department of Engineering for Innovation, University of Salento, 73100 Lecce, Italy</footnote>

<footnote><sup>p</sup>Center for Biomolecular Nanotechnologies, Istituto Italiano di Tecnologia, 73010 Arnesano, Italy</footnote>

<footnote><sup>q</sup>Mines Paris PSL University, Centre des matériaux, CNRS UMR 7633, Versailles, France</footnote>

<footnote><sup>r</sup>Dept. of Materials Science and Engineering, Saarland University, 66123 Saarbrücken, Germany</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email address: martin.mueser@mx.uni-saarland.de (Martin H. Müser)</footnote>
# 1. Introduction

As in other scientific and engineering fields, modeling plays a central role in tribology [1], which concerns solid friction, lubrication, and wear, or more generally, all processes that occur when two surfaces are in relative motion. Understanding and ultimately improving tribological systems benefits energy efficiency and sustainability. Nonetheless, applications are still significantly driven by trial-and-error procedures, partly because it is difficult to obtain direct experimental information on stresses, temperature, chemical composition, and related quantities near the tiny zones in a buried interface [2], where surfaces touch and rub. This is where modeling offers a distinctive advantage. It allows complex dynamics to be directly visualized by simulating the mechanical, thermal, and chemical processes locally. However, this is easier said than done, since the simulator ideally needs expertise in mechanics, materials, and chemistry. It also requires knowledge of quantum physics and electromagnetism—for example, when addressing triboelectricity [3]. In addition, incorporating experimental data into models often requires further proficiencies, such as handling open data or even integrating measurements on-the-fly during simulations, as in digital twins.

Fortunately, many interesting tribological problems exist where “single-physics approaches” allow a given question to be answered or a specific target property to be optimized. Examples include assessing the viscoelastic dissipation when an elastomer slides past a textured or randomly rough surface using continuum methods [4], or identifying oils whose viscosity changes as little as possible with temperature through a combination of molecular dynamics and artificial intelligence [5]. However, even these seemingly straightforward problems often become computationally prohibitive—primarily due to the challenge of reaching experimental time scales—or inaccurate due to the limited validity of material models. Matters become even more complex when different physical descriptions are optimal at different length scales. A paradigmatic example is the generation and propagation of a single crack in a tribological coating. Its modeling requires atomic resolution at the crack tip—achievable through molecular dynamics—while boundary conditions are most effectively determined at much coarser scales using continuum-mechanics-based techniques [6, 7].

Scaling up from a single crack to modeling wear in machine components—such as an entire bearing—and ultimately to a full engineering system like a vehicle with numerous moving parts necessitates careful methodological choices. Ideally, modeling should offer added value beyond simply reproducing the you-get-out-what-you-put-in pitfall. Despite the challenge of reaching meaningful time and length scales while maintaining high accuracy, a key advantage of simulations is their ability to examine a moving, buried interface—an aspect that is extremely difficult to examine in *in situ* experiments [2]. For instance, flash temperatures can, in principle, be determined in a sliding contact between non-transparent surfaces before the rubbed surfaces are exposed to air. Atomistic simulations can even address fundamental questions, such as whether the distribution of velocities adheres to the Maxwell-Boltzmann distribution—a prerequisite for defining temperature in analogy to an equilibrium system.

Figure 1 provides a comprehensive overview of the mechanisms governing tribological interfaces across multiple scales. The schematic integrates representations of length and time scale-dependent processes, the underlying mechanisms shaping macroscopic contact behavior, and the inherently multiphysical nature of the problem. At the smallest scale, it explicitly resolves atomic interactions between solid surfaces and any confined medium (e.g., gas, liquid), depicting key elements such as lubricants, additives, oxides, and tribofilms. Moving to the mesoscale, it highlights the role of defects, dislocations, and near-surface microstructure, along with local solid-liquid interactions such as adhesion, wetting, and lubricant film formation, all critical for large-scale behavior.

Figure 1 also serves as a framework for understanding the constellation of phenomena influencing interfacial performance, including damage localization, material transformation, tribolayer evolution, and additive-surface interactions. Additionally, it explicitly captures how external macroscopic influences, such as electric or thermal fields and variable applied loads, directly impact smaller-scale responses. A major challenge in tribology modeling lies in accounting for this vast array of mechanisms and selecting the appropriate tools to capture them. The key lies in determining which aspects must be explicitly included in models. While significant progress has been made in addressing many of these elements through advanced modeling techniques, Fig. 1 serves as a canvas to present the “big picture” by integrating the developments into a unified perspective.
<figure><img src="image_2.png" /><figcaption>Figure 1: A representation of the complexity of a tribological interface capturing multiple scales, mechanisms and external influences that control friction, wear, and engineering performance.</figcaption></figure>

In this review, we highlight recent advances in tribological modeling and summarize key insights gained over the past decade, building on a previous comprehensive review [1] co-authored by many of the authors of this work. We focus on areas where significant progress has been made. These range from AI-assisted real-time dynamic simulation of cross-scale coupling to the long-standing challenge of bridging atomistic simulations with continuum models, encompassing the modeling of bond-breaking and formation in tribo-chemistry, plasticity, or wear under realistic conditions. The main body of this review begins with a broad discussion of methodological advances in Section 2, covering both conventional modeling approaches and developments involving machine learning and artificial intelligence. Section 3 examines recent progress in the continuum mechanics of solids and lubricants—fields with long-standing foundations but ongoing debates on fundamental issues, such as how adhesion contributes to energy dissipation in sliding contacts or whether friction causes contact areas to shrink or expand. Some tribological challenges require materials- or chemistry-specific approaches. These are explored in Section 4, with a focus on the friction and wear of metals and tribochemical processes driven primarily by high local stresses. Application-oriented aspects of tribology are discussed in Section 5, particularly in relation to tunable interfaces, energy harvesting, and batteries. Finally, conclusions are presented in Section 6.

## 2. Modeling

Traditional models for tribological systems mold what we know about mechanics, materials, and chemistry into models suitable to answer specific tribological questions. This requires approximations and assumptions: for example, molecular models for predicting (tribo-) chemical reactions are carried out at some specific external pressure, but which pressure is relevant at molecular scales is determined by surface roughness, which varies widely. Validating a “single-physics” and a “single-scale” model typically requires bespoke experiments: for example, friction measured from atomic-force microscopy (AFM) could be matched with molecular calculations if AFMs had the ability to move at technically relevant velocities. However, for most engineering applications there is still a lack of coupling schemes for arriving at truly macroscopic, “multi-physics” and “multi-scale” predictive simulations. New flavors of models, in particular digital twins and machine-learned models, attempt to fill this gap by combining simulation data with measurements that can
update models on the fly. The hope is that such hybrid models can offer predictive power even in scenarios as complex as those encountered in most mechanical devices. This section highlights advances in the technical aspects of modeling in broader contexts than those covered in later sections. For example, we summarize which artificial intelligence (AI) methods have been applied in tribology, while later chapters provide a more detailed exploration of their use in addressing specific topical problems.

## 2.1. Single-physics models and beyond

Modern modeling tools allow tribologists to study processes across a wide range of scales, from the atomic to the macroscopic. Over recent decades, atomistic simulation techniques like molecular dynamics (MD) and density functional theory (DFT) have provided unique insights into various tribological systems and processes. Molecular simulations can now capture phenomena such as fluid rheology [8, 9, 10], tribochemistry [11, 12, 13], and the formation of rough interfaces [14, 15]. Atomistic models are particularly valuable for revealing molecular processes that often drive the macroscopic behavior of tribological systems [16]. A key limitation of molecular models, however, is the restricted length and time scales, and the chemical complexity [17], they can access. While length scales can be extended through parallelization—allowing simulations with tens of millions of atoms [18]—bridging time scales requires problem-specific, advanced sampling techniques, such as kinetic Monte Carlo [19]. Similarly, despite mounting evidence that tackling chemical complexity is becoming increasingly feasible with machine-learned potentials [20, 21, 22], even much less computationally demanding two-body potentials impose tight bounds on accessible system size and time scales.

Simulating realistic rough interfaces requires system sizes that are typically beyond the reach of molecular techniques. To address this challenge, the tribology community has a rich history of developing mesoscale simulation methods for studying roughness, particularly the boundary element method (BEM) [23, 24], which simplifies the problem by mapping the three-dimensional elastic response of a body with a (typically flat) interface onto the response of the interface alone. Recent advancements include extending BEM to nonuniform grids [25] and tackling nonlinear problems such as plastic flow occurring at and below contacting asperities [26, 27, 28, 29]. BEM has been employed to study wear [30] and the evolution of roughness in contacting interfaces [31]. The largest linear elastic BEM models in recent literature have discretized surface roughness using grids exceeding 65,536 × 65,536 nodes, amounting to over 4 billion degrees of freedom [32, 33, 34]. These high-fidelity calculations serve as reference standards for developing simplified analytical models [32, 33] and facilitate the study of the fractal nature of contacting geometries [34] (see Fig. 10a-c). In situations where molecular information in the boundary itself is still required, BEM models and models with atomistic resolution can be combined [35, 36], even for many-body potentials and finite-thickness boundaries [37].

Another recent advancement is combining BEM with a Griffith-like criterion to study wear particle formation. In this approach, wear particles are generated either when the size of the contact islands exceeds a threshold (see Fig. 10d, Refs. [38, 39]), or when the elastic energy released by removing an asperity contact—calculable within BEM—exceeds the energy required to create new surface area (see Fig. 10e-f, Ref. [40, 41]). When combined with simple geometric criteria for asperity removal, these models provide insights into wear processes. The dimensional reduction inherent in BEM can be extended for adhesive contacts, resulting in one-dimensional models of the contact line. These models enable simulations that bridge length scales over several orders of magnitude [42, 43].

In recent years, two exciting developments have emerged in mesoscale modeling of wear: molecular-dynamics-like models, where beads interact via custom-made pair interactions to create two-dimensional brittle solids at numerically accessible scales [44, 45, 46, 47], and phase-field fracture models [48, 49, 50]. Both types of models overcome a critical size limitation inherent to molecular techniques, which cannot access the length scale of the “process zone”—the region near the crack tip, where the material undergoes plastic deformation. The process zone can extend up to centimeters in ductile metals, and even in apparently brittle materials like glass, it spans micrometers—beyond the reach of molecular simulations. Calculations using a molecular-dynamics-like discrete element model (DEM) have confirmed that a Griffith-like criterion governs the transition between particle formation and asperity flattening: particles form only at length scales larger than a critical value determined by the material's mechanical properties and surface energy [45, 46].
<figure><img src="image_3.png" /><figcaption>Figure 2: Contact morphology and wear models. Contact islands in rough elastic contacts form fractal morphology with complicated branched structures. Panels (a) to (c) show contact islands of different size obtained from a multi-billion degrees-of-freedom boundary-element calculation reported in Ref. [34]. Such calculations have recently been used to formulate wear models. Panel (d) shows the selection of wearing contact islands from Ref. [30] based on a critical junction-size model, where red-islands have junction sizes that would wear. Panel (e) shows a contact morphology in combination with a color-code map (f) of energy released when a hypothetic half-sphere of diameter *D* is removed from the contact, also obtained from boundary-element calculations in Ref. [40].</figcaption></figure>

These DEM models enable the study of system sizes large enough to represent the evolution of self-affine roughness during wear [51]. Phase-field models are similarly suited for studying the formation of wear particles, but reported simulations have so far been limited to single-asperity junctions in elastic [48, 49] and elastoplastic [50] contacts. Finally, coarse-graining techniques utilizing artificial intelligence—such as a recent approach using graph neural networks in a biophysical context [52]—might help parameterize such higher-scale models from first principles.

Mesoscale models that incorporate the structure of surface roughness are also crucial for understanding lubrication. Most lubrication models are formulated in a height-averaged manner, typically leading to the Reynolds equation. Probably the most significant advance in this field over the last decade is the development of quickly converging and stable cavitation algorithms, specifically the Fischer-Burmeister-Newton-Schur (FBNS) algorithm developed by Woloszynski, Podsiadlo, and Stachowiak [53]. Other unconventional methodological developments include solution strategies that are agnostic to the actual constitutive behavior of the fluid [54], making them suitable for learning shear-thinning, shear-thickening, or other nonlinear constitutive laws directly from molecular dynamics simulations [55]. Additionally, the use of smoothed particle hydrodynamics for modeling the flow of lubricants [56, 57] and the material-point method for plastic flow of metals [58] represent significant advances, broadening the modeling tools available to tribologists. Another important aspect in coarse-graining tribological contacts can be the necessity to retain some stochasticity to reflect the coupling of scales. One example is the prediction of leakage rates when fluid flow is blocked by sparsely distributed critical constrictions [59]. Reproducing the desired heterogeneity can be achieved by introducing realistic randomness into the local constitutive laws relating stress and displacement [59, 60].

Simulations at the scale of a tribological component have historically employed continuum methods that do not resolve microscopic features of roughness, such as the finite element method (FEM). However, understanding a component often simultaneously requires precise chemistry (bond breaking and formation), large-scale nonlinear phenomena (plasticity, fracture), and structural geometries of industrial relevance (bearings, gears) subject to thermo-elastic loading [61]. This necessity strongly drives the extension of
available FEM technologies for nonlinear large-scale problems by adding microstructural details with novel discretization techniques. This involves embedding textures and roughness into interface elements, as proposed in Refs. [62, 63, 64]. If there is a scale separation between microstructural features of roughness and the structural elements of the component, concurrent coupling of FEM and the boundary element method (BEM) enables nonlinear structural mechanics simulations of real components with implicit FE solution schemes [65], while incorporating the effect of roughness in the cohesive law that describes contact in the FE simulation. Conversely, FEM can also be useful to turn linear BEM calculations into a full nonlinear solver, where FE calculations are only carried out in contact regions with large slopes where the largest deviations from linear elasticity are expected [66].

BEM can also be exploited to interpret phenomena occurring over multiple scales of observation. Such an approach was attempted by Marulli et al. [67], using elevations acquired by photogrammetry at the macroscale and detailed confocal profilometric height fields in the surface areas where the origin of wear traces needed to be investigated. We expect these and related coupled techniques to impact the description of phenomena like contact-induced fracture, where the solution of the contact problem affects crack growth. Similar to the mesoscale wear problem described above, crack growth (and wear) can be simulated using FEM formulations of phase-field problems [68, 69, 70, 71].

## 2.2. Digital twins

Digital twins are virtual counterparts of machines, systems, or processes [72], which evolve over time and are regularly updated based on measurements from the real-world system. Ideally, they persist throughout the entire life cycle of its physical counterpart. In a tribological context, sensor data on surface roughness gathered during operation can be used to update and improve the accuracy of the digital twin [73]. This is particularly useful for processes that are difficult to capture using conventional models, such as running-in and wear. Digital twins are already used in manufacturing and a wide range of industrial applications, including those where tribology is central to performance, such as wind turbines and aerospace engines [72]. In tribology, digital twins are still largely prospective, but they illustrate how modeling could evolve into continuously updated, predictive tools directly linked to experiments. We here report on a few examples of digital twins that we are aware of.

As a first example, Bucknall [74] recently reported a digital twin of a full tribological system, developed at bp/Castrol for modeling combustion engines. The twin is a virtual replica of the mini traction machine (MTM2), a ball-on-disk tribometer from PCS Instruments (Acton, UK) [74]. Their physics-based model utilized the classical, yet contentious, Greenwood-Tripp contact model [75, 76] for two rough surfaces. The topography of the contact surfaces started from three-dimensional surface profiles of the ball and disk, obtained with white-light interferometry. Repeated measurements of surface topography were then used to validate and update the digital twin for the effect of wear. Several other dynamic models were incorporated to describe the lubricant's density, temperature-viscosity and pressure-viscosity relationships, shear thinning, limiting shear stress, as well as flash heating and asperity friction. Thus, the model was sensitive to both base fluid and additive properties and was capable of reproducing the full Stribeck curve observed in the MTM experiments. Similarly, Hansen et al. [77] recently reported matched pin-on-disk experiments and fluid mechanics calculations based on the Reynolds equation. A key input was the post-experiment topography, which they used in continuum calculations assuming a mixture of hydrodynamic lubrication and a fixed friction coefficient for solid-solid contact. They then "calibrated" their model by smoothing the topography to reproduce the experimental Stribeck curve in the boundary lubrication regime.

Simpler models for digital twins have been developed to model wear progression, typically by applying Archard's law [78] and updating predictions using measurements of wear track depth. Regis et al. [79] applied this approach to a helicopter bush bearing system, where wear was predicted based on a Hertzian calculation of contact stresses that fed into Archard's wear model. This method can be combined with reduced-order dynamical models of the macroscopic device encompassing the tribosystem, allowing the digital twin to be updated based on vibration monitoring [80].

The usefulness of a digital twin can be questioned if the underlying model is fundamentally flawed. For instance, bearing-area models (such as Greenwood-Tripp) fail to account for the effect of long-range elastic deformation, which is just one of many factors that must be considered to accurately model mechanical
contacts [33]. Relying on insufficient models will prevent a digital twin from passing rigorous tests of its predictive ability.

## 2.3. Machine-learned models

Recent advances in computational power have made feasible the use of “non-parametric” models derived from purely statistical considerations, often referred to as machine-learned (ML) models. Physics-based models incorporate physical principles discovered over centuries of human research, such as energy and momentum conservation or the fact that, in the long-wavelength limit, the mechanical response of solids can be described by a few elastic constants. In contrast, data-driven models essentially interpolate between measured or simulated data points—or even extrapolate when the number of parameters is so large that it is impossible to construct data on a fine mesh [81]. When based on probability theory or random processes, these models can estimate uncertainties or even provide the full *a posteriori* probability distribution of potential model outcomes.

The predictive ability of a machine learning model depends sensitively on the feature vector—that is, the chosen input variables—and the database. Together, these determine the dimensionality of the input space and the available data points from which the model extrapolates. Training an ML model involves optimizing a loss function, which in the simplest case is the squared difference between the model’s predictions and the target values. While variants of neural networks (NN) are the most commonly used schemes, other machine learning approaches such as Gaussian processes also exist [82]. The remainder of this section will be structured around generic features of applying such ML models.

A common application of machine learning (ML) models is to serve as surrogates for physics-based models that are quicker to evaluate. The primary reason for using ML in this context is to replicate the physics-based model while saving computational resources. For example, Walker et al. [83] trained neural networks (NNs) on film thickness and friction forces obtained from classical elastohydrodynamic (EHL) simulations. The feature space, or input parameters, of the model consisted of geometric parameters of the machine element. The authors reported that the NNs were more computationally efficient than full EHL calculations.

In a similar vein, Rapetto et al. [84] and later Kalliorinne et al. [85] conducted elastoplastic BEM calculations of the contact between randomly rough interfaces with 256 × 256 nodal points. They trained neural networks (NNs) to predict the real area of contact from these calculations, using a feature vector consisting of standard scalar roughness parameters such as average surface roughness ($S_a$), root-mean-square roughness ($S_q$), skewness ($S_{sk}$), and others. The database comprised 6,144 calculations of the contact mechanical response of topographies with different statistical properties. This dataset provided distinctive outcomes for learning the area-load dependencies for various types of surface topographies, including self-affine randomly rough surfaces, Weierstrass patterns, bi-Gaussian distributions, and exponential models. Construction of such surrogate models for contact mechanics or EHL appears to be starting to become a sport in the tribology community [86, 87, 88, 89].

A major limitation of purely data-driven models is their reliance on large datasets to adequately capture and learn the underlying physical principles. This limitation can be mitigated by incorporating prior knowledge of the physics into the learning procedure. For example, when learning the pressure profile in a lubricated contact, the optimization process can be regularized by adding a penalty to the loss function for solutions that do not satisfy the Reynolds equation. This approach was suggested by Almqvist [90] and later implemented by others using neural networks [91, 92, 93]. Since the solution to the hydrodynamic lubrication problem is represented by the weights (or parameters) of the neural network rather than by values at specific grid points, such procedures are sometimes referred to as “meshless.” However, it should be emphasized that the Reynolds equation is used only for regularization and does not have to be strictly satisfied.

Machine learning models can be useful for recognizing patterns in topographic data of experimental tribological interfaces. In this context, much research effort is directed toward constructing relevant feature vectors, which are statistical descriptors of the rough topography. These feature vectors represent data in a high-dimensional vector space, where each dimension corresponds to a specific feature extracted from the
<figure><img src="image_4.png" /><figcaption>Figure 3: Coupling continuum scales with molecular scales using machine-learning. The continuum solver (a) is agnostic of the constitutive law and requires fluxes as a function of stresses in the system [54]. This relationship is obtained from (b) molecular dynamics calculations, whose results are interpolated with (c) Gaussian process regression. (d) New datasets are added where the prediction error is large. From Ref. [55]</figcaption></figure>

data, enabling quantitative analysis and pattern recognition by ML models. Techniques like principal component analysis (PCA) are often used to reduce dimensionality while retaining the most significant features. The field was pioneered by Stachowiak and co-workers [94, 95, 96], who introduced feature vectors based on discrete wavelet transforms for detecting wear during machine condition monitoring. They employed a linear Support Vector Machine (SVM) to classify adhesive and abrasive wear levels. More recently, Wolski et al. [97, 98] and Yesilli and Khasawneh [99] used simpler multidirectional geometrical parameters, such as mean curvature and fractal dimensions, to distinguish between abrasive, adhesive, and corrosive wear. Sanner and co-workers have defined derivatives (slopes, curvatures) of topographic data as a function of scale [100] and suggested performing classification based on such statistical descriptors [101].

Machine learning when used as surrogates is synonymous to interpolation between given data points. This interpolation can come with prediction for the error made when interpolating, for example when using Gaussian processes. The predicted error allows to generate new data points in regions when the error is high, a scheme that is often called active learning of kriging. Examples include the computation of constitutive relations to whenever the Gaussian process is uncertain about its interpolation between existing data. A paradigm example for the use of such a scheme in a tribological context is the active learning of viscosity laws and slip boundary conditions of lubrication from molecular dynamics calculations [54, 55]. These microscopic laws are incorporated into a macroscopic lubrication solver, yielding the multi-scale simulation framework shown in Fig. 3.

## 2.4. *Open Science*

A review of physics- and data-driven methods in the 2020s would be incomplete without discussing Open Science practices in the field [102]. Open Science seeks to enhance the accessibility of raw tribological data and ensure the reproducibility of experiments and simulations. The first attempt to generate a truly findable, accessible, interoperable, and reproducible (FAIR) [103] dataset for a full tribological experiment
was carried out by Garabedian and co-workers [104]. Their work is accompanied by a complete description
of the experimental workflow in the form of computer-readable data formats.

Documenting tribological research in this manner requires future standardization in how experimen-
tal protocols and topographic data are communicated. Topographic data can now be deposited in the
central repository contact.engineering [105], which provides standardized analysis workflows for such
data—such as computing the power spectral density or autocorrelation function [106]—as well as facilities
for data publication and citation. One recent example of topography data sharing is the Surface Topography
Challenge [107], in which more than 150 scientists collaborated to characterize stochastically identically sam-
ples. The study highlighted discrepancies in data from commodity instruments and developed best-practices
in how topographic is collected. It may serve as a blueprint for how the community can come together to
further develop standards and practices in data collection and sharing.

Another key form of data publication is the availability of open-source codes that adhere to established
quality control procedures, particularly test-driven development and continuous integration. In the field of
contact mechanics, high-quality codes like Tamaas [29] and the contact.engineering [105] ecosystem have
seen significant development in recent years. The commitment to open data and open-source codes not only
underpins reproducible computational tribology, but it also paves the way for collaborative innovation and
enhanced scientific discovery across the discipline.

3. Mechanics

Contact mechanics investigates how forces and displacements are distributed and transmitted at the interface between two bodies in contact. Early approaches neglected the effects of adhesion, interfacial friction, and roughness beyond deterministic shapes, such as flat punches and parabolic indenters. Time-dependent effects due to the presence of a lubricant, the viscoelasticity of the contacting bodies, or interfacial elastic multistability were also long overlooked. Centered on the framework of contact mechanics, this section highlights recent developments related to (1) interfacial adhesion, (2) frictional sliding, and (3) lubricated contacts.

3.1. Adhesion

When two non-conformal elastic bodies come into contact, surface energy is gained locally where true
contact occurs, albeit at the cost of elastic deformation energy, which counteracts contact formation and
reduces the perceived adhesion [108, 109]. Conversely, several hysteresis mechanisms, such as viscoelasticity,
can hinder contact breakage and effectively increase the perceived adhesion. While it is straightforward to
understand qualitatively that viscoelastic contact hysteresis grows with both a larger work of adhesion and
shorter-range surface interactions, quantitatively predicting interfacial dynamics—arising from the interplay
of adhesion, viscoelasticity, and roughness—remains a significant challenge and an active area of research.
This section highlights recent key works in this context and discusses open challenges.

The (dissipative) dynamics of adhesive interfaces can be divided into four categories, which will, however,
be discussed below in inverse order to their listing in this paragraph. First, jump-in and snap-off contact
instabilities [112, 113], which Tomlinson [112] argued are a primary cause of friction. Since the snap-
off occurs at a larger separation than the jump-in, hysteresis remains even when driving velocities are
extremely low. Fig. 4(c) shows a multi-asperity contact and its stress distribution just before the moment
of snap-off. The Tomlinson model complemented the earlier work by Prandtl [114, 115], who emphasized
the general significance of instabilities for solid friction by showing that sliding-induced instabilities lead
to Coulomb friction when thermal effects are negligible, but result in Stokes friction in the (creep) limit of
low velocities at finite temperatures. Prandtl’s model primarily focused on the in-plane motion of atoms,
which brings us to the second category: the quasi-discontinuous motion of contact lines caused by moderate
waviness [116], leading to adhesive hysteresis at low mean contact-line, or, depending on context, crack-tip
velocities. Fig. 4(b) shows contact lines during moments of arrest, while they are unstable in between.
The energy loss associated with a hysteresis loop occurs even in a perfectly elastic medium, because the
kinetic energy produced during the sudden advancement of a crack front is distributed among vibrations of
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Smooth, adhesive indenter continuously sliding over a viscoelastic substrate at varying velocities. The red line shows the stress, while the straight black line separates compressive (top) from adhesive (bottom) stress. Adopted from Ref. [110]. (b) Stable contact line configurations during the retraction of an indenter from a wavy contact; simulation (top) versus experiment (bottom). Adopted from Ref. [43]. (c) Topography of a multi-asperity indenter (top, color) and stress (bottom, color) during snap-off, blue indicating tensile and orange compressive stress. The remaining black-and-white images (d) compare actual contact (dark) and non-contact (light) between experiment (top) and simulation (bottom) at different times t during the retraction process. Adopted from Ref. [111].</figcaption></figure>

disparate/incommensurate frequencies, preventing the local energy density from ever returning to its initial high value without external work. In other words, there is an implicit assumption that the instable degrees of freedom are damped. Third, steadily moving interfacial cracks—or evolving contact geometries—as depicted in Fig. 4(a). They occur when a smooth, rigid, adhesive solid indents into or retracts from an elastomer but also during peeling [117]. Their description is similar to crack initiation, growth, and closure in viscoelastic media, a well-established topic [118, 119], but one that continues to be actively investigated because analytical work can only address relatively simple geometries, while numerical simulations face tremendous challenges due to the broad spectrum of time- and length-scales involved in viscoelastic adhesion. Last but not least, adhesion in hard-matter systems can lead to plastic rather than elastic deformation [120]. Except for this last point, we have largely avoided classifying contacts as soft and tacky versus hard and non-sticky, as this distinction is non-universal and breaks down at small scales and high excitation rates, as might become clearer below.

To help the reader better digest the remainder of this section, we re-iterate some key results on rough contacts and the competition between surface energy and elastic energy [109]. For contacts of micron size or larger to exhibit adhesion (or “stickiness”), a substantial relative contact area must develop without any external load. This requires that the work of adhesion, $\gamma$, be comparable to or greater than the elastic energy per unit area, $v_{\text{ela}}$, needed for the elastic body to conform closely to its counterface. (To be more concise, one needs to consider the combined contact compliance $1/E^*$ and the combined roughness of the two contacting bodies. As usual, $E^* = E/(1 - \nu^2)$ denotes the contact modulus, where E is the Young’s modulus and $\nu$ the Poisson’s ratio.) The cost of perfectly following a single sinusoidal surface undulation, $u(x) = u_q \cos(qx + \varphi_q)$, with wavelength $\lambda = 2\pi/q$, can be expressed as [109, 121, 122]:

$$v_{\text{ela}} = \frac{qE^*}{4} u_q^2 \qquad (1)$$

This relation, which can be seen as the working-horse for any efficient, Fourier-based boundary-element method but also for Persson’s contact mechanics theory [4], follows from simple dimensional analysis, apart, of course, for the numerical prefactor $1/4$. This is because the elastic energy is proportional to stiffness and scales quadratically with the amplitude. The result for $v_{\text{ela}}$ suggests that the critical roughness parameter for a single surface undulation is $h_{\text{adh},q} = qu_q^2/4$ and that generalizations involve summing individual roughness parameters over wave numbers or wave vectors for more complex surface shapes [43]. It is worth noting that the scaling $v_{\text{ela}} \propto q$ holds for a semi-infinite solid can be interpreted as a fractional Laplacian operator acting on the surface deformation [42, 43]. For contact with a membrane (under tension) the scaling relationship is modified to $v_{\text{ela}} \propto q^2$, while for thin sheet (plate) one finds $v_{\text{ela}} \propto q^4$ [123].
<figure><img src="image_6.png" /><figcaption>Figure 5: Canonical model of the power-spectral density (PSD) used in many theoretical works (black line), with a constant PSD from system size at wavevector $q_0$ to the roll-off at wavevector $q_r$, where the PSD crosses over to a power law at short wavelength. The power-law region often extends down to the atomic scale, described by the short-wavelength cutoff $q_s$. For a self-affine topography, the power-law region follows the scaling $q^{-1-2H}$ for the one-dimensional PSDs shown here, where $H$ is the Hurst exponent. The data shown is from a measurement of a microcrystalline diamond film that exhibits this idealized PSD. Black circles show the average over all individual measurement, shown by crosses. The roll-off wavelength $2\pi/q_r \approx 1$ µm corresponds to the film thickness. The experimental data (different colors refer to different techniques) is taken from Ref. [130].</figcaption></figure>

The height spectrum of most real surfaces $C(q)$, which is proportional to the (local) mean or expected value of $|\tilde{u}(q)|^2$, can be approximated as constant for small wave-numbers $q$ up to a roll-off wave vector $q_r$ and then as a power law, $C(q) \propto q^{-D-2H}$ down to very small scales (see Fig. 5 and Refs. [124, 125]). Here, $D$ is the interfacial dimension (i.e., $D=1$ for line and $D=2$ for areal height profiles) and $H$ is the Hurst exponent, which typically takes values below but close to one). Thus, when summing up individual contributions $h_q$ to a net (or scale-dependent) $h_{\text{adh}}$, the near-$q_r$ contributions dominate on both the roll-off and the power-law or self-affine branch of $C(q)$ – at least when $H > 0.5$ [124, 126, 127, 128]. Therefore, regularly stated quantities like mean or root-mean-square roughness, which have their dominant contributions at small $q$, are not useful descriptors in the context of adhesion. Our semi-quantitative discussion allows further conclusions to be drawn, most importantly that length but also time scales matter, e.g., the propensity of a small rigid body to stick to a viscoelastic counterface can decrease substantially with an increasing linear dimension $l$ when $l < 1/q_r$. Likewise, the frequency with which contact is formed and broken also matters since elastomers stiffen at high frequencies.

Inspired by earlier work of Anderson [129] on dislocation loops near crack tips, Sanner et al. [43] suggested an approach with which the local work of adhesion is reduced by the (semi-) local elastic energy needed to form local contact. Their idea is connected to the interpretation of a term like $\sqrt{q\tilde{h}(q)}$ in a Fourier representation as a real-space fractional derivative. Exploiting this analogy allows the modeler to map the effect of roughness on a spatially heterogeneous work of adhesion [42], which mitigates the need for fine spatial discretization in situations where most of the contact is conforming with the rough topography.

### 3.1.1. Hard adhesive contacts

The introductory discussion of this section revealed that even hard matter can be sticky when the particles are small or the interface smooth. However, it appears that sub-micron rubble particles constituting asteroids are still too hard to allow for this exception. Persson and Biele [131] realized that the observation of large asteroids spinning with a period of at least $T_{\text{min}} = 2.3$ h implies that adhesion forces are no significant addition to the (already pathetically small) gravitational forces in asteroids. $T_{\text{min}}$ is exactly the period at which the centrifugal and gravitational forces of an asteroid with radius $R > 300$ m and a small gravel particle balance each other out, assuming the mass density is similar to that of fused silica, i.e., $\rho \approx 2$ g/cm$^3$. Their analysis of the competition between surface energy and roughness led them to conclude that contact between rubble particles forms merely at one to three points, the latter number following from the requirement that
two particles squeezed together, even with a microscopically small force, must be torque-free. Similar considerations should apply to other granular matter like that used for 3D metal printing, where adhesion between granules might be best done by renormalizing the local radii of curvature rather than the surface energy. Wang and Müser [132] supported the findings by Persson and Biele with simulations showing that rough, adhesive contacts touch only in a single point during detachment when $\gamma/v_{ela} \lesssim 1/2$. (The torque argument does not apply because of periodic boundary conditions.) In separate work [106], they also provided a simple formula to estimate the scale- or size-dependent $v_{ela}$ from the expected height variance $\langle\{h(0) - h(\Delta\mathbf{r})\}^2\rangle$.

Thimons et al. [133] revealed that neither elastic nor rigid models could reproduce pull-off force measurements of spherical ruby probes retracted from ultranano- to microcrystalline diamond coatings. The contacts had been subjected to preloads leading to moderate nominal Hertz stresses of 135 MPa, but local contact pressures significantly exceeded common values for penetration hardness. After accounting for plastic deformation with a simple bearing-area model, measurements could be fitted assuming a rigid retraction and a generic cohesive-zone model, thereby highlighting the need to consider plasticity in (pre-loaded, i.e., non-asteroidal) hard-matter contacts. However, the deduced range of adhesive interaction of 5.6 nm leaves room for interpretation as to whether a more advanced elasto-plastic model is required for a truly quantitative description, or if capillary or electrostatic forces contributed to the adhesion. The fact that this range is compatible with other estimates employing different experimental methodology, e.g. atomic-force microscopy [134, 135], hints that there may be a fundamental gap in our knowledge about the forces acting across contacting interfaces.

Non-elastic effects can occur in small hard-matter contacts even without preload, as demonstrated in a combined experimental-numerical study by Vishnubhotla et al. [120]. They found adhesion-induced, reversible dislocation motion in self-mated platinum contacts, which caused significant hysteresis between loading and unloading as well as large deviations from expected linear contact mechanics, including the theory by Johnson, Kendall, and Roberts (JKR) [108]. Unfortunately, the extent of reversibility remains unclear, as the authors did not provide approach-retraction curves, and any residual dislocations may leave surface traces too subtle for clear experimental detection. Their expectation that these results have implications for probe-based microscopy and lithography is certainly valid for other solids that, like platinum, exert intense adhesive stresses on each other due to lacking passivation layers. However, the small adhesive stresses caused by the Van der Waals forces of passivated surfaces are small fractions of the Young's moduli of stiff solids and thus unlikely to cause significant deviations from (scale-free) elasticity. Nonetheless, an interesting aspect of the observed hysteretic dislocation motion is its loose analogy to (thermally activated) conformational changes in polymers, which are a primary cause of viscoelasticity and thus contact hysteresis.

### 3.1.2. Smooth contacts: continuously moving contact lines

Processes similar to those in ideal opening cracks (brittle fracture resulting in smooth surfaces) occur when indenters retract from a substrate, as well as at the trailing edges of rolling or sliding contacts, such as that represented in Fig. 4(a). Likewise, leading edges and approaching indenters can be seen as analogous to closing cracks. This is why traditional approaches to crack dynamics and their results pertain to the contact mechanics of bodies with deterministic or smooth surface profiles. In such contacts, continuous approach, retraction, sliding, or rolling imply continuously moving contact lines, which mark the boundary between closed (intact) and open (cracked) contact.

Closing contacts and opening them again generally entail hysteresis, as can be most easily seen from the following consideration. When an indenter is swiftly pushed into an elastomer, the latter will respond predominantly with its high-frequency modulus, so that not much contact is formed during the brief time of action. Keeping the indenter at its position for a long time then allows the elastomer to adapt to the shape of the indenter with its low-frequency modulus, leading to larger contacts. When the contact is broken again under swift motion, the elastomer will again respond with its high-frequency modulus, but now the contact is larger (and more relaxed!) than immediately after contact formation. The frequency-dependence of elastic moduli is also important, because a lower contact modulus leads to a steeper adhesive neck (as in the JKR limit of short-ranged adhesive, Hertzian indenters) than a larger modulus. Steep moving necks in turn imply large strain-gradient velocities and thus high dissipation. Effects of this kind have certainly
been understood for a long time; however, theoretical treatments like that by Persson and Brener [136], who elegantly exploited the observation that normal stress decays with the inverse square-root distance from a contact line but is bounded at small distances by the maximum interfacial tension, are involved and do not generalize easily beyond simple contact line motion. Numerical simulations are also challenging, because steep adhesive necks require high spatial resolution near the contact lines, forcing the simulator to compromise, e.g., to abandon efficient Fourier-based techniques, to have high resolution everywhere, or to focus on special geometries like rotationally symmetric [41] or singly connected contacts [42]. Alternatively, extending the range of adhesive interaction relative to the actual value while increasing the relaxation times of the elastomers accordingly, the correct resistance of a contact line to motion can be obtained also with coarse spatial resolution [137]. Qualitatively similar considerations as those mentioned above in this paragraph hold as well when the contact hysteresis originates from mechanisms besides viscoelasticity, such as crazing, cavitation, or plasticity. However, their quantitative description through meaningful theories is even more challenging and simulations typically require specialized modeling approaches beyond boundary-element methods.

Although crack-propagation theories are established, verifying their accuracy with rigorous numerical simulations that leave no room for tweaking parameters could only be achieved recently after a one-dimensional, short-range viscoelastic contact was discretized into $2^{18}$ surface elements [138]. Despite this fine discretization allowing for rather short-ranged adhesive interactions, the adhesion had to be labeled long-ranged with respect to the high-frequency modulus, which is why the detachment mechanisms transitioned from crack propagation to quasi-uniform bond breaking with increasingly large pull-off speed $v_p$. This occurred in agreement with the Persson-Brener model of uniform bond breaking at small scales [136], in which the product of contact area and maximum adhesive interfacial tension provides an upper bound for the pull-off force. The simulations [138] furthermore revealed local maxima in the viscoelastic work of adhesion $W_{\text{adh}}$ at intermediate $v_p$, which prompted the speculative suggestion that the enhanced energy supply at intermediate $v_p$ may contribute to increased pain felt when pulling off bandages at intermediate velocities. The relative maximum in $W_{\text{adh}}(v_p)$ is most easily rationalized for flat punches, for which the quasi-static work of adhesion $W_{\text{adh}} = 2\pi a^2 \gamma$ does not depend on $E^*$ so that this value also applies when the retraction speed is very high (given that the adhesion is still short-ranged). Viscoelastic dissipation can only add to the adhesion in between these two limits. Recently, Mandriota et al. [139] have observed a similar non-monotonic behavior of the effective work of adhesion in retracting axisymmetric contacts by tackling the problem with a specific numerical technique that can handle the limit of infinitely short range adhesive contacts, hence avoiding a very fine discretization of the contact. In the same study, they also confirmed that very high pull-off force can be achieved during fast retraction, at almost constant contact area (flat punch-like behavior) [140, 138]. This is similar in spirit to earlier work on (non-viscoelastic) adhesive contacts by Pohrt and Popov [141].

An open issue that has been recently addressed by simulations relates to an apparent “size-dependent” effect observed in some experiments of viscoelastic contact breaking, where the pull-off load was found to depend on maximum preload and thus on initial contact radius, e.g. [142], while other experiments, e.g [143] showed size-independence. Simulations showed that such an effective size-dependence emerges whether or not the viscoelastic solid was in its fully relaxed state or still in a transient [144, 145, 139]. The size dependence is hence really a history dependence of the contact, where an effective viscoelastic stiffening of the solid can promote cohesive breaking over edge propagation.

Interfacial contacts must form before they can break. The related process of crack closure has traditionally been studied by exploiting the $r^{-1/2}$ singular form of the stress near the crack tip [118, 119]. Such an approach suggests that viscoelastic losses increase the effective adhesion energy $G(v)$ of an opening contact by a similar factor as they reduce it for a closing contact at a given contact-line velocity $v$. Persson [146] proposed an energy-based method and arrived at similar $G(v)$ dependencies as traditional methods [118, 119], albeit without accounting for the stress field close to the crack tip. He argued that significant viscoelastic losses during crack closure must also occur near the crack tip in closing cracks, where traditional approaches would assume the material to be in a loss-free glassy state. This issue is now discussed controversially [110, 147, 145]. Resolving it could benefit from simulations in which contacts with
short-range adhesive even at high frequencies, while simulation cells are large enough to avoid finite-size effects.

In sliding and rolling motion, e.g., of a rigid object over a deformable body, contact closure and opening occur simultaneously. This unavoidably leads to retardation effects and thus dissipation even if there is no explicit interfacial friction and be it only because any body is made up of atoms with inertia, whose effect is abstracted away when postulating a high-frequency modulus. As can be seen in Fig. 4(a), the leading edge of a sliding tip has a smaller slope and thus less dissipation than the trailing edge. While it is easy to grasp that viscoelasticity and adhesion boost frictional losses in moving contacts, quantifying them is far from trivial [148, 149].

Recently, Carbone, Mandriota, and Menga proposed a new approach to moving, viscoelastic contacts [110, 150], which allows one to identify the contact area using energy balance principles. A peculiar prediction of the treatment is that the pull-off force assumes its maximum at intermediate velocities, as also observed in approach/retraction viscoelastic JKR-like simulations [139] and in qualitative agreement with the experiments carried out by Lakhera et al. [151] (see Fig. 4 therein). In Refs. [110, 150], the authors also show that the effective work of adhesion at the contact leading and trailing edges (closing and opening cracks, respectively) presents a non-monotonic trend with the sliding velocity (i.e., the excitation frequency), as previously suggested by Persson for opening cracks [152] in finite sized solids, thus entailing an additional friction peak related to adhesion hysteresis as experimentally shown by Grosch [153] and Roberts [154]. In fact, in the presence of a characteristic length *l* (either the contact length or periodicity), the frequency at which the material is excited far from the contact edge can be estimated as *v*/*l*, thus leading to a glassy (stiff-elastic) behavior at large velocity with no dissipation. The same authors also investigated the viscoelastic friction, which at relatively low (Fig. 4a top) values of *v*/*l* is mainly governed by viscoelastic induced adhesive hysteresis at the contact edges, whilst at relatively high (Fig. 4a bottom) values of *v*/*l* friction mainly depends on the viscoelastic bulk dissipation. Importantly, at the intermediate range of *v*/*l* both contributions take place, but cannot be linearly separated.

### 3.1.3. Wavy contacts: discontinuously moving cracks

When waviness is added to an otherwise smooth surface profile of an indenter, which is gradually retracted from or pushed into an adhesive, elastic body, the motion of the contact line will become modulated to the extend that the motion decomposes into stable and unstable segments, given sufficient waviness [116, 155]. Moments of contact-line arrest of such a serrated motion occurring in a singly-connected contact are depicted in Fig. 4(b). Arguing that each unstable segment dissipates energy (even during quasi-static operation, i.e., at velocities low enough where viscoelastic losses to be irrelevant), Guduru [116] proposed that waviness-induced, discontinuous contact line motion increases the work of adhesion. At the same time, the energy gained on approach is reduced, resulting in contact hysteresis similar in nature to that caused by viscoelasticity on smooth interfaces. It may be useful to note that the energy dissipated during an unstable segment can also be due to viscoelasticity. The important distinction is that the instability occurs on time scales much shorter than the external motion of the slider, in which case the detailed dissipation channel does not matter. Having two potential scenarios for energy dissipation has sparked some debate in the literature over the (predominant) origin of adhesive hysteresis in specific real or numerical experiments.

A seemingly strong argument in favor of viscoelasticity is the frequent observation of the Gent-Schultz law [156], which states that changes in the pull-off force scale with a (small) power law in pull-off speed $v_p$. Such a power law is a typical outcome of a smoothly moving contact line obtained readily with a standard linear solid. However, recent analysis of the Prandtl model [114, 115], where friction arises due to instabilities, reveals that a power-law response is also a natural outcome of instability-induced dissipation at finite temperature, with a logarithmic or Eyring-like dependence appearing only as a limiting case of weak pulling springs and low temperatures [157]. This suggests that care may be needed to pinpoint the origins of adhesive hysteresis.

In a theoretical treatment of the detachment process of a wavy sphere from a viscoelastic foundation, Ciavarella and Papangelo [158] suggest that the roughness-induced enhancement of the pull-off force identified by Guduru disappears at speeds near those where pure viscoelastic enhancement is comparable to that of elastic instabilities. This result aligns with expectations from the Prandtl model, assuming that explicit
damping occurs in the substrate rather than in the spring pulling the atom (or another generalized degree of freedom) over past the corrugated potential representing the substrate. The rationale in both cases is that damping smooths out instabilities. The damping found by Pérez-Ràfols et al. [159] for an elastomer with $E_∞/Ε_0 = 10$, did however not affect significantly the amount of dissipation, as the change in area was not unstable but still occurred at nearly constant load. Therefore, the pull-off load enhancement due to viscoelasticity and roughness was found to be roughly additive.

Recent work by Sanner et al. [43] correlates the stick-slip motion of contact lines observed experimentally (PDMS against nano-crystalline diamond) with that deduced from simulations in which the effect of roughness was mapped onto chemical heterogeneity. Fig. 4(b) compares the position of simulated and experimentally observed contact lines during moments of rest. The limit of short-range adhesion could be explored in the simulations through the use of a crack-front model for singly connected contacts [42], which only necessitates the contact line to be discretized rather than the entire surface. It remains to be seen to what extent previous studies may need reinterpretation, particularly those that matched only force-displacement curves without contrasting microscopic contact features.

Recently, numerical simulations have demonstrated that in the case of long-range adhesion, cavitation [160, 161], i.e., the rupture of a singly connected contact into two or more contacts separated by gaps, at a constant 'apparent contact size' [161], takes the role the elastic instability. Cavitation clearly complicates the contact problem, as now instead of a couple of advancing cracks there are multiple close-by crack fronts that require fine discretization, possibly interacting at small scales.

### 3.1.4. Rough contacts: jump-in and snap-off

Nominally flat contacts form at isolated sites, as in Fig. 4(c), with microscopic contacts forming and breaking during both normal and lateral motion. When the surfaces attract each other, asperities jump into contact and snap off again during sliding, as in cyclic normal loading. Tomlinson's suggestion that this process could be the main source of solid friction was reformulated over the years, notably by Fuller and Roberts [162], who replaced Tomlinson's atomistic model with a bearing-area approach and local JKR characteristics. Yoshizawa et al. [163] demonstrated that the friction coefficient of surfactant monolayer-coated surfaces closely correlates with adhesion hysteresis rather than with adhesion force or energy, as load, temperature, humidity, velocity, and holding times were changed. Their observation suggests that there are at least some systems in which adhesion hysteresis is the primary type of instability causing solid friction. Yet, numerical studies quantifying the contribution of adhesion hysteresis to friction in specific systems appear to be lacking.

The difficulty in modeling asperity jump-in and snap-off lies in the poor convergence of the energy hysteresis with decreasing adhesion range. Although the load-displacement relation of a receding adhesive Hertzian indenter closely corresponds to the JKR limit already at a Tabor parameter of $\mu_T = 5$, the indenter jumps into contact so prematurely relative to JKR that the adhesion hysteresis does not even reach 50% of the JKR limit at $\mu_T = 5$ [164]. Here, $\mu_T \equiv (R_c\Delta\gamma^2/E^{*2})^{1/3}/\zeta$ is a dimensionless measure for the (inverse) range of adhesion, where $R_c$ is the (local) surface curvature, $\Delta\gamma$ the work of adhesion, and $\zeta$ the actual range of adhesion. At $\mu_T \approx 100$, errors still exceed 10%, while adhesive hysteresis vanishes entirely once $\mu_T$ drops below $\mu_T \approx 1$. Thus, reproducing the short-range limit of adhesive hysteresis requires using small $\zeta$, which in turn necessitates a linear mesh size $\Delta a$ on the order of $\zeta^2$ to avoid lattice-trapping artifacts—i.e., jump-in and snap-off instabilities of individual discretization points, which rapidly increase energy dissipation [165]. (The $\Delta a \propto \zeta^2$ requirement arises from the stability that the curvature of the driving spring / local elasticity, which is proportional to $E^*/\Delta a$ must exceed the curvature of the interfacial interaction, which is proportional to $\gamma/\zeta^2$.) Since the error in adhesive hysteresis only diminishes with $\Delta a^{1/3}$, reducing it by a factor of 10 for a two-dimensional contact requires at least $10^6$ times the computational effort with standard methods. This highlights the need for new approaches, potentially beyond locally adaptive meshes, which are challenging to integrate with efficient Fourier-based boundary-element techniques.

As a numerical example, it may be instructive to see that adhesive hysteresis for a contact characterized by $R_c = 50$ nm, $\Delta\gamma = 50$ mJ/m$^2$ (typical for chemically passivated surfaces), and $\zeta = 3$ Å (typical for van der Waals interactions) yields a threshold value of $E^* \approx \Delta\gamma\sqrt{R_c/\zeta^3} \rightarrow 2$ GPa below which adhesive
hysteresis starts to kick in. This is of similar order of magnitude as the high-frequency modulus of many rubbers and elastomers.

To assess the origin of adhesive hysteresis, Müller *et al.* [111] successfully contrasted experimental and computed the approach and retraction process of flat punches to which single sinusoidal roughness was added. One of the studied topographies is depicted in Fig. 4(c). Not only the contact evolution, see Fig. 4(c), showed good agreement, but also force-displacement were found to be in good agreement for both surface patterns (triangular and hexagonal), after the viscoelastic nature of PDMS had been considered. In their system, the largest contribution to contact hysteresis was due to contact instabilities of minima and saddlepoints rather than those of the height peaks. For small preload and small operating velocities, these instabilities occured under compression but moved to tensile loads for large (compressive) preloads and “large” speeds of 25 µm/s. Near the moment of maximum tensile force, all height maxima in contact were under compressive load, as visualized in Fig. 4(c); tension originated exclusively from the saddlepoints and the ridges between height peaks. Pull-off forces were miniscule in comparison when only the height maxima formed contact during the moment of detachment. It thus seems that contact instabilities are not necessarily related to the highest asperities, as assumed historically, but due to the valleys in between the highest peaks.

## 3.2. Dry friction

### 3.2.1. Displacement coupling

The majority of works on elastic contacts focuses on purely normal loading in the absence of interfacial shear stress. However, a normal point force acting on a surface induces not only normal but also in-plane displacements just like an in-plane point force causes normal displacements in addition to in-plane displacements. This coupling, which Bentall and Johnson [166] emphasized to occur already in linear elasticity without exploring it much further at that point, arises from the resistance of an elastic body to density changes that are induced by normal or shear stresses acting on it. There are many consequences of coupling, both qualitative and quantitative, which have been explored mostly in single-asperity contacts. As evidenced by the lack of its discussion in the precedent community review on modeling in tribology [1], the importance of coupling is often underappreciated so that summarizing not only recent but also past achievements on this topic is certainly in place.

The linear coupling of in-plane and normal surface displacements of semi-infinite (!) solids does not need to be taken into account explicitly whenever the second Dundurs constant $\beta$ [167] vanishes. It is a measure for the relative difference of the contact compliance of two contact pairs, specifically, $\beta = \frac{(1-2\nu_1)/E_1 - (1-2\nu_2)/E_2}{2(1-\nu_1)/E_1 + 2(1-\nu_2)/E_2}$, where $E_{1,2}$ are the elastic moduli and $\nu_{1,2}$ the Poisson's ratios of the contacting materials. The constant vanishes when the solids are either elastically similar or both incompressible. This is because the displacements in response to in-plane and normal shear superimpose in that case. However, once $\beta$ is finite or a compliant solid has finite height, the superposition ceases to hold, because relative in-plane (normal) displacements arising due to normal (in-plane) stresses are no longer symmetry-related. In other words, an interfacial normal stress makes the two materials “want” to displace differently in the tangential direction. The subsequent motion or its suppression is counteracted by a frictional force, which would be absent in a contact in which both bodies have identical elastic properties and heights. Likewise, an interfacial tangential stress makes the dissimilar solids want to displace differently in the normal direction, which causes the normal stress or the interfacial separation to change compared to a frictionless contact, for which the effects induced by coupling are effectively taken care of by the contact modulus. In the following, (linear) coupling effects arising due to a finite second Dundurs parameter and due to finite or mismatching heights will be referred to as material and geometric coupling, respectively.

In one of the first quantitative studies on material coupling, Spence [168] found it to reduce the contact area in Hertzian frictional contacts at a given normal load compared to the classical, i.e., frictionless solution and to induce partial slip near the contact edge. This result implies that fretting can occur merely due to normal motion. To estimate the synergetic effects on fretting fatigue occurring due to combined normal-tangential loading, Nowell *et al.* [169] derived a numerical procedure accounting for coupling in contacts of Hertzian geometry. It allowed them to conclude that symmetry is broken in a way that coupling causes normal and in-plane stress to increase toward the trailing edge under gross slip conditions. These results
were supported by a semi-analytical approach [170], which later also indicated material coupling to decrease the critical tangential load needed to trigger gross sliding [171]. Recently, material coupling was also studied in sliding, rough contacts [172] and shown to change the quasi-linear relation of tangential force and the ratio of unstuck to total contact to one with a significant fraction of contact in sliding motion even at extremely small shear forces.

Since material interactions across an internal surface in a body follow the same rules as those across a “true” interface, coupling effects in layered [173] and graded [174] elastic contacts are intimately linked to those just described. For example, a pure normal impingement of a layered material can induce delamination inside the material.

In one of the first quantitative studies on geometric coupling, Bentall and Johnson [175] investigated a thin elastic strip passing through the nip of two identical and ideally counter-spinning wheels, one time as a free strip as in rolling (metalworking) and one time attached to one of the wheels like a tyre, which is also frequently encountered in processing machinery. Their study reveals that geometric coupling almost inevitably results in microslip at the extreme edges of the contact, particularly when the rollers are much stiffer than the processed strip. Nowell and Hills [176] overcame the restriction of ideally counter-spinning wheels by including the effect of the rollers' torque. In such an asymmetric set-up, a net interfacial tangential force is exerted on the film's surface. Ultimately, this increases the maximum normal pressures while skewing the normal stress to the contact's leading edge, i.e., to the opposite edge as under material coupling. Scheibert et al. [177] confirmed these trends through rigorous numerical models. They also tested the degree with which typical model assumptions—like linear elasticity, small-slope approximation, and Amontons's law at small scales—stand the comparison to an experiment, in which a glass cylinder indented a thin PDMS film resting on rigid substrate containing a MEMS-based device sensing normal and tangential stress fields. Predicted and measured stress fields were found to be within a few per cent.

Recently, Menga and various co-authors extended the analysis of coupling from partial slip single-asperity contacts first to gross slip of randomly rough line contacts [178, 179] and later to areal random roughness [180]. While some features are qualitatively similar as in single-asperity contacts, such as the skewing of stress distributions toward the leading or the trailing edge of individual contact patches depending on the (predominant) nature of coupling, new features arise. This includes the observation that coupling can destroy the validity of Amontons' law when it holds microscopically in a way that geometric coupling increases interfacial and viscoelastic friction, while material coupling has the opposite effect. Moreover, both geometric and material coupling reduce the propensity of interfacial fluid flow or leakage especially when sliding occurs normal to the pressure gradient (see Fig. 6a₂-6a₃), which, in parts can be related to the coupling-induced asymmetries of individual contact patches. Coupling moves the maximum of the in-plane tensile stress to the trailing edge of a contact and, potentially more importantly, increases its value substantially (see Fig. 6a₁). Thereby crack nucleation and propagation occur more easily, which can boost the risk of material failure.

Quite a few questions related to coupling have remained unexplored. For example, it has not yet been investigated how coupling affects the normal contact stiffness in areal contacts, which is the change of normal pressure with the mean interfacial separation (though, line contacts results suggest coupling-induced lower contact stiffness [178]). It is an important property, because it sets an upper bound for interfacial electrical conductance, e.g., in sliding electrodes [184]. An extra resistance comes from a thin oxide surface on metals and dendrites growth, which is inversely proportional to the true contact area, whose size and shape is again affected by coupling [185, 186]. At the same time, interfacial stiffness sets a lower bound for the heat conductance, as evanescent waves allow for heat conduction in addition to heat diffusion. Their relevance depends on the distribution of the gap, which can also be affected by coupling [187, 188]. Surprisingly unexplored is also the effect of coupling in (multi-) layered configurations, such as human skin and road tyres, even if Wang et al. [173] set the framework for such a study. It also seems unexplored how coupling affects adhesion in particular under tangential loading. However, even without coupling, it remains challenging to predict how frictional shearing affects the contact under a normal and/or adhesive force, as is discussed next.
<figure><img src="image_7.png" /><figcaption>Figure 6: (a) confinement induced linear coupling effects on in-plane maximum principal stress (a₁) and leakage anisotropy (a₂,a₃). Adapted from [180]. (b) Images of PDMS lens vs smooth glass contact in static (b₁) and steady gross slip (b₂) conditions, showing anisotropic contact shrinkage. Adapted from [181]. (c) Measured longitudinal in-plane stretch (c₁) and derived shear stress (c₂) distributions in shrunk sliding contacts. Adapted from [182]. (d) Measured sliding contact area (red) compared to static one (dashed) (d₁) and analogous nonlinear FEM simulations results for Neo-Hookean (d₂) and Mooney-Rivlin (d₃) models highlighting the effect of finite deformations/displacements in contact shrinkage (green area). Adapted from [183].</figcaption></figure>

### 3.2.2. Contact shrinkage

Even in the absence of elastic coupling between normal and tangential displacement at the interface, combined tangential and normal loading of mechanical contacts may induce changes in the contact size, as often observed in experimentally [181, 189]. Despite significant efforts made over several decades, the origin of this is still discussed, which is why a summary of past and very recent achievements in the field may be in place. Before cutting to the chase, we wish to specify that we are not concerned with the friction arising due to the viscoelastic dissipation that a sliding (stiff) indenter induces in an elastomer in the absence of explicit interfacial friction. That issue is dealt with in Sect. 3.1 as the attentive reader certainly remembers. Here, we discuss explicit, interfacial friction, which is often modeled using a local constitutive equation of the form $\tau = \tau_0 + \alpha p$, where $\tau$ is the shear stress, $\tau_0$ an offset, $\alpha$ a dimensionless parameter, and $p$ the normal compressive stress, e.g., that induced by adhesive attraction, squeezing the surfaces against each other.

Savkoor and Briggs [190] were the first ones to observe a slight reduction of contact area when tangential loading is combined with a squeezing normal force in macroscopic contacts. Since the contact shape stayed axisymmetric during the application of the tangential load, they related the observation to adhesive contact mechanics. Specifically, they postulated friction to allow for stuck interfacial conditions (leading to singular shear stresses at the contact periphery and to mode II and III brittle crack propagation) and argued that the tangential fields increase the elastic energy stored in the body and, in turn, the energy release rate, which is the reversible energy per unit contact area variation available at the interface to break the adhesive bonds. Since the work of adhesion per unit area $\gamma$ stays constant, this requires the contact area to reduce at equilibrium. However, their prediction overestimated the contact radius reduction. Inspired by works on fracture mechanics by Hutchinson [191, 192] while aiming at empirically correcting Savkoor and Briggs explanations, Johnson [193, 194, 195] argued that, for a contact shrinking under the actions of normal and tangential loads, only a portion of the released elastic energy is available for reversible transformations (effective energy release rate). In fact, Johnson's theory focused on the role of partial slip conditions, so that a stuck inner circle (with bounded shear stresses) is surrounded by an annular region where slip occurs under constant shear stresses, i.e., $\tau = \tau_0$ - in contrast to Cattaneo [196] and Mindlin [197] who used $\tau_0 = 0$.
Since in-plane displacements at the contact interface also depend on the contact size [198], Johnson assumed that part of the elastic energy released during the contact shrinkage is dissipated in the annulus by frictional relative sliding. However, the rationale behind this picture was only vaguely identified in possible surface heterogeneity able to reversibly store energy (other than the Dupré work of adhesion) through 'spatially high frequency slips' and 'evolution of the interface-dislocation structure at the atomic level' [195], or by 'injecting surface microstructures into the material due to locally heterogeneous deformations' [199]. Only recently, as discussed below, Barras et al. [200, 201] have tried to quantify the effective energy release rate, though in the much more complicated framework of nonlinear elastodynamic FEM calculations. On the contrary, the phenomenological Johnson's approach exploited 'rather arbitrary' [194] functions and experimentally fitted coefficients to describe the 'mixed-mode' effective energy release rate. As a consequence, the model had limited predictive utility, being potentially able to either recover the Savkoor and Brigg's overestimation of the contact area reduction or to indicate no reduction at all.

Nonetheless, in some cases (mostly at the small scale) contact area changes were experimentally found to be marginal to non-existent, as shown by Israelachvili's group [202] for atomically smooth mica cross-cylinder contacts, by Carpick et al. [203] for AFM tips coated with mica, and by Vorvolakos and Chaudhury [204] even for polymeric contacts up to sliding velocity of order $10^{-3}$ m/s (with shear stresses amounting to almost 10% of the elastic modulus). The independence of the adhesive contact area on the tangential loading was also analytically predicted, in linear elasticity, by Menga et al. [205, 206] in the case of a rigid sphere sliding past an elastic solid in the presence of uniform interfacial shear stress [198]. In the same work, they show that shear stress spatial fluctuations may lead to a reduction in the contact area, as the energy release rate increases globally. Importantly, detailed in-situ observations of the contact interface rarely report a contact area increase during sliding, even when accounting for factors such as surface roughness [207], which enhances viscoelastic adhesion [208], and elasto-plasticity in spherical contacts [209].

Later experimental observations on basically smooth contacts found more significant anisotropic effects, mostly shrinkage along the sliding direction (the trailing edge has a greater radius of curvature than the leading edge), and potentially a slight growth in the orthogonal, in-plane direction, as shown in Fig. 6b. More importantly, the contact size was also reported to drop below the Hertzian value [204, 189, 181], suggesting that a proper model should also account for nonlinear elastic effects. This was also confirmed by Chateauminois's observations in polymeric sliding contacts [210, 182], reporting strain values localized at the contact edges well beyond the linear elasticity threshold (see Fig. 6c). Indeed, FE could reproduce contact behavior including finite deformations and rotations. Recently, these have been investigated in detail by Stupkiewicz's group [183] showing, as reported in Fig. 6d, that the anisotropic contact shrinking observed in Ref. [189] is a non-adhesive phenomenon mostly ascribable to 'local contact lifting' occurring in the contact region and induced by nonlinear coupling in neo-Hookean elasticity. On similar grounds, experiments and nonlinear FE calculations showed that the contact area can even slightly increase in the partial slip regime before gross slip occurs due to the enlargement of the leading edge induced by the local compression and consequent material bulging [211]. Recent FE calculations also confirmed the role of finite elasticity in frictional contacts, even in purely normal indentation, showing that hysteresis can be observed in nominally uncoupled contacts and partial slips nucleate somewhere in the contact region (depending on the local indenter slope) rather than at the contact edges [212].

Also, the classic assumption of quasi-static contact evolution lacks physical meaning. Since the frictional rupture is a path-dependent process, the residual stress and displacements at the interface, and the energy dissipation, and, in turn, the effective energy release rate associated to a contact area variation crucially depends on elastodynamic bulk effects (e.g., wave radiation from the frictional interface to the surrounding material) [200, 201]. Moreover, a series of studies by Fineberg's group [213, 214, 215] have shown that macroscopically observable quantities, such as the contact area and the frictional sliding stability, in rough sliding contacts suffering from confinement coupling are strongly affected by the interaction between the high-speed elasto-dynamic propagation of spontaneous rupture fronts, the degree of confinement, and a postulated velocity-dependent friction law. Despite many promising works dealing with the question of how combined tangential and normal load affects contacts and their mechanics, it remains challenging to quantify the synergistic effects in real tribological systems.
### 3.3. Lubrication

Lubrication is the process of reducing friction and wear in mechanical contacts by introducing a substance—such as oil, grease, or even water—that forms a protective film between two surfaces. The sliding motion between these surfaces generates hydrodynamic lift, which facilitates film formation and thus separates the sliding bodies. Lubricant research in modeling encompasses diverse aspects, most importantly the theoretical or numerical modeling of fluid flow at moving interfaces, the engineering of surface textures, but also the molecular simulation of bulk liquids to aid the modeling-supported design of new lubricants. Key challenges are to properly account for the interaction between fluid flow and solid mechanics, but also to reach meaningful shear rates in molecular simulations as well as sufficient accuracy as it comes to making chemistry-specific predictions [216]. Regarding this latter point, we note exemplarily that effectively modeling the chemical binding between friction-modifying additives—such as gallate molecules in aqueous lubricants—and surfaces is infeasible with density-functional methods, inaccurate even with advanced conventional force fields, but apparently satisfactory when employing machine-learned potentials [217].

#### 3.3.1. Bulk Lubricant Properties and Design

Lubricants are crucial in determining the performance and durability of tribological contacts. Designing lubricants is complex because base lubricants are typically combined with additive molecules that perform various functions, such as modifying viscosity or protecting sliding surfaces. In recent years, there has been a surge in molecular methods employed for designing bulk properties of lubricants, such as their viscosity. However, the proper design of surface-active species—such as for wear protection—has largely remained a trial-and-error process. Designing new lubricants begins with the ability to predict their constitutive behavior—such as viscosity, shear thinning, and how these properties change with temperature and confinement—from the molecular structure of the lubricant.

Discussions about the best functional forms as well as the physical origin of constitutive relations remain active. Jadhao and Robbins [10] observed that squalane is essentially fully aligned with the shear flow at intermediate shear rates $\dot{\gamma}$ when the effective viscosity $\eta_{\text{eff}} = \sigma/\dot{\gamma}$ has only dropped by a factor of three compared to the equilibrium viscosity $\eta_0$, in a broad range of temperatures $T$ and pressures $p$. At low $T$ and large $p$, the ratio $\eta_{\text{eff}}/\eta_0$ can become many orders of magnitude smaller than one third, although the molecules do not further align parallel to the flow with increasing shear rate. Thus, attributing the reduced effective viscosity to alignment is often premature and risks confusing correlation with causation. Viscosity and shear thinning appear to be determined by how energy barriers are overcome in a lubricant during enforced shear flow and most notably by the energy difference at the onset of a structural instability and the point, where a local structure has found a new stable position [157].

The conceptual pitfall of confusing correlation with causation may also appear in the context of pressure dependence and free-volume interpretations of the equilibrium viscosity $\eta_0(p)$. While most liquids certainly increase their steric repulsion with increasing pressure—whereby free-volume theories may appear quantitative—fluids with typical shear thinning can display negative piezoviscous coefficients, indicating a decrease in equilibrium viscosity with increasing pressure [218, 219]. For example, tetrahedral network formers such as water have an increased number of coordination defects with increasing pressure, providing a mechanism for how reducing (free) volume can enhance mobility, as the latter is carried by these high-density coordination defects [220]. Correlating lubricant molecular structure with viscosity appear to be particularly challenging when viscosity is time dependent as in viscoelastic fluids. This can increase load-bearing capacity and reduce friction [221]. Assessing the lubrication performance of viscoelastic fluids is, *inter alia*, crucial in polymer processing, where viscoelastic models may be applied to polymer melts and solutions [222], or biomechanics, where biological fluids, including e.g. the synovial liquid [223], are marked by a strong viscoelastic rheology.

Despite frequent challenges in correlating fluid viscosity with simple geometric parameters of lubricant molecules, there are certainly exceptions. For example, the insight that the hydrodynamic radius and the distance a molecule must be displaced before its motion becomes diffusive can be quite useful for estimating viscosities from short molecular runs [224]. Exploiting such correlations can benefit the simulation-assisted design of new lubricants, as equilibrium viscosities are notoriously difficult to determine from molecular simulations.
A paradigm example for the simulation-assisted design of base lubricants was presented by Kajita, Kinjo, and Nishi [5]. They used a “genetic” algorithm that randomly combines different fragments of lubricant molecules into new candidates to solve this design problem, where minimizing the change of viscosity η with temperature was the design target. At a given temperature, η was estimated from very short molecular dynamics simulations by linking short-time shear-stress autocorrelation functions to equilibrium viscosities within a modified shoving model [225]. Their work [5] is a first demonstrator that automatic lubricant design is possible, however, the design of complex additive combinations requires further research into accelerated simulation methods.

### 3.3.2. Thin-film flow and the percolation problem

Macroscopic models for lubricated tribological systems, as well as mesoscopic models that explicitly represent surface roughness, often begin with simplified hydrodynamic descriptions. In the classical Reynolds’ partial differential equation (PDE),

$$ \nabla \cdot \mathbf{J} = 0, \qquad (2) $$

J is the flow rate, flow is assumed to be of a Newtonian fluid, and stick-boundary conditions are considered, leading to

$$ \mathbf{J} = -\frac{h^3}{12\eta}\nabla p + \frac{h}{2}\mathbf{U}. \qquad (3) $$

Here, h denotes the interfacial separation, η the viscosity, U the entrainment velocity, and p the fluid pressure. Since it is generally unfeasible to resolve small-scale details of h when solving the Reynolds equation, a number of homogenization methods have been proposed.

One of the earliest homogenization approaches goes back to Bruggeman [226], whose work allows the leakage in seals to be homogenized, i.e., when the shear-induced flow rate or Couette flow $J_c = hU/2$ term is null. This limit matters not only theoretically but also concerns the efficiency of components in industrial applications [227]. For leakage problems, the flow rate reduces to $J = \bar{\sigma}\nabla p$, where $\bar{\sigma}$ is the effective fluid-flow conductance at coarse scales, accounting for the microscopic surface roughness. For isotropic roughness profiles, $\bar{\sigma}$ can be estimated from the microscopic conductance $\sigma = h^3/12\eta$ by solving the self-consistent equation

$$ \frac{1}{\bar{\sigma}} = \int \mathrm{Pr}(\sigma) \frac{D}{\sigma + \bar{\sigma}(D-1)} d\sigma, \qquad (4) $$

for $\bar{\sigma}$ [228]. Here, $\mathrm{Pr}(\sigma)$ is the probability density for the local conductivity $\sigma$ to occur, while D is the spatial dimension of the interface.

Eq. (4) is known to yield astoundingly accurate results even for highly complex contact geometries [229] even close to the percolation threshold [230, 231], in particular when D is replaced with an effective interfacial dimension that reproduces the correct percolation threshold for the relative contact area, e.g., $a_c^* \approx 0.42$ for D = 2 and random roughness, which implies $D_{eff} = 1.72$ to be used instead of D = 2 in Eq. (4) in that case – deterministic profiles can lead to different percolation thresholds [230]. This outcome is also fully consistent, in terms of $a_c^*$ with other contact mechanics results obtained on purely deterministic way and available in literature [232]. The accuracy of the theory is demonstrated in Fig. 7, which compares full Reynolds flow calculations conducted on the gap topography from the contact-mechanics simulation of a randomly rough surface to the results of a calculation based on Bruggeman’s effective-medium theory.

Persson [233] generalized Eq. (4) to anisotropic surface roughness that can be described via the Peklenik number, i.e., for surfaces whose topography results from isotropic roughness by compressing one in-plane direction relative to the other. Wang and Müser [234] later corrected a typo in the final equation and showed that the extension to anisotropic surfaces was accurate and that the percolation threshold of $a_c \approx 0.42$ for two-dimensional surfaces remained unaffected within numerical accuracy. A remarkable reason for investigating anisotropy is related to what happens in sliding viscoelastic contacts: here, even in presence of isotropic interfaces, the solution in terms of contact areas, displacements, and pressure distribution is anisotropic due the different relaxation mechanism between leading and trailing edges [235, 236]. A relation between percolation and anisotropy, as previously described, might also potentially arise also from coupling between normal and tangential tractions [180].
<figure><img src="image_8.png" /><figcaption>Figure 7: Fluid pressure during fluid flow through a contact formed by a randomly rough surface and a linearly elastic body (a–d), where fluid pressure decreases from blue to green. Panel (e) shows the mean current density *j* through the contact, where letters reference the previous panels. Panels (c) and (d) highlight a situation near the percolation threshold, where the fluid pressure drops almost entirely near a critical constriction. The inset of panel (e) supports the hypothesis that fluid leakage in an elastic contact disappears with the exact power law $(L - L^*)^{3.45}$ as a function of load *L* when the percolation load $L^*$ is approached [230]. Full lines depict the Bruggeman effective theory [228], while symbols show different disorder realizations. Adopted from Ref. [231].</figcaption></figure>

For sliding contacts, the full Reynolds thin-film equation needs to be regularized so that Couette flow is considered in addition to Poiseuille flow, potentially including the effect of anisotropic surface roughness. This can be achieved again with the Bruggemann effective medium [237] but also using perturbation [238] or renormalization-group based approaches [239], or an approach [240] based on Persson's contact mechanics theory [4]. It turns out that the shear stresses between the fluid and the solid surfaces can induce a non-negligible, anisotropic flattening of asperities belonging to initially isotropic surfaces [241]. This can greatly impact the magnitude of fluid-solid interactions, just like the orientation between anisotropy and sliding direction. Specifically, when an O-ring slides steadily over a grooved, randomly rough surface, mixed friction decreases when the grooves are oriented perpendicularly to the sliding direction. This reduced friction is attributed to the enhanced coupling between fluid and asperities, leading to increased fluid pressure and reduced asperity-asperity contacts. [240]. Fluid-induced asperity flattening might also cause the occurrence of a secondary (micro-)elastohydrodynamic friction regime in the Stribeck curve, as demonstrated in the case of a grooved rough soft contacts [242].

Another key problem is that constitutive behavior may no longer be known for highly confined and pressurized fluids, or that the constitutive behavior starts to explicitly depend on the gap separating the sliding interface, for example because of epitaxial ordering [243], fluid layering [244], or jamming of the fluid. This high-dimensional parameter space can be explored using active learning, e.g. by using Gaussian process regression to interpolate between molecular dynamics simulations and decide in which regions of the parameter space to run new calculations, see Sect. 2.3 and in particular Fig. 3 or Ref. [55].

### 3.3.3. *Surface texture and functionalization*

The Reynolds' equation works extremely well in the limit of hydrodynamic lubrication, where interfaces are well separated and no solid-solid contact occurs. This limit is achieved at high sliding velocities and with viscous fluids. Optimization of sliding systems often aims to extend this hydrodynamic regime to slower velocities and lower viscosities, the latter often introduced to increase energy efficiency. One way to achieve
<figure><img src="image_9.png" /><figcaption>Figure 8: (a) Parameteric [248] and (b) nonparametric [250] optimization of surface textures. The “adjoint” textures in panel (b) are obtained using PDE-constrained optimization employing the adjoint method. Adapted from Refs. [248] and [250].</figcaption></figure>

this optimization is to pattern sliding interfaces with dimples, an approach pioneered by Etsion [245]. In addition to increasing the load-bearing capacity, dimples are also considered beneficial for capturing wear particles, effectively removing them and their adverse effects from the contact. Furthermore, dimples directly participate to the friction reduction, promoting cavitation when the fluid expands in the micro-hole, and, thus, a remarkable reduction of the shear stresses [246, 247]. The latter are, in fact, much lower when the fluid lubricant is replaced by a gas and contribute to the overall friction reduction in the lubricated system.

Figure 8 shows examples of surface textures simulated and manufactured for pin-on-disc experiments. Putignano et al. [248, 249] used predetermined patterns of dimples (Fig. 8a) and concluded that the main geometric factor influencing performance is the fraction of the surface occupied by dimples. While optimizing dimple geometries is typically carried out in this parametric manner, that is, by designing dimple patterns and testing them in numerical calculations, Codrignani and co-workers [250] used PDE-constrained optimization based on the Reynolds equation, which does not impose predefined patterns, to optimize the geometry of a sliding contact. This optimization strategy yielded the geometries shown in Fig. 8b, where the step-like geometries on the left are those predicted through completely free-form optimizations to most efficiently increase load-bearing capacity. The dimple patterns in Fig. 8b were obtained by constraining the optimization to specific locations of the dimples. Codrignani [250] also points out that the ideal patterns depend on the operating conditions (load, sliding velocity, etc.) of the tribological contact.

### 3.3.4. *Soft lubricated contacts*

In Sect. 3.1, we discussed sliding on dry viscoelastic solids, which has been extensively studied using theoretical techniques over the last decade. Viscoelasticity of the contacting bodies also plays an important role when the contact is lubricated, as is often the case in soft lubrication and biomechanics. The scenario commonly found in elastohydrodynamics (EHL) changes dramatically in terms of both film thickness and pressure distributions [251, 252, 249]. This leads to significant variations in the classic Stribeck relation between fluid speed and friction coefficient [251]. In fact, friction is produced not only by fluid losses but also by hysteretic viscoelastic losses. As shown experimentally by Sadowski and Stupkiewicz [253], this means that the classic sequence of boundary, mixed, and elastohydrodynamic lubrication regimes is altered by the occurrence of a bell-shaped viscoelastic hysteretic friction contribution to the Stribeck curve. Consequently, Putignano has shown that when dealing with visco-elasto-hydrodynamic lubrication (VEHL), the Hersey number—related to the fluid entrainment velocity, fluid viscosity, and load—is not the only physical parameter governing the process. It must be complemented by the two dimensionless speeds of the two viscoelastic bodies, leading to three different and independent dimensionless groups that characterize the VEHL regime [252].

The above considerations make it clear that the time-dependent dissipative characteristics of a viscoelastic contact surface significantly influence the dynamics of the lubrication process, especially under non-steady-
state conditions, where transient effects [254] are amplified by solid viscoelasticity [255]. This is exemplified in the case of lubricated squeezing of viscoelastic bodies: before reaching a Hertzian-like pressure distribution, the system exhibits a distinctly non-elastic evolution, with viscoelastic-induced pressure peaks at the contact edges [256]. Additional effects arise from various configurations of viscoelastic solids, such as thin layers [257].

Approximate yet effective methods have been developed for multi-field problems, particularly for fully coupled thermo-viscoelastic contact problems, with applications to layered materials and coatings as well as the simulation of frictional and viscoelastic heating [258]. The main challenges related to soft lubrication—including the impact of material viscoelasticity, large deformations, and non-Newtonian effects in the fluid—are further discussed in [259].

### 3.3.5. Aqueous systems

The pursuit of sustainable and energy-efficient lubricants has sparked interest in water-based lubrication for technical applications. Much of the expertise in water lubrication comes from research in winter sports, which often involves contact with ice and aims to improve sports equipment performance. [260]. However, the mechanisms of ice friction remain a topic of considerable debate. Two common explanations are typically cited: first, that frictional heating melts the ice, creating a lubricating water film; and second, that due to the unique properties of water (shared with other tetrahedral liquids [261]), pressure lowers the melting point, also causing the ice to melt. Recent molecular dynamics simulations, however, have shown that ice can undergo shear melting even both at temperatures far below and relatively close to freezing temperatures [262], similar to the solid-state amorphization processes observed in other tetrahedral solids [263, 264].

Water in its liquid form is also the primary lubricant in biological systems. Joints and other body parts are lubricated by synovial fluid, which consists of water enriched with a variety of biomolecular 'additives' such as lubricin. Biological sliding interfaces are often composed of cartilage tissue, a hydrogel that requires hydration. Water is expelled from the hydrogel during compression, but Burris and co-workers [265] recently showed that the formation of a lubricant film during sliding helps rehydrate cartilage. Putignano and co-workers [266] used an elastohydrodynamic lubrication (EHL) model combined with a mean-field model for porous-media flow through cartilage to demonstrate that sliding induces a net flow through the cartilage, rehydrating the hydrogel. This mechanism appears to be essential for maintaining cartilage stability under severe frictional conditions and supports the advice from physiotherapists to keep joints in motion and, recently, has been artificially replicated on hydrogels [267].

## 4. Materials

### 4.1. Metals

Metals are typically lubricated in engineering applications, except in the case of sliding electrical contacts in which lubricants can affect conductivity. Nonetheless, extreme conditions can lead to the squeeze-out of lubricant and even the removal of the (protective) oxide layer. Modeling friction and wear in metals, particularly at the atomistic level, therefore often assumes that any oxide layer has worn away during sliding, leaving the tribological response dominated by cold-welded metal junctions, adhesion, and plowing. Moreover, soft metals, for example gold, are even used as solid lubricants in vacuum bearings for space applications [268, 269], which adds to necessity to study direct metal-metal contacts. The material's microstructure—such as grain boundaries and defects—and surface topography then influence friction, wear, and scratch resistance [270] of the metal contact.

#### 4.1.1. Friction and wear of pure metals

Plasticity in metals is governed by two dominant bulk deformation mechanisms, dislocation-mediated plasticity and grain boundary sliding, and these mechanisms are also responsible for the tribological response [271]. Dislocation-mediated plasticity dominates when grain sizes are above approximately 10 nm, as intragranular deformation—through the nucleation and motion of dislocations—becomes energetically more favorable than intergranular deformation for larger grains. However, in refined grain structures with sizes below about 10 nm, grain boundary sliding dominates and can lead to low friction. Molecular models
<figure><img src="image_10.png" /><figcaption>Figure 9: Schematic summary of molecular dynamics (MD) simulations and experimental observations for a gold–tungsten (Au–W) sliding contact. Panel (a) shows MD snapshots illustrating the formation of a loosely packed interlayer (IL) of Au (yellow) at the interface with W (gray) as a function of increasing W-Au distance. This IL is associated with a significant reduction in friction. Panel (b) presents a stylized schematic of the high-resolution transmission electron microscopy (HRTEM) setup used in the accompanying experiments. Panel (c) compares the frictional response in interfaces with and without an IL. Figure recreated and reinterpreted based on Ref. [274].</figcaption></figure>

for friction in pure metals connect the tribological response to the behavior of cold-welded commensurate (epitaxial) or incommensurate junctions.

Understanding the relationship between microstructure and friction in metals allows for tailoring frictional properties in practical applications for both bulk metals and surface films. For instance, recent experiments have demonstrated that a single-crystal nickel (Ni) surface with the 011⟨211⟩ orientation preferentially develops an ultra-nanocrystalline near-surface layer when subjected to shear against a Si₃N₄ ball [272]. This transformation leads to a reduction in the coefficient of friction from 0.64 to 0.32, resulting from a transition from dislocation-mediated plasticity to grain-boundary sliding. Another approach to forming ultrafine grains that lead to low friction via grain-boundary sliding is through the careful selection of layer thickness in gold-nickel (Au-Ni) multilayer systems [273]. In these systems, mechanical mixing under tribological conditions in ultra-high vacuum for layer thicknesses below approximately 50 nm led to an increase in grain boundary density, lowering friction by a factor of about two.

High-fidelity studies often combine atomistic modeling with high-resolution microscopy to directly validate simulation outcomes. As an example, He et al. [274] combined in-situ high-resolution transmission electron microscopy (HRTEM) and atomic force microscopy with molecular dynamics simulations to observe interfacial structure formation during frictional processes in atomic-scale Au-W contacts, see Fig. 9. This approach enabled real-time visualization of the interface microstructure, overcoming limitations of indirect or ex-situ techniques. They found that a loosely packed, disordered interfacial layer acted as a lubricant between tungsten (W) and gold (Au) asperities, resulting in friction coefficients below 0.01 in nanoscale contacts, similar to structural lubricity. The relative motion between the Au and W asperities was mediated by diffusion of Au atoms rather than by rigid shear, yielding, or junction fracture, as would occur in high-friction scenarios. The intimate connection between microstructure and friction explains why emerging modeling efforts are focusing on understanding aspects of wear, plastic deformation, and scratch resistance of metallic surfaces.

Wear can be distinguished by various mechanisms including abrasion, fatigue, adhesion. Two-body wear involves the direct interaction between two surfaces, resulting in material exchange and removal due to
<figure><img src="image_11.png" /><figcaption>Figure 10: Surface fragmentation and three-body formation in a discrete-particle simulation of sliding wear. The simulations are configured so that brittle fracture occurs at length scales smaller than the simulation box—something not achievable in atomic-scale models. Colors indicate distinct fragments that detached from both sliding surfaces and eventually agglomerated to form a third body. Cold welding between surfaces is prevented by dynamic detection of surface atoms that are assigned a weaker adhesion energy, mimicking the effects of oxidation and surface passivation in metals. From Ref. [278]</figcaption></figure>

abrasion and adhesion. Three-body wear involves abrasive particles caught between two surfaces, resulting in wear from the rolling and sliding actions of the particles. The material properties, shape, and size of these particles influence the degree and mechanism of three-body wear. Recent modeling advancements enable realistic simulations of mechanical wear, leading to a better understanding of wear processes at different scales. Molecular dynamics (MD) and discrete element methods (DEM) can describe wear mechanisms at small scales, while finite element (FEM) and boundary element methods (BEM), as well as phase-field models have been used to study subsurface stress evolution and damage in sliding contacts. In the following, we will highlight some of these studies.

While MD simulations excel at handling complexities associated with severe deformation at small scales, they struggle with phenomena at larger scales like crack shielding and amplification or debris formation. The latter has been recently simulated using coarse-grained potential based discrete element simulations [44, 46], that are similar in spirit to molecular dynamics but operate on coarse-grained elements that can deform and interact via forces and torques not necessarily derived from an underlying energy function. The calculations revealed formation of wear particles (Fig. 10) and the existence of a critical junction size between two colliding asperities where the failure mechanism upon sliding transitions from smoothing by plastic deformation to subsurface fracture and debris formation—a mechanism postulated more than 60 years ago by Rabinowicz [275]. These works are notable for being able to treat brittle propagation and interaction of cracks, leading to fracture. This is in contrast to molecular dynamics results, that in most cases show ductile behavior of the contacting bodies and blunting of pre-existing cracks. Indeed, the critical-size idea explains why this happens: The length scales accessible in atomistic calculations are much smaller than the size of the critical junction (or crack-tip process zone) expected in most materials, with the potential notable exception of brittle crystals [7, 276] or amorphous carbon [277]. Studying the effect of material heterogeneity on asperity failure, it was shown that the critical junction size formula works for hard materials until a ductility threshold is reached, leading to mixed behavior between shear localization and mode II crack opening [47]. This reveals that simple homogenization schemes are insufficient for capturing all underlying aspects of deformation, as local fluctuations increase ductility, but also create brittle areas prone to crack formation, resulting in competing mechanisms.

Follow-up works [279, 280] identified three distinct mechanisms for material removal during adhesive wear at the asperity level: atom-by-atom removal (atomic attrition) [281, 282, 283, 284] (see also Sec. 4.2.1), plasticity (dislocation-based or amorphization-induced material removal) [285, 286], and fracture-induced fragment detachment [287]. Mechanistic wear models have been recently developed based on the critical junction size model [41, 39], that connects the macroscopic wear rate, particle sizes, and emission rates. Simulating the formation process of third-body particles and the wear they induce on sliding surfaces presents numerous challenges. For example, Ref. [288] showed that detached fragments can first act as roller bearings while transitioning into a shear-band-like state via adhesive bonding, suggesting that proximity between particles and size relative to surface roughness affect the stability of the rolling regime. However, Aghababaei
demonstrated [278] that interfacial adhesion also affects the particle growth rate by controlling the particle
formation probability and determining the energy per unit volume of removed material (see also Fig. 10).

These studies emphasize the significant impact of wear particles on friction and wear, even at the
nanoscale, highlighting the need to incorporate additional complexities such as the coupling between surface
roughness and wear particle dynamics in continuum simulations. Recent work [289, 290] investigated the
influence of microscale features such as surface roughness on the wear behavior of engineering systems by
coupling semi-analytical models into Archard's model-based finite-element framework. The micro-contact
analysis shows a considerable decrease local pressure and temperature with time due to wear. The poten-
tial of multiscale approaches to model wear at the mesoscale has also been recently demonstrated with a
concurrent multiscale framework, in which a coarse-grained discrete approach [44] was coupled with finite
elements to simulate fracture-based material removal during adhesive wear [291]. A related multiscale effort
to study surface finishing based on MD and the material point method (MPM) [292, 293] found evidence
of material removal rates and surface roughness parameters scaling linearly with the abrasive size over four
orders of magnitude, allowing predictions regarding the effects of abrasive blunting on the surface quality of
the work piece.

Scratching is considered a subset of two-body wear, particularly when abrasion predominates the process of material removal. Large-scale MD simulations of single-crystalline W and Cu [294, 295] demonstrated that the connection between adhesion and abrasion in material removal is strong at shallow scratching depths but weakens as the scratching depth increases. In fact, there exists a critical adhesive strength, which is a function of the scratching depth at which the material removal mechanism transitions from atom-by-atom to wear fragment removal. This study confirms that the contribution of the atom-by-atom removal mechanism is suppressed when plowing dominates the process of material removal. Scratching polycrystalline Fe [296] showed that reducing the mean grain size changes the plastic deformation mechanism from a dislocation-mediated mechanism to one involving grain lattice rotation and grain boundary sliding. As the scratching depth increases to a certain level, radial, median, and lateral cracks can initiate, and there exists an optimal indenter half-apex angle, which depends on the scratch depth and indenter radius, where it is most difficult to initiate radial and median cracks [297].

Most recent atomic-scale scratching research has focused on FCC and BCC metals, whereas HCP struc-
tures, such as titanium (Ti) and magnesium (Mg), have received less attention, despite their technological
relevance, e.g., in lightweight components. Atomistic scratch simulations of Ti single- and bi-crystals [298]
revealed that softer Ti grains exhibit shallower residual scratch grooves than initially harder ones, and
this could be attributed to fragmentation-induced internal interfaces, see Fig. 11. An in-depth analysis of
scratching across a grain boundary revealed grain boundary bending (panel a), asymmetric pile-up formation
(b and c) and the formation of nano-grains immediately at the grain boundary (panel d). It was generally
shown for HCP metals that the basal surface has lower hardness than the prismatic surface, on which plastic
activity is more pronounced and from where dislocation loops are emitted into the substrate [299].

Studying the scratching process of a deformable substrate with a rigid asperity led to the formulation of
an analytical model that predicts wear profiles and substrate penetration depths as a function of load and
asperity geometry [300]. Interestingly, the experimentally validated model revealed a nonlinear relationship
between wear volume and force, differing from the classical Archard wear model. This emphasizes the
crucial role of asperity geometry in material removal rate and mechanism. Recent studies [301, 302, 303]
have shown that wear volume in single-crystal materials varies greatly with crystal orientation and scratch
direction, challenging the notion of a unique wear coefficient for a given material. This can be understood
as plastic hardening and size effects are strongly influenced by the crystal orientation and the direction of
applied forces [304]. Variations in these factors significantly affect the material's response to wear, leading
to differences in wear behavior. Furthermore, it has been demonstrated that while smaller grains result
in increased scratch hardness, they have negligible impact on the pile-up topography [305]. This suggests
that although finer grains enhance the material's resistance to scratching, they do not significantly alter
the deformation pattern that forms around the scratched area. These observations suggest that indentation
hardness alone may be insufficient to characterize the wear response of crystalline materials using the Archard
wear model, as it fails to account for crystal anisotropy. Alternatively, scratch hardness has shown promise
as a more effective measure for estimating wear responses that are both size-dependent and anisotropic,
<figure><img src="image_12.png" /><figcaption>Figure 11: Visualization of how scratching across a grain boundary in hcp Ti leads to persistent grain boundary (GB) bending (a), the asymmetric pile-up formation in front of the scratch tip (b) further evidenced by the atomic displacement vectors (c), as well as the subsequent formation of nano-grains (NG) near the grain boundary (d). Adapted from [298].</figcaption></figure>

proving useful at both the continuum level [306] and the atomic scale [307]. These insights provide a novel and efficient method for indirectly characterizing material properties, such as elastic modulus, hardness, hardening coefficient/exponent, and fracture toughness through scratch testing [308, 309].

### 4.1.2. Friction and wear of metal alloys

While the considerations discussed above for pure metals can provide deep insights into friction and wear mechanisms, most metals in real-world applications are alloys (and sometimes even amorphous) leading to failure mechanisms that can differ considerably from those in elemental metals. Chandross and Argibay [310] proposed a model for predicting the shear strength of polycrystalline metals and alloys based on the stress required for deformation via amorphization, which does not occur in pure metals, thus eliminating the need for adjustable parameters. They further extended this concept to predict the ideal strength of metallic glasses using stress-activated transformations with a regular solution model, offering a possible pathway to alloy design at low computational cost [311]. These methods show promise for designing and optimizing high-performance alloys and composites for maximum strength by optimizing the readily-calculated heat of fusion rather than performing resource-intensive simulations.

Going beyond the properties of single crystals by modeling alloys at the nanoscale under tribological loading conditions has yielded valuable insights into their interfacial properties. As a prominent example, multi-principal element alloys (MPEAs) or complex concentrated alloys (CCAs), often referred to as high- and medium-entropy alloys (HEAs and MEAs), have shown potential as robust materials for extreme environments, with increased strength and creep resistance compared to traditional alloys [312]. These materials also offer the possibility of being produced from scrap [313]. Atomistic modeling techniques advanced the understanding of the tribological behavior of these alloys, for which the reader is also referred to a recent review dedicated to HEAs [314]. Most of the particular studies outlined in the following used ball-on-flat contact geometries, that may be interpreted either as single-asperity contacts, nanoindentations or scratch tests [315], or even in the context of nanomachining.

Recent MD simulations revealed that a decrease of the chromium content to 22% in an FCC NiCoCr MEA leads to lower coefficients of friction against a diamond sphere, but a notably increased hardness when compared to the traditional equiatomic MEA composition [316]. In this case, the rearrangement of atoms at the sliding interface leads to a repulsive layer that is enriched in nickel, improving the tribological properties of the alloy. Other atomistic simulations with polycrystalline CuCoCrFeNi HEAs demonstrated
<figure><img src="image_13.png" /><figcaption>Figure 12: MD-simulated grain orientation (EBSD) tomographs of polycrystalline CuNi25 after sliding against a rough rigid counterbody (gray) for 280 nm at a normal pressure of 0.7 GPa at different velocities $v_x$. The three panels show evidence of (a) dislocation mediated plasticity, (b) fluid-like flow (note the high contact conformity), and (c) relocalization of the sliding interface to liquid contact regions at ballistic speeds. Images adapted from [325].</figcaption></figure>

that stresses induced at grain boundaries under cryogenic conditions (77 K) can reduce stresses within the grains, enhancing wear resistance [317]. While single-crystal samples exhibited the lowest friction, columnar polycrystalline samples with a grain diameter of 10.7 nm showed the highest wear resistance as well as the best overall performance for friction and wear. Nanoscale modeling has also provided insights into the tribological performance of HEAs used as coatings. An additively manufactured AlCoCrFe coating with a thickness of 3.2 nm on pure aluminum increased the Young's modulus of the sample by 50%, enhanced its hardness by a factor of ten, reduced the wear scar volume by an order of magnitude, and lowered friction by 75% compared to the uncoated Al substrate [318]. Similarly, a FeNiCrCoCu coating on copper exhibited protective properties through stress absorption and subsequent release [319]. The stress release leads to the generation of stacking faults in the coating, that then act as a barrier to lattice damage in the Cu substrate. In nanotwinned Cu/FeCoCrNi nanolaminates, the coefficient of friction, wear rate, and surface roughness exhibited relative minima at a sliding speed of 100 m/s, attributed to the formation of V-shaped dislocations and stacking fault tetrahedra [320].

Additional alloys, such as Invar (FeNi36), AZ91 (simplified to Mg₁₇Al₁₂), and cupronickel (CuNi) have been recently studied for their tribological behavior using atomistic modeling. Studies of the polishing of the Invar alloy have focused on optimizing conditions to reduce subsurface damage and improve surface quality [321, 322]. The authors found that a rolling speed of the abrasive particles that corresponds to pure rolling is most detrimental to both surface roughness and subsurface damage, while some slip (in either direction) produced better-quality polished surfaces. Based on the MD results, an analytical model was developed to predict the minimum uncut chip thickness (an indicator of a material's machinability) as a function of the elastic modulus, Poisson ratio, shear strength, and hardness of the processed material as well as geometrical considerations. Another study investigated the influence of continuous and discontinuous precipitates on scratch-induced wear in the AZ91 magnesium alloy through experiments and atomistic simulations [323]. The orientation of the BCC precipitates on the surface significantly affected the alloy's deformation behavior; a longitudinal orientation with respect to the cutting direction led to a ductile-brittle transition, while a transverse orientation yielded a fully ductile response. Discontinuous precipitates exhibited higher wear resistance, but also higher friction regardless of orientation compared to continuous precipitates, a conclusion that was supported by an analysis of planar defects and dislocations. Finally, a series of large-scale MD simulations on polycrystalline CuNi alloys have explored deformation regimes and created deformation mechanism maps [18] that can be considered an enhancement of Ashby charts (a materials selection tool) for tribological systems. The validity of the obtained maps was assessed via comparisons to SEM images of worn samples from ball-on-plate experiments, allowing researchers to identify different deformation mechanisms (twinning, grain refinement, shear-induced grain growth, etc.) as a function of alloy composition, normal load, temperature [324], and sliding velocity. At extreme sliding velocities (~320 m/s), a regime characterized by fluid-like behavior of the softer sliding partner with a highly conformal interface was observed [325], see the differences between panels (a) and (b) in Fig. 12. Further investigation at ballistic speeds (> 640 m/s) revealed a rise in contact temperature, the relocalization of the sliding interface
to liquefied metal contact patches (see panel c), as well as a decrease in plastic deformation and interfacial friction, that could be aligned with experimental observations by Bowden and Freitag [326].

Friction and wear modeling of single-crystal metals at small scales has reached a mature state, with well-established mechanisms and tribological laws governing contact, friction, and wear. However, numerical studies must now address more complex systems—such as polycrystalline metals and high-entropy alloys—where the detailed contributions of subsurface defects are critical. For high-entropy alloys, future atomistic modeling should also investigate dislocation mechanisms and interactions with lubricants or 2D materials like graphene [327]. Recent research [328, 329] further highlights the potential of tribology-assisted material characterization, particularly through scratching techniques, where combined topographical and force measurements facilitate the extraction of diverse material properties and microstructural features.

To upscale this understanding, the development of concurrent and sequential multiscale models is essential to bridge from nano- to macro-scale information. These models will enable a more comprehensive description of friction and wear mechanisms across different length scales. To ensure their accuracy and reliability, systematic small-scale experiments examining factors such as contact area, pressure, surface topography, and lubrication conditions are necessary for rigorous validation. Addressing these challenges will pave the way for more robust numerical models, ultimately advancing materials design, durability, and component lifespan.

## 4.2. Tribochemistry

Tribochemistry is the study of chemical reactions occurring at the interface between surfaces in relative motion. Tribochemical reactions are caused by interfacial stresses and/or frictional heating and in return lead to changes in friction, wear, and lubrication. While mechanical alloying or even the motion of dislocations at interfaces described in the previous sections is a type of tribologically-induced reaction, because existing chemical bonds are broken while new ones are formed, mechanochemistry encompasses a broader range of stress-induced chemical reactions, including those occurring in bulk solids, liquids, or gases due to mechanical stresses caused by grinding, milling, or shearing but also those triggered by forces on individual molecules. Whether a change of bond topology should be classified as tribochemical, plastic, or more generally mechanochemical may not always be unambiguous. This is also because the degree with which chemical bonds are directional even in elemental solids changes almost continuously throughout the periodic table.

One of the tribochemical reactions, which frequently affects the lives of car engines, is the degradation of engine oils. At the same time, effective lubrication by engine oils relies on tribochemical reactions such as those leading to the formation of surface tribofilms from antiwear additives and friction modifiers. The question how to tribochemically “synthesize” in situ tribofilms that protect surfaces from wear and reduce friction in boundary lubrication has become increasingly important. The goal is to replace expensive coatings of limited life span with tribofilms that are self-generating and self-repairing [330, 331, 332]. The classical example is the anti-wear film formation of zinc-phosphates, which is still controversially discussed and hence reviewed below along with the growth of low-friction carbon tribofilms, which has also enticed much recent interest.

While the interpretation of tribochemical pathways is much assisted by atomistic simulation, attempts are also made to develop a formal theory for tribochemical reactions. The field is still in an early state, perhaps also because, with rare exceptions [333], most theories [334] overlook that stress is not a number but a (symmetric) tensor of rank two and that different stress-tensor elements can trigger different chemical reactions (see Fig. 13a-d). Consequently, formal theories such as the von Mises yield criterion in solid mechanics, which are independent of coordinate system choices, have yet to be formulated.

### 4.2.1. Generic models for tribo- and mechanochemistry

According to Eyring’s transition-state theory [337], the central quantity affecting the rate of gas-phase reactions $v$ from reactants to products is the activation energy $\Delta E$, as $v$ depends on it exponentially

$$v \propto \exp(-\beta\Delta E), \qquad (5)$$
<figure><img src="image_14.png" /><figcaption>Figure 13: Upper row: Educt (a) and products (b-d) of strain-induced reactions in a model ZDDP decomposition molecule. Imposed strains differed in shape and thereby lead to different hydrostatic pressures (p) and shear-stresses (τ), where the reactions took place, but also to different bulk moduli B, as indicated in the panels. Adopted from [335, 336]. Lower row: Tight-binding MD simulation of the tribochemical composition of a glycerol molecule between two ta-C surfaces [12]. The glycerol molecule chemisorbs simultaneously to both surfaces (e). Tribologically induced molecular decomposition combined with interface cold welding leads to the formation of a H- and O-containing amorphous carbon sliding interface (f). Localization of the plastic shear deformation results in the formation of a passivated, partially aromatic sliding interface (g), which lies at the root of superlubricity. The polycyclic structure of the surfaces can best be seen in panel h, which is a top view of a the atoms in 2.8-Å-thick region of the bottom surface.</figcaption></figure>
with $\beta = (k_B T)^{-1}$. Due to a directed force $F$ pulling a bond in a molecule, the reaction barrier to bond dissociation $\Delta E$ can be reduced to the extent that it approaches the thermal energy $k_B T$ or that $\Delta E$ even becomes negative, in which case the reaction, e.g., the breaking of bonds, happens (quasi-) spontaneously. In leading order, it is described through a linear relation in the so-called Bell model [338], though previously proposed by Evans and Polanyi [339], Kauzmann and Eyring, [340] and others

$$ \Delta E = \Delta E_0 - F\Delta x, \qquad (6) $$

where $\Delta E_0$ is the energy barrier in the absence of an external force, while $\Delta x$ could be coined a reaction length.

In a bulk material, stresses rather than forces affect the rate of reaction, so that $\Delta E$ is a function of pressure or stress. In mechanochemistry, one possible and frequently used generalization of Eq. (6) reads

$$ \Delta E = \Delta E_0 + \alpha \Delta H, \qquad (7) $$

where $\Delta H$ is the reaction enthalpy. However, other generalizations are possible. To this end, assume the force $F$ in Eq. (6) is applied to an initially cubic box with volume $V = a^3$ containing an isotropic material on the surface normal to the $z$ direction yielding a tensile stress of $\sigma_{zz} = F/a^2$. The most general expansion of $\Delta E$ in the stress-tensor invariants $I_n$ yielding linearity in $F$ is [341]

$$ \Delta E = \Delta E_0 + \sum_{n=1}^{3} \alpha_n I_n^{1/n}. \qquad (8) $$

Here, $I_1$ reduces to $\sigma_{zz}$, in which case the hydrostatic pressure becomes $p \equiv -I_1/3 = -F/(3a^2)$, while the von Mises or equivalent shear stress $\sigma_{vM} \equiv \sqrt{3}I_1^2 - I_2$ would be $\sigma_{vM} = F/(3a^2)$.

For a general loading condition, there is no reason for $\sigma_{vM}$ to be equal in magnitude as $p$. A proper generalization of the linear Bell relation would read

$$ \Delta E = \Delta E_0 + pV_1 - \sigma_{vM}V_2 - \sqrt[3]{I_3}\alpha_3, \qquad (9) $$

where the $V_n$ are expansion coefficients, which one may want to call free volume in the case of $V_1$ and associate with a volume needed for a reaction to take place. For non-isotropic systems, up to six volumes can be defined [342], yet, most studies focus on $V_1$ and associate it with an activation volume. A difficulty for the description of a tribochemical reaction is that stress tensor elements and the $V_n$ coefficients change discontinuously across the interface so that they are generally undefined right at a tribological interface. In this case, one may want to replace $p$ with the normal load per unit area and $\sigma_{vM}$ with the frictional shear stress [333], and ignore the effect of additional or higher-order strain tensor elements on the reaction rates.

Classic evidence for stress-assisted bond breaking in sliding contacts comes from the observation of wear during sliding atomic-force microscopy tips, such as by Park et al. [343] for the formation of wear tracks in ionic solutions or by Gotsmann and Lantz [344] for wear of the tip itself under dry conditions. Gotsmann and Lantz characterized wear indirectly from tip radii extracted in pull-off experiments. They were the first to show that the worn volume upon sliding tips for hundreds of meters on a polymeric substrate follows the Bell model rather than the phenomenological law of Archard [78] stating that wear is proportional to the work done by friction. This was later confirmed in direct electron-microscopy measurements of the tip geometry by Jacobs and Carpick [342], albeit for wear rates orders of magnitude larger than the previously observed atom-by-atom attrition [344]. We now understand that not only wear, but also tribofilm growth [345] is governed by such stress-assisted activated processes. As many other works, all of the works cited in this paragraph [344, 342, 345] consider only the pressure and ignore the other terms in Eq. (9), although Jacobs and Carpick explicitly recognize that as an approximation.

The need to incorporate a more rigorous description of the local stress field in the models is exemplified by the lack of consensus on which particular scalar pressure value (e.g. normal or shear interface pressure) determines the atom-by-atom wear process. For example, Gotsmann et al. [344] assumed that the energy barrier for the wear process is reduced by the interface shear stress that stretches the bonds at the interface
between tip and substrate. In contrast, Jacobs et al. [342] use the compressive normal contact stress to fit the experimental wear data, but note that the shear stress used in previous works scales linearly with the compressive stress. Although these studies are unable to identify the precise reaction pathway for atom-by-atom wear, Jacobs et al. [342] propose that each single wear process could be decomposed into two steps. In the first step, compressive stress at the contact between tip and substrate promotes the formation of covalent bonds between an atom of the tip and one of the substrate. In this case, this is the rate-limiting step. In the next step, the bonds connecting the first of the two atoms to the tip are broken upon sliding of the tip on the substrate.

The idea of decomposing the atomic wear process into a bond-formation and bond-breaking (or wear) step proposed by Jacobs et al. [342] has been recently resumed and elaborated in works aimed at modelling the load [346] and speed dependencies of atomic wear [347]. While the work of Raghuraman et al. [347] focuses on the removal of oxygen atoms from graphene by means of a silicon tip, and the wear step is a surface hopping process that is biased by the tip lateral velocity, the study of Wang al. [346] is closer to the idea of Jacobs and Carpick. The authors assume that the energy barrier for the bond-formation step is lowered by a compressive stress, while the barrier for the bond-breaking step is reduced by the local shear stress. The energy barriers for the two reactions are the only fitting parameters as the authors directly estimate the work done by the compressive and shear stress components using simple models for contact deformation energy and bond breaking, respectively, without introducing activation volumes. Evaluating the wear volume boils down to finding a suitable expression of the normal-load dependency of the real contact area. Comparison with the results of MD simulations on DLC/DLC contacts with different load dependences of real contact area are accurate and show that this dependency determines the wear law.

In their two-step atomic wear model, Wang et al. [346] use simple estimates for the mechanical work that lowers energy barriers and avoid introducing an activation volume, which they claim is a fitting parameter with no explicit physical meaning. In most of the other studies mentioned in this section, key discussions revolve around the interpretation of the volume term $V_1$ in Eq. (9). It is sometimes referred to as the (negative) activation volume. Interestingly, its value for a given interface can take values ranging from 6.7–115 Å³ [342, 344, 348, 349, 347], making it difficult to assign a local interpretation to $V_1$. One proposed explanation is that $V_1$ does not only have an interfacial but also a bulk contribution, where the latter can depend on the radius of curvature of the scratching tips or the stiffness of the material close to the interface [348]. However, this suggestion violates the principle that chemistry is generally local. It may well be that the locality principle can be reestablished by including all relevant stress terms into a theory, e.g., the $I_3^{1/3}$ term or by properly deducing $\sigma_{vM}$, which requires the generally neglected in-plane stresses to be deduced, for example, through a proper contact-mechanics analysis.

### 4.2.2. Anti-wear tribofilms

Many studies on the formation of anti-wear films focus on the rate of tribofilm formation but ignore that the products themselves depend on the stress that acted when the films were produced, which violates the assumptions made for Eq. (7). For example, decomposition products of anti-wear additives, most notably of zinc dialkyldithiophosphates (ZDDPs), form substantially stiffer polyphosphate glasses with longer chains on top of asperities, where stresses had been previously high than in the valleys between asperities with more moderate stress history [350]. Nonetheless, it is certainly interesting to learn that zinc-phosphate tribofilms can grow without direct asperity contact provided that the shear stress is sufficiently high: under identical normal stress and shear rates, Zhang and Spikes [351] found film formation with a high-viscosity but not with a low-viscosity base oil. Fang et al. [352] corroborated these results. By blending together high-viscosity and low-viscosity base oils, they systematically controlled both the shear stress and normal stress. In agreement with Zhang and Spikes [351], they found that increasing the temperature and shear stress increase tribofilm growth exponentially; however, increasing the normal stress had the opposite effect.

In another recent study, Zhang et al. [353] compared the tribofilm growth rates of ZDDPs with a range of alkyl substituents. They showed that the ZDDPs with long linear chains had a stronger temperature dependence and thus a larger activation energy, while those with branched and bulky groups had a stronger stress dependence, and thus larger activation volumes. Zhang et al. [354] also showed that ZDDPs with secondary alkyl groups formed tribofilms much faster than those with primary ones under a wide range
of temperature and shear stress conditions. Both ZDDPs had similar a activation energy and activation
volume and the main reason for the larger reactivity of the secondary ZDDPs was a larger pre-exponential
factor [354].

As well as for its growth, mechanochemistry also plays a role in the mechanical strengthening of ZDDP tribofilms. Mosey *et al.* [11] studied the pressure-induced changes in zinc phosphates using DFT-based MD simulations. They showed that, at very high hydrostatic pressure, the coordination of zinc atoms increases from tetrahedral four coordination to six coordination and back to four, albeit in a seesaw geometry. The irreversible changes lead to cross-linking, network formation, or alternatively to longer phosphate chains [336], and an increase in elastic modulus. Such transformations could be important to the durability of ZDDP tribofilms. It was argued that the pressure required for this transformation (17 GPa) was beyond that typically found in tribometer experiments or engineering components [351]. However, a follow-up study by Shakhvorostov *et al.* [355] showed that the pressure required for pressure-induced cross-linking was actually much lower than the initial estimate (7 GPa). Moreover, a recent study by Sukhomlinov and Müser [335] showed that stress anisotropy (e.g shear stress) further reduces the required network formation pressure to levels expected inside rubbing steel contacts (1-2 GPa), thereby strengthening the hypothesis that two and potentially even more terms on the right-hand-side of Eq. (8) matter. They also showed that the tribofilms formed under compression are stiffer than those formed under shear, which suggests that although normal stress may inhibit ZDDP tribofilm growth [352], it may be important to obtain durable tribofilms.

A combined experimental and multiscale modelling study by Ruiz *et al.* [356] showed how contact stiffness can be another important consideration in ZDDP mechanochemistry studies. They used macroscale tribometer experiments, contact mechanics calculations, and tight-binding (TB)-DFT-based NEMD simulations to study the performance of ZDDP on diamond-like-carbon (DLC) contacts of varying hardness. They showed that an optimal combination of ultralow friction and negligible wear can be achieved with ZDDP-containing lubricants by using hydrogen-free tetrahedral amorphous carbon (ta-C) with moderate hardness. Softer coatings exhibited low wear but higher friction, while harder ta-C underwent severe wear and subsurface sulphur contamination. Similarly, Peeters *et al.* [357] used a combination of AFM experiments and DFT-based NEMD simulations to study the adsorption and decomposition of ZDDP on lightweight metallic substrates. The AFM experiments showed that ZDDP gives poor wear protection on aluminium and magnesium, but performed better on aluminium-magnesium alloys [357]. The DFT calculations showed that this was due to the higher surface energy and affinity of the sulfur atoms in ZDDP to magnesium atoms in the alloy surface [357]. This provides anchor points on the opposing sliding surfaces, enabling efficient transmission of stress to the molecules through the sliding surfaces, and thus higher mechanochemical reactivity. This is consistent with previous TB-DFT-based NEMD simulations of another type of additive, organic friction modifiers on sliding ta-C surfaces [12].

Another notable study in this area is by Boscoboinik et al. [358], who used a combination of AFM experiments in ultrahigh vacuum and quasi-static DFT calculations to investigate the mechanochemical decomposition of methyl thiolates on copper surfaces. They showed that the normal-stress dependency of the activation energy was in excellent agreement between the two techniques [358]. In a follow-up paper, Rana et al. [359] studied the mechanochemical decomposition of other alkyl thiolates on copper surfaces in ultrahigh vacuum. A buckling theory analysis of the effect of a normal stress on an adsorbate that is oriented perpendicularly to the surface that reacts by tilting suggests that a critical value of the stress should be required to initiate a mechanochemical reaction. This concept was verified by using DFT calculations to simulate the effect of compressing a homologous series of alkyl thiolate species on copper [359]. This predicts that a critical stress is indeed needed to initiate methyl thiolate decomposition, which has a perpendicular C-CH₃ bond. In contrast, no critical stress is found for ethyl thiolate with an almost horizontal C-CH₃ bond, while a critical stress is required to isomerize propyl thiolate from a trans to a cis configuration. These predictions were tested by measuring the mechanochemical reaction rates of these alkyl thiolates on a copper substrate by sliding an AFM tip over the surface [359]. In accord with the DFT predictions, a critical stress of 0.43 GPa for methyl thiolate and 0.33 GPa for propyl thiolate was found, but no evidence of a critical stress for ethyl thiolate.
### 4.2.3. Low-friction tribofilms

Lubricants generally contain other surface-active additives in addition to anti-wear additives. An important example are friction modifiers (FM) that reduce friction in boundary lubrication through the formation of low-friction tribofilms. In recent years, numerous studies combining experiments and atomistic simulations have focused on low-friction tribofilms made of carbon. The formation of carbon tribofilms is often driven by tribochemical reactions involving organic FM [12, 360, 361], molecules in the gas phase [331, 362] or film transfer from carbon surfaces [363]. The stable chemical passivation and the smoothness of amorphous carbon tribofilms make them ideal candidates to achieve superlubricity [12, 360].

Molecular modelling of carbon tribofilms from FM focuses on three aspects: the surface adsorption of FM [364, 365], the tribochemical processes leading to the film formation [12, 360, 361] as well as the relationship between friction and the chemical structure of the tribofilm [366, 365]. Broadly speaking, two different tribofilm structures are identified in molecular modelling studies. On the one hand, there are molecular films composed of amphiphilic molecules (e.g. oleic acid) that chemisorb onto steel surfaces via their polar group [366, 365]. On the other hand, there are polymeric [332] or amorphous carbon tribofilms that form through the tribochemical degradation of the FM [12, 360, 361].

Quantum-mechanical or reactive MD simulations show two different mechanisms for the surface passivation of amorphous films. The first mechanism is the saturation of surface dangling bonds with chemical species (e.g. H, OH, O) originally present in FMs [12]. The other passivation mechanism is the formation of aromatic surface structures through the coalescence of aromatic molecules [361] or through plastic deformation of amorphous carbon films containing impurities like O, N and H from the tribochemical fragmentation of FM [12, 360, 367] (Fig. 13e-h).

While atomistic simulations have helped shed light on many tribo-chemical mechanisms involved in the formation of tribofilms, their application to complex tribo-chemical systems is not without its limitations. The tribofilm structure depends on several factors, like load, surface chemistry, synergies and competitions between additives. A very recent experimental study by Song et al. [368] shows that even a slight difference in the position of a carbonyl group in two isomers (glycerol monostearate and stearyl glycerate) can result in different tribofilm structures (amorphous or brush-like) and hence in different friction coefficients.

Yet, in molecular modelling, the tribofilm structure seems to be strongly influenced by initial assumptions and simulation setups. For example, a frequent hypothesis is that FMs diffuse through the lubricant and reach the surface where they physi- or chemisorb with no need for mechanical activation. Adsorption energies and configurations can be obtained with quantum-mechanical simulations (adsorption enthalpy) [365] or free energy calculations based on classical force fields (physisorption free energy) [364]. Based on the resulting adsorbed configurations, molecular tribofilms with a given surface coverage (usually a free parameter) are “constructed” and then used in non-reactive MD simulations to study structure-friction relationships [366, 365]. These simulations, however, cannot provide information about the chemical stability of the molecular tribofilms (or their oligomerization) under tribological load.

Conversely, amorphous tribofilms nucleation is often modelled as a sequence of tribochemical FM fragmentation reactions followed by the shear-induced plastic deformation of the nascent amorphous film. This is usually the result of reactive MD simulations in which a small, arbitrary number of molecules (often one in DFT simulations) are squeezed between two reactive surfaces in relative motion and cannot escape the periodic, flat contact [12, 360] (e.g., Fig. 13a-d). In systems for which the accuracy of reactive interatomic potentials is sufficient, this approach can be applied to dozens of molecules for several nanoseconds. For example, this approach provided insight into the tribochemical steps involved in the nucleation of carbon-based tribofilms by shear-induced oligomerization of hydrocarbons [330, 369]. In particular, these works emphasize the role of shear stress in enabling tribochemical reactions at relatively low temperatures. The chemical iteration between molecule and surface plays a key role in the activation, shear-driven deformation, and dehydrogenation of hydrocarbons, which undergo fragmentation and react with other molecules to form oligomers that can act as tribofilm nuclei.

Ideally, to minimise the influence of the initial assumptions and be able to cope with complex but essential aspects like competition/synergy between additives or molecular size effects, it would be necessary to merge the modelling strategies outlined above. This would imply having models with the efficiency and
wide range of chemical elements of force fields and the chemical accuracy of quantum-mechanical methods, at least in localized regions of the system. Recent attempts have been made to apply hybrid QM/MM (quantum mechanics/molecular mechanics) methods to tribochemistry simulations [370]. However, the large size of the QM regions, also caused by the non-trivial handshaking between the QM and MM regions [371], does not allow for significant efficiency gains compared to traditional QM simulations. A possible way to achieve chemical accuracy more efficiently could be provided by the rapid development of machine-learning interatomic potentials that the user can train for specific applications [20, 21, 22]. The trade-off between the accuracy and efficiency of these methods in tribochemical applications that are chemically complex and far from the equilibrium (i.e. the standard training conditions) are now under scrutiny [372].

In spite of recent developments in the field of interatomic potentials, modelling the growth of a tribofilm usually requires time scales that are beyond those currently allowed by MD. Alternatives based on the combination of MD and kinetic Monte Carlo have recently been proposed to model the thermal growth of polyphosphate films on iron surfaces [373]. However, the development of analogous methods for non-equilibrium systems with complex tribochemical reaction paths seems unlikely. An alternative “reverse engineering” approach is to focus on the final tribofilm structure instead of the formation mechanisms. Recently, Zarrouk et al. [374] have proposed a scheme to infer the structure of oxygen-rich amorphous carbon films starting from experimental X-ray photoelectron spectra (XPS). The scheme combines grand-canonical Monte Carlo, a machine-learning model for prediction of XPS spectra, and a machine-learning interatomic potential to generate O-containing amorphous carbon structures that match the experimental XPS spectra and are “low in energy”.

The overview proposed in this section is obviously far from complete but highlights some general open questions. Generic models for tribochemistry and the interpretation of tribochemical reactions may benefit from more rigour. For example, models for tribochemical reactions that emphasize the role of the stress tensor ignore the presence of surfaces and chemical heterogeneity, while those including surfaces generally assume in-plane strain, which may not always be a justified approximation. Aside from questions about boundary conditions and accuracy of atomistic techniques under mechanically and chemically complex tribo-chemical conditions, there are also open questions on the assumptions made in preparing representative atomistic systems. These include: the thickness of the lubricant film in which the tribo-chemical reactions take place; its composition and the state of the additive molecule (e.g. solvated in oil, physisorbed or chemisorbed on the surface); interactions between different surface-active additives; the chemistry and passivation of the surface itself (e.g. metal, one of its oxides, a previously formed tribofilm). Overall, it is evident that approaching complex tribo-chemical problems requires a multidisciplinary and multiscale approach. This is being recognised both from a modelling and experimental perspective and put into practice in an increasing number of studies that are making progress by combining tribometry, chemical and topographical surface characterization with fluid dynamics, contact mechanics and molecular modelling (e.g. Ref. [330, 362, 360, 361]).

# 5. Applications

The last section of this review illustrates a few modern application domains of tribology. We have selected examples that present interesting opportunities for the modelling community.

## 5.1. Tunable interfaces

Friction is often considered to be predefined by surface chemistry and structure. Yet, tuning friction, especially dynamically, has many applications ranging from clutches to haptic displays. Over the last decades, some interesting approaches towards tunability have emerged. While many of those have been experimentally realized, an understanding of the fundamental dissipation processes, as well as support from modeling for further optimization, is often lacking. Tuning friction traditionally takes two routes: either through morphology or chemistry. Recently, dynamic control with external electric or magnetic fields has also been pursued [375]. As we discussed in the lubrication section above, surface texturing is a traditional way to use morphology to tune the frictional behavior of lubricated interfaces.
<figure><img src="image_15.png" /><figcaption>Figure 14: Interface with strain-induced surface patterning. Lateral strain modifies the geometry and wetting behavior, offering a potential route to tunable tribological properties. Such interfaces have received little attention in modeling. Images adapted from [380].</figcaption></figure>

The tunability of interfaces, for example, in the concept of switchable adhesives achievable via the “manipulation of near interfacial chemistry and dissipation, control of contact area, or through mechanical or structural switching of adhesive compliance or geometry” [376], is related to adhesive failure and governed by fracture mechanics. Even simple geometric structures, bumps on a surface can be use to tune the friction force vs normal load response of soft interfaces [377]. Other means to control friction make use of biomimetics, for example, to design kirigami-activated textured surfaces on highly deformable materials that mimic the scales on the skin of snakes in achieving directional frictional properties [378]. In this case, geometry rather than the constitutive properties of the material dictates the frictional performance. Since geometry can be adapted on-the-fly, interfacial properties of such textured surface are in principle switchable. This is the case of morphing interlockers, enabling switchable mechanical adhesion [379]. Figure 14 shows an example of a surface with switchable wetting properties [380].

Similar principles apply to the class of architected materials or metamaterials that add functionality to materials through their geometry rather than their material properties. An interesting application in this direction is the use of friction within the internal features of a metamaterial unit cell to dissipate energy [381]. Reversible energy dissipation can be achieved in these metamaterials even with metals and ceramics. Analytical and finite element models have been shown to trace the hysteresis observed in experimental results, but the use of constant friction coefficients to capture “internal” contacts in these models appears to be insufficient in accurately matching the results. Another way of frictional control stems from the structural altering of contact interface pairs using combinations of low and high friction materials. Such an alteration could be achieved by mechanical actuation such as, for example, inflation in soft robots, which can ensure heterogeneous friction across the macroscopic contact interface and promote a directional motion, as was done for worm-like locomotion [382]. Miniaturization of such a frictional control can lead to a smart friction control in a more continuum sense. In contrast to asymmetric friction mentioned above, a heterogeneous frictional control can ensure locomotion in any direction.

The aforementioned contact-related architected materials are relevant even without internal frictional dissipation. Even the inclusion of frictionless internal contact interactions, resulting in strong anisotropic asymmetry of elasticity, provides interesting emerging phenomena including wave-control, filtering and absorption [383]. Other interesting implementations of the contact interaction inside the material include hierarchical woven structures [384], interpenetrating lattices [385], entangled “granular” structures [386], and topologically interlocking materials [387, 388].

Yet another way of friction control comes from electromagnetic field operating at mechanisms from atomic to structural scales and spanning triboelectric generators, geological drilling, automotive braking and efficiency, spacecraft systems, biological systems, and magnetic spintronics and others; 25 interesting examples are provided in a review paper [375]. New developments include the realization that voltage-induced fluctuations in the electronic properties can enhance energy barriers [389] and a more detailed understanding of how electrostatic fields can induce dipole reorientation thereby allowing dynamical switching between smooth sliding and stick-slip motion [390]. However, one may want to address the question in this but also related cases, if electrostatic fields of 4 GV/m or greater are meaningful given that the breakdown field
strength of air (3 MV/m) or even that of diamond (2 GV/m) are clearly exceeded. Finally, the vibrational control of friction, known since the 50s [391, 392], to the best of our knowledge did not bring breakthrough discoveries but some micro- and structural-scale models were improved [393, 394].

The last decade has also seen the development of several solutions based on hydrogels, which have often been adopted to obtain dynamic friction control by manipulating the interplay between the polymer network, its composition and surface properties, and the surrounding fluid [395, 396, 397]. It has been shown that tunability can be obtained e.g. through different gel's networks and functionalisations and by manipulating external fields [398, 399] and that surface properties, such as roughness and pillar geometry [400], significantly affect the contact area and the mechanical response, leading to complex friction behaviors. Models have been successfully developed to explain some of the experimentally-observed mechanisms, e.g. pH-tunable friction [401], or the reduction of friction and wear in composite hydrogels for cartilage mimics [402]. However, much of the modeling work does not address the intimate coupling between stimuli and function modulation, therefore failing to provide universal understanding and better design strategies for hydrogel-based solutions; progress in this area could certainly unlock their full potential for several soft contact applications.

Switchable properties can also be achieved in microgels, which are small particles of hydrogels. Understanding and control of their tribological properties saw much progress in recent years. In Ref. [403], the authors developed dual-responsive microgels, surface-functionalized with PEG or PDMAEMA, which allow for independent tuning of swelling behavior and surface properties in response to thermal and pH changes, which demonstrate a unique decoupling of the two attributes. Modelling and simulations of such microgels could provide extra insights into their behavior and enable tailoring these materials for specific use-cases. Similarly to hydrogels and microgels, the use of polymer brushes has been recently advocated to manipulate interfacial friction, especially in boundary lubrication in the presence of soft contacts [404]; however, the use of coarse-grained MD simulations to predict their behaviour has resulted in discrepancies between computational and experimental observations, which can only be resolved by producing models that aim to simulate experimentally imposed velocities and obtainable grafting densities [405].

Liquid crystals (LCs) have also been adopted to tune friction, most often by inducing viscosity modulation in fluidic layers via the addition of nematic and rod-like LCs to the lubricating fluids [406, 407]. The mechanisms for tunable friction include shear-induced alignment, i.e. when the shear flow can cause the LC molecules to align with the direction of the flow, and control through external fields (see e.g. [408]), whereby the LC can be oriented by applying an external electric or magnetic field, which affects viscosity and, consequently, friction. Despite a number of attempts to use theory and simulations to shed light on LCs behaviour and their effect of various lubrication regimes [409, 410], and the fact that molecular simulations have provided fundamental insights into how tunable friction works with LCs, limited research has been done using MD simulations to investigate the surface induced alignment of LCs and the individual and synergistic effect of topography and chemical functionality of a surface [411]. Progress in such simulations would provide unique insight and allow exploring properties that are difficult to measure experimentally.

The integration of LCs as active constituents in liquid crystal polymers (LCPs) has enabled remarkable control over surface characteristics. Specifically, LCPs can exhibit dramatic and reversible surface topographical changes when exposed to light [412] or temperature fluctuations [413], consequently allowing for tunable friction modulation. Expanding on this principle, numerous other stimulus-responsive polymers (SRPs), or smart polymers, are detailed in the literature. These advanced polymeric materials are characterized by their ability to undergo adaptive physicochemical transformations in response to various environmental perturbations. These perturbations include temperature, solvent ionic strength, modified solvent conditions (e.g., via co-solvent/co-non-solvent addition), UV-vis light irradiation, redox potential, and the application of electric fields. When anchored to surfaces, these SRPs impart dynamic control over interfacial friction and adhesion [414]; being able to model these complex interactions, which is yet unexplored in the tribology field, would provide us with better tools and more opportunities to develop new disruptive solutions to actively control friction.

Much progress was achieved in recent years on haptic feedback technologies for different kinds of human-machine interfaces. The tribological community contributes considerably to this progress with experiments and modeling. Modern technology does not reduce the feedback to mechanical haptic perception both force
generated [415, 416] and nerve stimulated [417] but can also include thermal effects [418, 419]; combined feedback is termed multi-mode haptic feedback. More relevant reviews can be found on this topic [420, 421], while numerous tribological studies on finger-surface interaction can also be found in the literature [422, 423, 424, 425]. Controlling friction due to the electroadhesion effect was studied experimentally and theoretically in [426, 427] and revealed the key role played by the skin surface roughness and the nature of the touchscreen coating. It was argued that reducing the effective thickness of the coating can drastically enhance the human tactile sensing experience.

## 5.2. Energy harvesting

In the past decade, there has been a significant surge in investigations of tribological interfaces in the presence of applied electric fields as well as in the exploitation of tribological interactions to harvest energy (e.g. with the proliferation of triboelectric nanogenerators - TENGs) and for the control of critical interfaces in both energy devices and components in which electric fields can be manipulated to achieve active tuning of friction. This section reviews recent advances made in these areas with focus on key fundamental discoveries and progress achieved using advanced theories and simulations. A critical overview of how advanced tribological models can help us obtain new disruptive solutions for energy devices is also provided.

An area that has seen an exponential growth in research activities in the last decade is triboelectricity (TE), which is the generation of static electricity through contact or rubbing. The phenomenon has supposedly been known much before Thales of Miletus (6th century B.C.), who is usually but falsely credited with its discovery [428]. While electron transfer clearly is the front-runner charge-transfer mechanism between metals, TE of insulators is frequently said to be poorly understood, despite a plethora of suggested theories [3]. Interestingly, Ernst's [429] generic but rigorous argument was forgotten for a long time: gradients in temperature, stress, and chemical composition and even more so their discontinuities at (tribological) interfaces all lead almost unavoidably to the flow of ions if present. An important point that may have been overlooked in this regard is that the nature and the number of ions produced from originally neutral (zinc-phosphate) molecules can sensitively depend on all three invariants of the stress tensor that acted during a (tribo-) stress induced chemical reaction [336]. Almost all TE theories only consider the properties of unstressed substances as reference. Thus, the question to be researched might not be what single mechanism explains all TE between dielectrics but which mechanism dominates under what circumstances. This offers great opportunities for the modeling of charge transfer using molecular simulations. However, popular reactive force fields fail to describe redox reactions and more generally electron transfer. Even DFT suffers from this problem as it equalizes the chemical potential of electrons throughout the simulation. Consequently, contact-induced electron-transfer between insulators could only be simulated so far using a generic, proof-of-concept potential to which the oxidation states of individual atoms was added phenomenologically as a discrete, internal atomic degree of freedom [430].

TE plays a crucial role in processes like printing, photocopying, electrospraying, and can lead to issues like electrical damage to microelectronics, disruptions in pharmaceutical processes, and increased friction and energy losses in industrial applications. It has only very recently (since 2012) been adopted in nano- to micro-scale devices to convert mechanical energy into electricity with the advent of TENGs [431]. Triboelectric nanogenerators are a promising technology that harnesses tribocharges for energy, enabling self-powered sensors, portable devices, and more, contributing to sustainable energy sources and the Internet of Things [432]. The fourth industrial revolution requires multifunctional arrays of sensing systems from air quality to in vitro health monitoring [433]. TENGs are ideally suited to provide a reliable and sustainable electricity supply for these systems by harvesting the mechanical energy that exists all around us. Networks of TENGs have already been successfully developed to harness energy from a wide range of sources, including: human walking (through devices sewn into clothing) [434], wind gusts (with oscillating flags) [435], and sea waves (using floating buoys) [436]. Many different working modes, material combinations and geometries have now been tested, with different TENGs showing better performance for various forms of mechanical energy [433].

Despite advancements in controlling charge output, the optimization of TENGs has been largely heuristic due to the lack of a solid theoretical foundations. Modelling the generation of triboelectric (nano-)currents, both in solid-solid and solid-liquid interfaces, is still based on the adoption of experimentally-found tribocharges surface densities. This is mainly due to the complexity of the phenomenon under investigation
and its dependence on the system being studied (solid-solid or solid-liquid), the modality of generation of charges (contact, sliding, wetting), the materials at play (and their position in the triboelectric series [437]), and the characteristics of the surfaces and their evolution, to name a few. Furthermore, the mechanisms governing charge generation, even in the simplest contact-separation modality whereby consecutively pressing and releasing two optimally paired surfaces produces charges and alternating currents (see Figure 15(a)), require fundamental research and models and theories developed at the smallest scales, while the influence of topography and materials microstructural properties require mesoscale models; the information obtained at the smaller scales must then be used as inputs to generate accurate models at the continuum scales aimed at predicting the behaviour and optimising the design of TENGs.

While significant progress has been made in some of the areas identified above, they have been system-dependent and sporadic. At the smallest of scales, researchers have attempted to explain why phonon-assisted tunneling plays a role in the generation of electromotive forces. The fundamental quantum mechanical mechanism that leads to an electromotive force between all co-moving bodies, suggesting that this is the origin of triboelectric charging, potentially paves the way for the development of a microphysical theory to predict the tribo-charges displacement [439]. Other fundamental studies have suggested a mechanism that attributes triboelectricity to the flexoelectric effect, which is present in all materials. [440]. This study showed that flexoelectric potential differences induced by inhomogeneous strains at nanoscale asperities drive tribocharge separation, which has opened new routes to investigate the triboelectric effect and extensions [441]. These, combined with other flexoelectricity studies (e.g. [442, 443]), have the potential to impact much of the research in this area and its exploitation. These theoretical models, which also include the thermoelectric effects [444], have explained some aspects but have not yet provided predictive capabilities for scalable engineering solutions, hampering the development of practical devices.

Even in the case of polymers, whose triboelectric behavior has been extensively studied, there is ongoing debate about whether the mechanism involves ionic transfer, material transfer, or electronic transfer. This uncertainty also applies to polymers like polytetrafluoroethylene (PTFE), known for its high triboelectric output. Recent *ab initio* investigations have shown that PTFE's triboelectricity is closely related to its electronic structure and defluorination [445]. As shown in Figure 15(b-c), this can, in the presence of negligible backflow and stuck charges, be directly applied to explain macroscopic behaviour and to demonstrate how defluorination enhances the triboelectric output of a TENG by an order of magnitude [438].

Looking at larger scales models, much effort has been recently devoted to the incorporation of surface roughness at the micro- and meso-scales. In particular, a unified model (for dielectric-to-dielectric TENGs) that adds consideration of surface roughness to the established distance-dependent electric field model has been recently developed and further extended [446, 447]. The model is applicable from first touch to nearly complete contact, provided that deformation remains elastic, and has obvious advantages compared to load-independent approaches. Incorporating the effect of surface roughness has also been shown to provide better prediction of TENGs performance using multiscale experiments coupled to theoretical analysis [448].

Research has also boomed when considering liquid-solid contact electrification and the capacity to collect energy from water-based TENGs. In this case, charge generation is associated with the formation of the electric double layer, with a multi-step process recently rationalised by looking at the electronic structure and electron transfer mechanisms at the smallest scales [449, 450, 451]. Many other models have been investigating charge transfer at the larger scales, with the majority of efforts focusing on surface design to maximise hydrophobicity to optimise charge transfer. A recent contribution has also looked at the intrinsic links between wetting (and de-wetting) dynamics and triboelectrification; using multiscale modelling tools and rigorous analyses, the authors have been able to link electric outputs with the surface structure and explain the mechanisms that govern the performance of water-based TENGs [452]. This type of modelling opens important routes for surface and materials design for practical applications in this context and the use of TENGs as sensor for wettability and contamination.

The generation of tribocharges, as a consequence of dynamic contact, can also be adopted to generate localized nanocurrents suitable to precision medicine applications. Indeed, triboelectrification has been used to fabricate neuromorphic tactile system, where a biological mechanoreceptor was mimicked by using a TENG comprising a self-powered pressure sensor and a bistable resistor neuron [453], as well as implemented to control cell regeneration and proliferation, or to provide bacteria suppression [454]. Electrical fields and
<figure><img src="image_16.png" /><figcaption>Figure 15: Schematics of working principle and possible governing mechanisms of contact-separation (CS) mode TENGs. (a) Diagram of the pushing/releasing mechanism to generate the triboelectric signal. The proposed working principles of a PTFE/Au CS TENG are also depicted: during approach (b) the electrostatic barrier lowers and enables the forward flow of electron; during the contact stage (c) charging is governed by the competition between the larger electron transfer prompted by an increasing defluorination and a backflow generated by a smaller interface separation. Images adapted from [438].</figcaption></figure>
electrical signaling are indeed essential in cellular processes. Triboelectric stimulation can thus induce reprogramming and proliferation of fibroblasts for wound skin healing, or it can enhance the proliferation and neural differentiation of mesenchymal cells, to cite a few [454]. Electrical stimulation is also one of the non-drug methods used to inhibit and control the growth of microorganisms, such as under soft tissue [455]. The above mentioned applications of triboelectricity-enhanced tissue interactions would strongly benefit from a theoretical understanding of the underpinning physics, currently lacking.

Another area of interest in the context of the exploitation of the triboelectric effect are tribotronic phototransistors. Researchers have shown that the coupling of a MoS2 phototransistor and a TENG in sliding mode can be used to improve the photoresponsivity of the device for photodetection [456]. Other efforts have more recently been focusing on the use of TENGs to produce flexible and self-powered photodetectors [457], large-area tribotonic transistors [458], and multibit tribotronic nonvolatile memory for intelligent instrumentation and self-powered wearable devices [459]. Such applications would certainly benefit from more in-depth research and advanced modelling efforts to unravel the coupling mechanisms and optimise these devices.

Turning now to the explicit use of electric fields (EFs) to control friction, the emerging trends embracing electrification and the acceleration of the adoption of electric vehicles has propelled advancements in lubricant technology for new operational environments under electric fields. Research has been conducted over decades on the effect of EFs generated by the EVs on lubricant performance, with a field referred to "Triboelectrochemistry" [460]. Research has particularly blossomed in this area, recognising the considerable opportunities for using applied EFs both to promote desirable and to suppress unwanted lubricant interactions with rubbing surfaces. The quest for better performance driven by nanoscale phenomena has led to various attempts to investigate electrotunable friction and lubrication using molecular dynamics simulations, both using canonical aqueous systems [461] and ionic liquids [462, 463]. Sustained efforts have been made to develop consistent methodologies to explain the key mechanisms driving the surface response (both in terms of reaction and wall slip) and identifying key challenges, such as preventing electrical shortcuts and developing optimal mixtures of the most promising compounds with organic solvent. In this field there is a clear need to link the research at the molecular scales with larger scale simulations (including the effect of surface roughness and changes at the fluid-solid interfaces) and explore the effect of alternating voltage in evolving tribological interactions to maximise the impact in terms of friction control and wear reduction, so that new compounds and formulations can be promptly used to promote the energy transition.

### 5.3. Batteries

The performance of all-solid-state batteries is critically influenced by the dynamic behavior of active interfaces, such as the contact between the lithium metal electrode and the solid-state electrolyte. The evolution of this interface during charge/discharge cycling involves multiple mechanical and chemical processes, spanning from the macroscale of the contact down to the nanoscale. These processes—including chemical, electrochemical, and chemo-mechanical phenomena [464]—depend strongly on battery operating conditions, such as current density and stack pressure [465], and quantitatively govern ionic charge transfer at the interface [464]. Solid-state electrolytes (SSEs) are generally classified into two main categories: i) inorganic solids and ii) solid polymer electrolytes, each associated with distinct multiscale mechanisms of ionic charge transfer [466]. The origin of interfacial phenomena in such systems is primarily determined by two interrelated factors: (i) the chemical and electrochemical compatibility between the electrodes and SSEs, and (ii) the stability of intimate (mechanical) contact between the mating surfaces.

A poor wetting between inorganic SSEs and molten lithium metal, as well as electrode surface roughening induced by lithium stripping/plating cycles, drastically affects ionic charge transfer due to the reduced (active) contact area where such transfer occurs [467]. More generally, the contact dynamics at the electrode/SSE interface play a key role in all-solid-state batteries, as achieving stable and nearly full intimate contact between active interfaces is essential. At the cathode/SSE interface, volumetric changes occur during (de)lithiation of the cathode, leading to cyclic and long-term evolution of surface micro- to nano-roughness and electro-mechanical characteristics. While similar effects may occur at the anode/SSE interface, the mechanical complexity is further exacerbated on the cathode side because conventional cathodes are composites of active materials, conductive additives, binders, and interface coatings. This composite nature results in
generally anisotropic and inhomogeneous bulk electromechanical properties [464]. For example, the role of binder rheology in the contact mechanics of cathode active particles has been investigated under the assumption of smooth contacting surfaces. In particular, plasticity has been shown to play a key role during the calendering process of the active cathode layer, strongly influencing its overall mechanical properties, such as differing unloading stiffness behaviors in compression and tension [468].

At the lithium metal/SSE interface, the anode can shift by several micrometers during charge/discharge cycling [467]. Maintaining contact under such interface roughening requires a relatively large effective work of adhesion. In practice, however, the effective adhesion between lithium anodes and SSEs is usually low. This arises from the large volume changes of the lithium metal during cycling—which increase the elastic strain energy at the interface and act effectively as a repulsive force—as well as from poor interfacial wetting. The resulting reduction in true contact area increases the interfacial impedance and, consequently, the local current density. Together with other critical factors (such as the SSE elastic modulus, the viscoplasticity of lithium metal, and compressive stresses at the interface), these conditions promote dendrite growth through the SSE, which can ultimately lead to short circuits [467].

In rough contacts, the true contact area can often be increased by applying a relatively soft adhesive coating, which lowers the elastic strain energy at the interface and raises the effective work of adhesion. By analogy, poor electrode/SSE contact can be mitigated using lithophilic polymer interlayers (e.g., polyethylene oxide [464]) [466]. However, thermal expansion mismatch typically induces thermal strains during cycling, which can cause warping, peeling, or cracking of the electrolyte [466]. In addition, even at open-circuit potentials, side reactions may occur and form passivating interphase layers at the electrode/SSE interface, whose thickness influences the mechanical contact area [464].

The contact pressure applied between the electrode and SSE during stacking is a key factor, as predicted by contact mechanics. Experiments on Li/LLZO interfaces [465] show that the true contact area (measured, for example, by acoustic impedance) depends nonlinearly on pressure during cycling. A threshold stack pressure is required to mitigate continuous void formation, although it does not completely suppress microstructural evolution [465]. More recently, ceramic inclusions in polymer interlayers of solid composite electrolytes (SCEs) have been shown to enhance battery performance [186]. This improvement was attributed to the evolving contact area, influenced by polymer viscoelasticity and dendrite growth at the anode.

Theoretical understanding of the interface science occurring in electrode/SSE interactions is typically achieved on two different length scales of observation. On one side, first principle calculations have been widely used to study the chemical and electrochemical stabilities of electrode/SSE interfaces, however a continuum description of such interactions is still limited [467]. Equally important, modeling approaches, from electronic and atomic scale to phenomenological (continuum) models need to be coordinated with each other in hybrid formulations in order to make quantitative predictions of the contact [469]. As an example, the fundamental physical, chemical and electrochemical properties (such as electron tunneling, ionic conductivity, adhesion, etc.) predicted from *ab-initio* models could be incorporated in large-scale homogenized formulation of the contact mechanics to provide a more predictive macroscopic picture of electrode/SSE contact and its impact on the overall battery performance.

At nanometer scale observation, the interface ionic resistance can be predicted with the adoption of the interface potential map model [470], based on bulk, surface and defect properties determined with density functional theory (DFT), even taking into account changes in the cathode with the actual state of charge. A high ionic resistance at the interface is indeed considered the battery technology bottleneck, due to its detrimental impact on the power performance, Coulombic efficiency, and short cycling life [470]. On the continuum level, contact mechanics has been adopted to investigate the interface impedance stability when a given contact pressure is applied to the stack, with a deterministic [471] or multiscale [472] approach, showing that lithium diffusivity is actually enhanced by mechanical stresses. The diffusion of lithium within the SSE has also been analyzed by including the chemical potential gradient as driving source (of diffusion) [437]. Interestingly, a marked asymmetry is found in the interface evolution characteristics during plating/stripping as a consequence of the reaction heterogeneity, mostly due to the non-uniform pressure and temperature distributions [473].
The closed-form frequency response function of a chemoelastic half plane has been only recently derived [474], and applied to the cylinder contact geometry. The theory could be adopted as well to build a chemoelastic contact mechanics based on Persson's formalism. Nevertheless, the contact interface in a generic electrode/SSE interaction might be far from being described by the typical roughness small slope assumption, such as during the occurrence of dendrites on the anode/SSE side. Interestingly, also the stripping process can be characterized by a random void accumulation on the anode [475], however, with an anodic void accumulation resembling the bubble production process in a liquid phase [475]. As such, the evolution of a Li anode void morphology has been modelled by the nucleation and growth theory [475].

A significant challenge inherent in the aforementioned numerical approaches is the comprehensive experimental acquisition of both bulk and localized properties (e.g. at grain boundaries in SSE). The related data is needed for a robust parametrization and subsequent validation of pertinent models. Standard experimental methodologies encompass electrochemical characterization, advanced imaging techniques (X-ray computed tomography, neutron and magnetic resonance imaging, in-situ (S)TEM), spectroscopic analysis, mechanical and thermal property assessment, and full-scale battery performance evaluation [476].

From a tribological perspective, the biggest challenge [477] is the prediction of the anode/SSE contact mechanics during dissolution, whose fundamental picture is shown in the schematic of Fig. 16a for a Li metal/SSE interface [478, 476]. In particular, the bottom-left schematic illustrates sustained interfacial contact, which occurs when vacancies created by Li stripping diffuse more rapidly than Li⁺ ions migrate under the influence of the anodic current. Conversely, the bottom-right schematic demonstrates void growth. The latter arises when the anodic current density surpasses the rate of vacancy replenishment, either by diffusion or mechanical deformation, leading to the nucleation and inward propagation of voids within the Li bulk. Recent *in-operando* X-ray computed tomography investigations, illustrated in Figure 16b [478, 476], have enabled the temporal mapping of Li/SSE interfacial contact. Colored pixels delineate contact regions, while black pixels indicate a lack of contact (corresponding contact area variations are presented in 16c, yellow points). The contact area is shown to decrease as time increases, due to a lithiation-induced surface roughening, in agreement with the aforementioned fundamental picture.

The findings just discussed show that multiscale contact mechanics modeling is imperative for simulating the impact of microstructural and material property evolution on interface degradation. This degradation is intrinsically linked to three primary chemo-mechanical challenges: Dendritic growth originating from Li-metal anodes, structural instability within composite cathodes, and solid electrolyte degradation driven by the formation of unstable interphases. The standard multiscale methodology typically integrates DFT at the nanometer scale, often off-line coupled with various continuum electro-chemo-mechanical formulations, predominantly employing FEM and/or Phase-Field Models (PFM) [479].

A recent multiscale, three-dimensional, time-dependent contact model delineating the evolution of the Li/SSE interface under stack pressure incorporates critical factors such as the surface roughness of both Li and SSE, Li elasto-plasticity, Li creep, and the dynamics of Li metal plating/stripping [480]. The study of that model established a correlation between contact elasto-plasticity and the ratio of Li yield strength to surface roughness, highlighting its significant influence on the formation and progression of Li dendrites and voids. Notably, the contact mechanics problem was addressed using fast Fourier transform-based numerical methods, a technique commonly employed in tribological research. The Butler-Volmer equation was utilized to quantify the effect of pressure on the stripping rate. A key finding was the observation of significant Li metal creep at the Li/SSE interface under pressure, even at ambient temperatures [480].

In Ref. [481], a transient, three-dimensional finite-element model of the evolution of a lithium anode due to stripping and plating is developed with *Sierra/Solid Mechanics*. The contact model incorporates both elastic and elastic-plastic deformation, alongside single-scale roughness. The interface between the separator and lithium is represented as a liquid electrolyte in this case. Charge and species transport within all materials are coupled at the lithium/separator and lithium/electrolyte surfaces using a modified Butler-Volmer equation. Evaluation of the galvanostatic, steady-state electrochemical response across the anode-separator gap revealed a strong dependence of interfacial reaction kinetics on applied pressure and assumed lithium yield strength.

A more sophisticated interface model, detailed in Ref. [482], integrates multiscale phenomena—specifically interface interactions, vacancy hopping, and plastic deformation—through a hierarchical hybrid workflow.
<figure><img src="image_17.png" /><figcaption>Figure 16: a), Schematic of the anode/SSE contact mechanics during dissolution, for a Li metal/SSE interface. b) In operando X-ray computed tomography results in terms of Li/SSE interfacial contact [478]. Colored pixels delineate contact regions, while black pixels indicate a lack of contact. c), contact area variations corresponding to b), with yellow points). Image adapted from [476].</figcaption></figure>