Is Ultralow Friction on Graphite Sustainable in Contaminated Environments? (2024)

1 Introduction

Structural lubricity¹ refers to an ultra-low friction state that arises from the systematic annihilation of interfacial lateral forces due to lattice mismatch.In this context, interfacial atoms can be treated as the smallest contact units, with their energy states, driven by thermal motion, distributed stochastically across a potential energy landscape.Stokesian damping², which describes a linear increase in friction with sliding velocity, results from small perturbations of these interfacial atoms.This leads to primary energy dissipation as the atoms vibrate around their equilibrium positions within the lattice, without engaging in long-range elastic or quasi-discontinuous motion under shear forces.A scaling argument, based on the law of large numbers, dictates a sub-linear relationship between friction and contact area^{3, 4}: $F\propto A^{n}$, where the exponent $n$ is capped at 0.5 and varies depending on the degree of lattice registry and geometric nature of the contact, including contact lines^{1, 5, 6}.While this low-friction state is highly desirable, it can become unstable when the elastic correlation length is exceeded^{7, 8}.Beyond this threshold, structural defects, such as dislocations, can cause the sliding surfaces to interlock, significantly increasing shear stresses up to the Peierls stress limit⁷.

Under ambient conditions, airborne molecules such as water and short alkanes tend to adsorb onto solid surfaces, forming an orderly arranged monolayer that aligns with the underlying solid lattice⁹.With prolonged exposure, this monolayer can develop into a nanometer-thick multilayer film¹⁰, as indicated by an oscillatory density profile resulting from the liquid’s wavenumber-dependent compressibility¹¹.The anisotropic behavior of confined liquids, characterized by high in-plane molecular ordering, enables a load-bearing ability to a certain degree while still allowing molecular diffusion^{11, 12, 13}.Interestingly, the presence of such an adsorption layer does not necessarily impede superlubricity^{14, 15, 16}, with friction only mildly increasing alongside observations of rejuvenation, aging, and friction switches¹⁶.This suggests that earlier nano-manipulation experiments may have overlooked the existence of these layers, as contamination is unlikely to be fully eliminated, even under nominal ultra-high vacuum conditions.Unlike the constant velocity gradient typical of Couette flow, where molecules must overcome energy barriers due to steric hindrances, the shear plane in the boundary lubrication regime is located at heterojunctions with the lowest shear strength, exhibiting solid-like behavior.Consequently, when sliding occurs over such a viscoelastic medium, the frictional dependence on sliding speed is expected to be sub-linear¹⁷.

Building on our previous studies^{6, 18}, we focus on highly oriented pyrolytic graphite (HOPG) due to its superior performance as a solid lubricant.The aim is to examine the extent to which superlubricity is maintained in the presence of linear alkanes, specifically $n$-hexadecane (HEX), as a representative adsorption contaminant.To replicate the experimentally observed morphology of parallel stripes of HEX, with their longitudinal axes aligned parallel to the zigzag direction of HOPG, we reparameterized the interatomic interaction between HOPG and HEX based on density-functional theory (DFT) calculations.To validate this new parameter set, we calculated binding and desorption energies and compared them with experimental measurements.We then investigate how topography influences molecular in-plane ordering during liquid confinement.Finally, we study energy dissipation both in area-filling buried interfaces and at contact lines within the boundary lubrication regime.

2 Methodology

Molecular desorption, liquid confinement, and boundary shearing are investigated using molecular dynamics (MD) simulations, employing model systems that incorporate Au(111) and graphite (0001) (Gr) as solid surfaces, with $n$-hexadecane (HEX) as the adsorbent.To ensure clarity, detailed geometric information and operational conditions for each model system are provided in the relevant sections.The interactions of Au, Gr, and HEX themselves are described using the EAM¹⁹, AIREBO²⁰ potentials, and the L-OPLS^{21, 22} force-field, respectively.For cross-interactions, Morse²³ and Lennard-Jones²⁴ potentials are utilized to describe Au-Gr and Au-HEX interactions, respectively.The interaction between Gr and HEX is reparameterized, with further details provided in the Results section.Unless stated otherwise, the system temperature is maintained at 300 K using a Langevin thermostat with a damping factor of 0.1 ps.The thermostat is applied in all directions except during sliding simulations, where it is restricted to the $y$-axis, perpendicular to the sliding direction.The simulation timestep is set to 1 fs for all cases.All MD simulations are carried out using the open-source code LAMMPS²⁵.To determine the optimal MD data output frequency, a time autocorrelation function (ACF) analysis is performed, with mean values and standard errors calculated based on uncorrelated data sets extracted.

3 Results and discussion

3.1 Interfacial Potential Development

The interfacial interaction between graphite (Gr) and $n$-hexadecane (HEX) is described using the Buckingham potential²⁶, formulated as:

$$E_{\rm int}=A{\rm e}^{-r/\rho}-\frac{C}{r^{6}},$$

(1)

where $A$, $\rho$, and $C$ are fitting parameters optimized through force matching (FM).This potential provides a more precise representation of short-range repulsive interactions with an exponential form, derived from the interpenetration of closed electron shells, making it more suitable for modeling boundary shearing^{27, 28}.Reparameterizing these potential parameters, rather than relying on mixing rules^{29, 30} for the non-bonded Lennard-Jones (LJ) parameters from the OPLS force field, is crucial because the latter, designed for bulk phase properties, fails to accurately capture the interfacial interactions, leading to a significant underestimation of the energy corrugation barrier.Two sets of parameters for ${\rm C_{Gr}}-{\rm C_{HEX}}$ and ${\rm C_{Gr}}-{\rm H_{HEX}}$ were determined by reproducing DFT-calculated static molecular forces ($f$) and binding energy ($\Delta E$, defined in Eq.4) using a single-molecule-on-graphene model system.Sixteen configurations ($N_{\rm conf}=16$) were considered, encompassing a range of both in-plane and out-of-plane relative positions between the two substances.

First-principles DFT calculations were performed using the Gaussian and plane-waves methodology^{31, 32} within the Quickstep module of the open-source software CP2K^{33, 34}.A plane-wave cutoff of 900 Hartree was applied to achieve a force convergence tolerance of 10^-5 atomic units (a.u.).$\kappa$-point sampling was employed along with periodic boundary conditions in all three dimensions.To minimize empirical influences, a van der Waals density functional (vdW-DF)³⁵ was used, combining exchange at the generalized-gradient approximation (GGA) level, correlation at the local density approximation (LDA) level, and a non-local correction for van der Waals interactions.Carbon and hydrogen atoms were described using a triple-$\zeta$ Gaussian basis set³⁶, paired with Goedecker-Teter-Hutter (GTH) pseudopotentials^{37, 38}.The results closely align with those obtained using other widely-used methods, such as Perdew-Burke-Ernzerhof (PBE)³⁹ with Grimme’s DFT-D3⁴⁰ corrections for dispersion interactions.

In FM, the discrepancy between DFT and force-field (FF) calculations was quantified using a penalty function $\chi^{2}$, defined as:

$$\chi^{2}=\frac{1}{D}\frac{\sum_{\alpha=1}^{D}\omega_{\alpha}\chi_{\alpha}^{2}}$$

(2)

where $D~{}(D=4)$ represents the dimensionality including three force components plus energy, and $\omega$ (with $\omega_{f_{x}}=\omega_{f_{y}}=5\times\omega_{f_{z}}=0.1\times\omega_{\Delta E}$) denotes the weighing factors determined based on their respective value spans.The term $\chi_{\alpha}^{2}$ is defined as:

$$\chi_{\alpha}^{2}=\frac{1}{N_{\rm conf}}\sum_{i=1}^{N_{\rm conf}}\frac{(\psi_{$$

(3)

where $\Delta=\underset{i,\alpha}{\rm max}(|\psi_{i,\alpha}^{\rm DFT}|)$ serves as a denominator adjustment to mitigate the impact of small values on fitting stability, and $\psi=\{f,\Delta E$} represents the physical properties being compared.The minimization of $\chi^{2}$ was executed using Monte Carlo-based simulated annealing⁴¹.The fitting parameters were accepted upon meeting the convergence criterion of $\chi^{2}\leq 10^{-3}$, where $\chi^{2}=0$ indicates perfect agreement between FF and DFT calculations.

Figure1 compares the results from DFT and FF calculations, showing that the newly parameterized Buckingham potential (red solid circles) outperforms others in accurately reproducing both molecular forces and energies.The corresponding potential parameters are listed in Table1.Notably, the out-of-plane forces ($f_{z}$) estimated with the mixed OPLS-LJ potential were significantly underestimated (blue open circles in Fig.1c), especially in the repulsive regime.As shown in Fig.1d and supported by our previous study⁴², the normal force component influences the energy barrier, thereby affecting lateral motion.Predictions based on the AIREBO potential (black open circles) are notably less accurate in this context.

$i-j$	$A_{ij}$ (eV)	$\rho_{ij}$ (Å)	$C_{ij}$ (eV$\cdot$Å6)
$\rm C_{Gr}-C_{HEX}$	1.75	0.09	57.13
$\rm C_{Gr}-H_{HEX}$	132.66	0.34	0.00

Reproducing collective properties like energy differences and molecular forces is relatively straightforward compared to reproducing per-atom forces, which require much complex functional forms that we were unable to identify.Specifically, we recognize our limitations in developing bond-order many-body potentials, including parameterizing the embedded-atom method (EAM)⁴³ and Steele⁴⁴ potentials, to account for possible many-body effects resulting from CH/$\pi$ interactions⁴⁵ that induce electron density changes upon approach.At this stage, we disregard the electrostatic interaction between Gr and HEX, assuming the graphene sheet is electronically neutral.We did not pursue machine-learning techniques, despite their promise, as the results might not align with established physical principles.Although imperfections exist, our approach achieves the initial objective of developing simple yet accurate means for modeling interfacial dynamics on a large scale with minimal computational cost.

3.2 Physisorption of $n$-hexadecane on Graphene

To validate the fitting parameters, we compared the binding energies ($\Delta E$) required to dissociate a single HEX molecule from a Gr substrate with experimental data^{46, 47}.These energies were determined based on two equilibrium energy states, formulated as:

$$\Delta E=\langle E_{\rm bound}\rangle-\langle E_{\rm sub}+E_{\rm HEX}^{\rm des}\rangle$$

(4)

Here, $E_{\rm bound}$ represents the potential energy of the system when the target molecule is adsorbed on Gr, while $E_{\rm sub}$ and $E_{\rm HEX}^{\rm des}$ correspond to the potential energies of the substrate (including any remaining HEX molecules) and the desorbed molecule in isolation, respectively.The notation $\langle\rangle$ denotes a time average taken over a minimum duration of 10 ns.The coverage density ($\it\Gamma$) is determined by calculating the total in-plane occupation of HEX molecules, denoted as $mA_{\rm HEX}^{\rm single}$, where the projected area of a single HEX molecule ($A_{\rm HEX}^{\rm single}$) is approximately 105 Ų⁴⁸.Here, $m$ represents the total number of molecules, which ranges from 1 to 30, corresponding to a $\it\Gamma$ varying from 0.03 to 0.77 monolayers (ML).

The MD-predicted binding energy for a single-molecule system ($m=1$) is 1.08 eV, closely matching the value of 1.11 eV obtained from DFT calculations⁴⁹, where the influence of thermal effects appears to be negligible.As additional molecules aggregate in close proximity ($m>1$), they arrange themselves in an orderly fashion, forming stripes with an average intermolecular spacing ($\Delta l$) of approximately 4.5Å, as examples illustrated in Fig.2(a) and (b).These adsorbed molecules adopt straight, all-$trans$ configurations, favoring an orientation along the zigzag direction of the graphite substrate, which aligns with experimental observations⁵⁰.Such a structural arrangement is energetically favorable, given that the periodicity along the graphene armchair direction is 4.26 Å, in close agreement with $\Delta l$.

Histograms showing the molecular orientation during dynamic equilibrium for two representative systems are presented in Fig.2(c), where the orientation vectors are determined as the end-to-end vector of each molecule along the longitudinal direction.An increase in binding energy of approximately 0.5 eV is observed as the number of molecules $m$ increases from one to numerous, as shown in Fig.3.Beyond this point, $\Delta E$ remains constant as $\it\Gamma$ further increases, regardless of the initial desorption position of the target molecule.This suggests that the binding energies obtained here correspond to the desorption of a molecule from the domain edges, which are preferred sites due to the low desorption energy barrier caused by structural asymmetry.

The binding energies presented in Fig.3 (red solid circles) are consistent with the MD-modeled energy barriers obtained from continuous pulling of a single molecule from an adsorption state to an isolated state along the surface normal while maintaining a constant center-of-mass (COM) displacement rate of 0.05 m/s.During this process, no reattachment of detached fragments is observed, owing to the high persistence length of these short-chain molecules.The correlation between free energy and the number of detached fragments suggests minimal configurational isomerism from $trans$ to $gauche$ conformations.These values are also in agreement with the desorption barrier ($\Delta E_{\rm des}$) extracted from temperature-programmed desorption (TPD) experiments (represented by the black dashed line in Fig.3), where the linearly interpolated peak desorption temperature ($T_{p}$) is approximately 308 K^{46, 47}, close to the 300 K used in MD simulations.Our simulation results support the experimentally observed first-order desorption process, where the desorption barrier is independent of coverage density, indicating relatively weak and short-ranged intermolecular interactions among the molecules.According to DFT calculations, the presence of oligomers on graphene could cause deviations in the electronic structure of graphene from the Fermi level observed in pristine graphene, resulting in an energy gap of up to 6 meV⁴⁹.

3.3 Nanoscale Confinement and Shearing

3.3.1 Solvation Forces

When liquids are confined within nanoscale slits, they often form layered structures^{51, 52} as molecules self-organize with significant in-plane order due to the perturbative effect of solid surfaces.A solid surface can induce in-plane molecular ordering that extends through several molecular layers, manifesting in oscillatory density profiles that decay with distance from the solid-liquid interface.The in-plane ordering of HEX on Gr is notably more pronounced than on a Au(111) surface, as shown in Fig.4.Specifically, the peak density of the HEX layer in immanent contact with Gr is nearly 2.5 times that on a Au(111) surface and 6.2 times that in a bulk phase.This structured, densely packed, anisotropic arrangement of molecules across neighbouring layers can withstand substantial normal loads on the order of 100 MPa without necessitating a phase transformation from liquid to solid¹¹.Additionally, the prolonged relaxation time observed in confined liquids, which is inversely related to shear viscosity, leads to noticeable hysteresis in stress relaxation¹², underscoring the complex dynamics at play in nano-confined liquid systems.

As a tip progressively approaches a substrate surface, the layer-by-layer squeezing out caused by in-plane molecular ordering leads to oscillatory solvation forces, with their magnitude influenced by the loading rate.As shown in Fig.5, the peak normal stresses ($\sigma_{zz}^{\rm p}$) during nominally steady-state compression ($v_{z}=0.2$ m/s) increase exponentially as the confining distance ($d$) decreases, with a wavelength comparable to the interlayer distance.The slightly negative stresses at the valleys correspond to the transition where molecular layers are expelled from $n$ to $n-1$.Our results suggest that pushing the last layer of molecules (at $d<0.9$ nm) outside the contact interface to achieve solid-solid contact requires substantial effort and may even be impossible, especially when the tips are blunt, resulting in a large aspect ratio.The heterogeneity of liquids under nanoscale confinement and the resulting stress anisotropy arise naturally from the liquid’s wavelength-dependent compressibility, which aligns with bulk-liquid density autocorrelation functions (ACF)^{53, 54}, rather than the specific nature of the solid surfaces.However, detailed molecular arrangements are associated with surface topography; thus, $\sigma_{zz}$ on Gr is notably higher due to the orientation-dependent higher Gr-HEX commensurability, which induces a higher corrugation barrier, making lateral displacement more difficult.

3.4 Shear Stress

Shear stresses ($\tau$), calculated as the ratio of lateral force ($F$) to the apparent contact area ($A$), are evaluated when a Au(111) slab slides against a HEX film situated atop a Gr substrate under an adhesive load.This loading condition is commonly employed in nanoisland manipulation experiments conducted in contact mode^{14, 42, 16}.With a monolayer of HEX ($\it\Gamma\approx$ 0.8), shear stresses exhibit a strong correlation with the lattice orientation of the confining solids.Specifically, as illustrated by the two leftmost blue bars in Fig.6, $\tau$ is approximately five times higher when the gold’s $[11\overline{2}]$ direction vector aligns parallel to graphene’s armchair direction (the sliding direction) compared to when its orthogonal vector, $[1\overline{1}0]$, aligns in the same direction.Note that the stresses reported in Fig.6 are from area-filling contacts, which exclude the energy dissipation at the contact lines⁶ or due to spontaneous cluster rotation¹⁸.

In the presence of a single monolayer of HEX molecules, the orientation-dependent shear stresses exhibit solid-like behavior, with the HEX monolayer effectively mimicking the graphene lattice.The notable lattice congruency between Au(111) and Gr(0001), as depicted in Fig.7b, facilitates the formation of Moiré patterns⁶, which, with extended superlattice wavelengths, promote longer-range instabilities at small lattice mismatches characterized by an elevated stress-strain gradient.Although the presence of a HEX monolayer amplifies $\tau$ by nearly a hundredfold compared to gold sliding directly against graphite⁶, it does not cause a breakdown of superlubricity.One reason is that this modeled large friction discrepancy would be significantly reduced at experimentally²²2Refers specifically to nanomanipulation experiments conducted with AFM realistic sliding speeds, which are typically nine orders of magnitude lower, due to different friction-velocity dependencies in the presence and absence of intermediate contaminants (discuss next).Additionally, the much larger contact area of nanoislands, on the order of 10${}^{3}\sim$ 10⁵ nm² in experiments, will mitigate the effects of energy dissipation at contact lines caused by structural discontinuities⁶ or molecular plowing⁵⁵.

In contrast, the orientation effect is not discernible in the Au-HEX-Au systems, as shown by the two rightmost blue bars in Fig.6.This can be attributed to HEX’s low binding energy to Au(111)⁵⁶, which fails to foster the same degree of molecular ordering observed on graphite, even when the two solid slabs are commensurate, as evidenced by the less pronounced density peaks in Fig.4.The orientation effect also diminishes as the HEX film thickness increases from one to six molecular layers, where viscous damping becomes increasingly dominant due to the decay of in-plane molecular ordering with distance from the solid surfaces.As a result, the reduced interlayer lattice congruency renders the solid’s lattice orientation insignificant, leading to comparable shear stresses (shown by the red bars in Fig.6) between systems with low and high commensurability.However, shear stresses arising from elastic instability continue to surpass those from viscous stresses under conditions equivalent to shearing a bulk liquid in a Couette flow¹⁸.In the latter case, although molecules align along the streaming direction according to their anisotropic radius of gyration¹⁸, their spatially-resolved density remains evenly distributed normal to the sliding plane rather than being oscillatory.

With a finite geometry, contributions from both inside and outside the contact are taken into account.Interestingly, as shown in Fig.8, shear stresses exhibit different dependencies on sliding speed depending on the coverage density.In the absence of contaminants (dry contact, $\it\Gamma$ = 0), the stresses (black circles) are the lowest, with damping being Coulombic, as indicated by the power-law exponent $n$ being nearly zero.When HOPG is lightly contaminated, meaning the island remains in direct contact with HOPG but with HEX molecules in close proximity ($\it\Gamma$ = 0.5), the shear stresses (blue triangles) are the highest, mainly due to the energy dissipated at the leading edge as HEX molecules are displaced, overcoming the accumulated corrugation barrier.At low speeds, molecules have sufficient time to refill the sliding-induced vacancies at the trailing edge, while at high speeds, they cannot quickly flow back to the sliding trace due to hysteresis, leading to a quasi-Stokesian dependency ($n=0.67$).In an extreme case, this can result in a drop in shear stress (at $v_{x}=50$ m/s) during reciprocal sliding when molecules cannot return in time at the leading edge.When a nearly full monolayer of contaminants forms such that the island floats on it, the shear stresses (red squares) remain low and exhibit a weak dependency on sliding speed, consistent with experimental observations⁵⁷.This behavior indicates that the monolayer is in a solid-like state, with closely packed and aligned molecules, persisting superlubricity or to say quasi-structural lubricity.

4 Conclusions

Under ambient conditions, the adsorption of airborne molecules on graphite (Gr) is inevitable, leading to the formation of orderly arranged molecular monolayers that impact interfacial mechanical and tribological properties.This study specifically examines $n$-hexadecane (HEX), a linear alkane, as a representative contaminant.To enable reliable molecular dynamics predictions, the interfacial interaction between HEX and Gr was firstly parameterized using the Buckingham formalism, based on DFT-calculated molecular forces and binding energies.This approach enhances predictions of force components, particularly in the normal direction, and thereby the corrugation energy barriers encountered during the lateral motions.The aim is to develop a simple yet accurate means for predicting interfacial behavior, while implicitly accounting for variations in electrostatic properties upon contact.The robustness of the developed potential parameters was next validated by reproducing the coverage density-dependent binding energy derived from temperature-programmed desorption experiments⁴⁷, with errors $\leq 5\%$.Furthermore, the study confirms that HEX molecules adopt a straight, all-$trans$ configuration aligned parallel to the graphite zigzag direction upon adsorption, consistent with experimental observations⁵⁰.

When a HEX film with multi-layer thickness forms on Gr, the pronounced layering-like structure, evidenced by the oscillatory density peaks along the surface normal, is significantly more distinct than on Au(111).This distinction arises from the high in-plane molecular ordering induced by the lattice similarity between the parallel-aligned HEX molecules and Gr(0001).Such structural anisotropy can lead to a discontinuous, layer-by-layer expulsion during steady-state compression, with the last monolayer of HEX being particularly resistant to removal.This monolayer, directly in contact with Gr, exhibits solid-like behavior despite its diffusive nature, as indicated by the orientation-dependent shear stresses observed during boundary shearing.This behavior contrasts with systems where the commensurability is low, such as on Au(111), where shear stresses are lower and unaffected by detailed surface topography or lattice orientation, and the shear plane locates typically at such heterojunctions.

The monolayer of adsorbents can serve as a protective film for substrate surfaces, where superlubricity does not necessarily break down when accounting for the length and speed scaling from the MD scope to practical engineering applications, i.e., from $\mathcal{O}(1~{}{\rm nm}^{2},1\sim 10~{}{\rm m/s})$ to $\mathcal{O}(10^{2}\sim 10^{4}~{}{\rm nm}^{2},10^{-1}\sim 10^{2}~{}{\rm\mu m/s})$.As the coverage density increases from $\it\Gamma$ = 0 to a full monolayer, the dependencies of shear stresses on sliding velocity transition non-monotonically from Coulomb to quasi-Stokesian, and then to quasi-Coulomb, highlighting the variation in energy dissipation channels within and outside the contact.Among these scenarios, displacing correlated molecules from the leading edge requires considerable energy, which increases nearly linearly with sliding speed up to a threshold, beyond which shear stress drops due to the lack of molecules reflow into the sliding trace during reciprocal sliding as a result of hysteresis.Conversely, when floating on a monolayer of HEX, shear stress increases only slightly, with comparable dissipation occurring both inside and outside the contact.As the thickness of the HEX film increases, the impact of lattice orientation diminishes due to the increasingly viscous nature of the film, shifting the lubrication scenario from the boundary lubrication (BL) to the elastohydrodynamic lubrication (EHL).

{acknowledgement}

We thank Martin Müser, Mehmet Baykara, Wengen Ouyang, and Wai Oo for useful discussions.This research was supported by the German Research Foundation (DFG) under grant number GA 3059/2-1.

References

• Müser 2004Müser,M.H. Structural lubricity: Role of dimension and symmetry. Europhysics Letters 2004, 66, 97
• Adelman and Doll 1976Adelman,S.A.; Doll,J.D. Generalized Langevin equation approach for atom/solid‐surface scattering: General formulation for classical scattering off harmonic solids. The Journal of Chemical Physics 1976, 64, 2375–2388
• Müser \latinetal. 2001Müser,M.H.; Wenning,L.; Robbins,M.O. Simple Microscopic Theory of Amontons’s Laws for Static Friction. Phys. Rev. Lett. 2001, 86, 1295–1298
• Dietzel \latinetal. 2013Dietzel,D.; Feldmann,M.; Schwarz,U.D.; Fuchs,H.; Schirmeisen,A. Scaling Laws of Structural Lubricity. Phys. Rev. Lett. 2013, 111, 235502
• deWijn 2012deWijn,A.S. (In)commensurability, scaling, and multiplicity of friction in nanocrystals and application to gold nanocrystals on graphite. Phys. Rev. B 2012, 86, 085429
• Gao and Müser 2022Gao,H.; Müser,M.H. Structural lubricity of physisorbed gold clusters on graphite and its breakdown: Role of boundary conditions and contact lines. Frontiers in Chemistry 2022, 10
• Sharp \latinetal. 2016Sharp,T.A.; Pastewka,L.; Robbins,M.O. Elasticity limits structural superlubricity in large contacts. Phys. Rev. B 2016, 93, 121402
• Monti and Robbins 2020Monti,J.M.; Robbins,M.O. Sliding Friction of Amorphous Asperities on Crystalline Substrates: Scaling with Contact Radius and Substrate Thickness. ACS Nano 2020, 14, 16997–17003, PMID: 33226231
• He \latinetal. 1999He,G.; Müser,M.H.; Robbins,M.O. Adsorbed Layers and the Origin of Static Friction. Science 1999, 284, 1650–1652
• Annett and Cross 2016Annett,J.; Cross,G. L.W. Self-assembly of graphene ribbons by spontaneous self-tearing and peeling from a substrate. Nature 2016, 535, 271–275
• Gao and Müser 2020Gao,H.; Müser,M.H. Why liquids can appear to solidify during squeeze-out – Even when they don’t. Journal of Colloid and Interface Science 2020, 562, 273–278
• Gao 2022Gao,H. History-Dependent Stress Relaxation of Liquids under High-Confinement: A Molecular Dynamics Study. Lubricants 2022, 10
• Huang \latinetal. 2023Huang,X.; Li,T.; Wang,J.; Xia,K.; Tan,Z.; Peng,D.; Xiang,X.; Liu,B.; Ma,M.; Zheng,Q. Robust microscale structural superlubricity between graphite and nanostructured surface. Nature Communications 2023, 14, 2931
• Cihan \latinetal. 2016Cihan,E.; İpek,S.; Durgun,E.; Baykara,M.Z. Structural lubricity under ambient conditions. Nature Communications 2016, 7, 12055
• Deng \latinetal. 2018Deng,H.; Ma,M.; Song,Y.; He,Q.; Zheng,Q. Structural superlubricity in graphite flakes assembled under ambient conditions. Nanoscale 2018, 10, 14314–14320
• Oo \latinetal. 2024Oo,W.H.; Gao,H.; Müser,M.H.; Baykara,M.Z. Structural Lubricity and Molecular Contamination: Rejuvenation, Aging, and Friction Switches. 2024; \urlhttps://arxiv.org/abs/2407.03360
• Müser 2020Müser,M.H. Shear Thinning in the Prandtl Model and Its Relation to Generalized Newtonian Fluids. Lubricants 2020, 8
• Gao and Müser 2024Gao,H.; Müser,M.H. On the Shear-Thinning of Alkanes. Tribology Letters 2024, 72, 16
• Zhou \latinetal. 2004Zhou,X.W.; Johnson,R.A.; Wadley,H. N.G. Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers. Phys. Rev. B 2004, 69, 144113
• Stuart \latinetal. 2000Stuart,S.J.; Tutein,A.B.; Harrison,J.A. A reactive potential for hydrocarbons with intermolecular interactions. The Journal of Chemical Physics 2000, 112, 6472–6486
• Price \latinetal. 2001Price,M. L.P.; Ostrovsky,D.; Jorgensen,W.L. Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the OPLS-AA force field. Journal of Computational Chemistry 2001, 22, 1340–1352
• Siu \latinetal. 2012Siu,S. W.I.; Pluhackova,K.; Böckmann,R.A. Optimization of the OPLS-AA Force Field for Long Hydrocarbons. Journal of Chemical Theory and Computation 2012, 8, 1459–1470, PMID: 26596756
• dela Rosa-Abad \latinetal. 2016dela Rosa-Abad,J.A.; Soldano,G.J.; Mejía-Rosales,S.J.; Mariscal,M.M. Immobilization of Au nanoparticles on graphite tunnels through nanocapillarity. RSC Adv. 2016, 6, 77195–77200
• Pu \latinetal. 2007Pu,Q.; Leng,Y.; Zhao,X.; Cummings,P.T. Molecular simulations of stretching gold nanowires in solvents. Nanotechnology 2007, 18, 424007
• Thompson \latinetal. 2022Thompson,A.P.; Aktulga,H.M.; Berger,R.; Bolintineanu,D.S.; Brown,W.M.; Crozier,P.S.; in ’t Veld,P.J.; Kohlmeyer,A.; Moore,S.G.; Nguyen,T.D.; Shan,R.; Stevens,M.J.; Tranchida,J.; Trott,C.; Plimpton,S.J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications 2022, 271, 108171
• Buckingham and Lennard-Jones 1938Buckingham,R.A.; Lennard-Jones,J.E. The classical equation of state of gaseous helium, neon and argon. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 1938, 168, 264–283
• Kong \latinetal. 2009Kong,L.-T.; Denniston,C.; Müser,M.H.; Qi,Y. Non-bonded force field for the interaction between metals and organic molecules: a case study of olefins on aluminum. Phys. Chem. Chem. Phys. 2009, 11, 10195–10203
• Müser 2022Müser,M.H. Improved cutoff functions for short-range potentials and the Wolf summation. Molecular Simulation 2022, 48, 1393–1401
• Lorentz 1881Lorentz,H.A. Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Annalen der Physik 1881, 248, 127–136
• Berthelot 1898Berthelot,D. Sur le mélange des gaz. Compt. Rendus 1898, 126, 15
• VandeVondele \latinetal. 2005VandeVondele,J.; Krack,M.; Mohamed,F.; Parrinello,M.; Chassaing,T.; Hutter,J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Computer Physics Communications 2005, 167, 103–128
• Lippert \latinetal. 1997Lippert,G.; Hutter,J.; Parrinello,M. A hybrid Gaussian and plane wave density functional scheme. Molecular Physics 1997, 92, 477–488
• Hutter \latinetal. 2014Hutter,J.; Iannuzzi,M.; Schiffmann,F.; VandeVondele,J. cp2k: atomistic simulations of condensed matter systems. WIREs Computational Molecular Science 2014, 4, 15–25
• Kühne \latinetal. 2020Kühne,T.D. \latinetal. CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. The Journal of Chemical Physics 2020, 152, 194103
• Dion \latinetal. 2004Dion,M.; Rydberg,H.; Schröder,E.; Langreth,D.C.; Lundqvist,B.I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401
• VandeVondele and Hutter 2007VandeVondele,J.; Hutter,J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. The Journal of Chemical Physics 2007, 127, 114105
• Goedecker \latinetal. 1996Goedecker,S.; Teter,M.; Hutter,J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710
• Hartwigsen \latinetal. 1998Hartwigsen,C.; Goedecker,S.; Hutter,J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. Phys. Rev. B 1998, 58, 3641–3662
• Perdew \latinetal. 1996Perdew,J.P.; Burke,K.; Ernzerhof,M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868
• Grimme \latinetal. 2010Grimme,S.; Antony,J.; Ehrlich,S.; Krieg,H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. The Journal of Chemical Physics 2010, 132, 154104
• Kirkpatrick \latinetal. 1983Kirkpatrick,S.; Gelatt,C.D.; Vecchi,M.P. Optimization by Simulated Annealing. Science 1983, 220, 671–680
• Oo \latinetal. 2023Oo,W.; Baykara,M.; Gao,H. A Computational Study of Cluster Dynamics in Structural Lubricity: Role of Cluster Rotation. Tribology Letters 2023, 71, 115
• Daw and Baskes 1984Daw,M.S.; Baskes,M.I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Phys. Rev. B 1984, 29, 6443–6453
• Steele 1973Steele,W.A. The physical interaction of gases with crystalline solids: I. Gas-solid energies and properties of isolated adsorbed atoms. Surface Science 1973, 36, 317–352
• Tsuzuki and Fujii 2008Tsuzuki,S.; Fujii,A. Nature and physical origin of CH/$\pi$ interaction: significant difference from conventional hydrogen bonds. Phys. Chem. Chem. Phys. 2008, 10, 2584–2594
• Paserba and Gellman 2001Paserba,K.R.; Gellman,A.J. Kinetics and Energetics of Oligomer Desorption from Surfaces. Phys. Rev. Lett. 2001, 86, 4338–4341
• Gellman and Paserba 2002Gellman,A.J.; Paserba,K.R. Kinetics and Mechanism of Oligomer Desorption from Surfaces: n-Alkanes on Graphite. The Journal of Physical Chemistry B 2002, 106, 13231–13241
• Londero \latinetal. 2012Londero,E.; Karlson,E.K.; Landahl,M.; Ostrovskii,D.; Rydberg,J.D.; Schröder,E. Desorption of n-alkanes from graphene: a van der Waals density functional study. Journal of Physics: Condensed Matter 2012, 24, 424212
• Kamiya and Okada 2013Kamiya,K.; Okada,S. Energetics and Electronic Structures of Alkanes and Polyethylene Adsorbed on Graphene. Japanese Journal of Applied Physics 2013, 52, 06GD10
• McGonigal \latinetal. 1990McGonigal,G.C.; Bernhardt,R.H.; Thomson,D.J. Imaging alkane layers at the liquid/graphite interface with the scanning tunneling microscope. Applied Physics Letters 1990, 57, 28–30
• Horn and Israelachvili 1981Horn,R.G.; Israelachvili,J.N. Direct measurement of structural forces between two surfaces in a nonpolar liquid. The Journal of Chemical Physics 1981, 75, 1400–1411
• Krass \latinetal. 2016Krass,M.-D.; Gosvami,N.N.; Carpick,R.W.; Müser,M.H.; Bennewitz,R. Dynamic shear force microscopy of viscosity in nanometer-confined hexadecane layers. Journal of Physics: Condensed Matter 2016, 28, 134004
• Fisher and Wiodm 1969Fisher,M.E.; Wiodm,B. Decay of Correlations in Linear Systems. The Journal of Chemical Physics 1969, 50, 3756–3772
• Nygård \latinetal. 2016Nygård,K.; Sarman,S.; Hyltegren,K.; Chodankar,S.; Perret,E.; Buitenhuis,J.; vander Veen,J.F.; Kjellander,R. Density Fluctuations of Hard-Sphere Fluids in Narrow Confinement. Phys. Rev. X 2016, 6, 011014
• Flater \latinetal. 2007Flater,E.E.; Ashurst,W.R.; Carpick,R.W. Nanotribology of Octadecyltrichlorosilane Monolayers and Silicon: Self-Mated versus Unmated Interfaces and Local Packing Density Effects. Langmuir 2007, 23, 9242–9252
• Wetterer \latinetal. 1998Wetterer,S.M.; Lavrich,D.J.; Cummings,T.; Bernasek,S.L.; Scoles,G. Energetics and Kinetics of the Physisorption of Hydrocarbons on Au(111). The Journal of Physical Chemistry B 1998, 102, 9266–9275
• Oo \latinetal. 2024Oo,W.H.; Ashby,P.D.; Baykara,M.Z. Structural Superlubricity at High Sliding Speeds under Ambient Conditions. 2024; \urlhttps://arxiv.org/abs/2407.06971