1 Introduction

Although zinc dialkyldithiophosphates (ZDDPs) have been used as antiwear additives for more than 80 years, they remain one of the most critical ingredients in engine lubricants [1]. When subjected to high temperatures and/or large stresses, ZDDPs decompose [2,3,4,5,6], and eventually form heterogeneous, patchy films, whose height increases with rubbing time [1, 7,8,9,10]. Not only the pad shapes turn out heterogeneous but also their mechanical properties: maximum values for elastic modulus and hardness are typically located on the most highly loaded regions of the films [11,12,13,14]. Despite significant progress, the antiwear-film-formation pathways are not fully established, and no consensus has been reached on the reasons for why the morphology and the mechanical properties of the tribo layer are so diverse [1, 6, 9, 15].

It seems clear that the film formation consists of several critical steps, which do not necessarily occur in the sequence mentioned in the following. In one step, ZDDP needs to decompose into its zinc phosphate-rich active products (ZnPs) [16,17,18] and the remaining sulphur as well as alkyl and aryl groups, whose content in (good) tribo-films tends to be relatively small [18, 19], although sulphur may certainly be critical to bind the films to metal surfaces [20, 21]. In another step, the ZnPs must become the dominant species on the surface of the growing film. It may be difficult to say with certainty, which of these two steps is first or if there is a well-defined sequence at all. However, it was argued that ZDDP adsorbs to (polar) surfaces before it decomposes [19, 22, 23]. In either case, the ZnPs must be attracted to the surfaces, e.g., via long-range Coulomb interactions and mirror charges induced in metals, or, to polar ionic terminations and/or to the highly polarizable preexisting ZnP film. In yet another step, initially disconnected ZnPs form networks, which are sufficiently stiff to protect the surfaces from counter faces but also compliant enough to be sacrificial under extreme rubbing conditions [7, 24, 25]. This step has traditionally been argued to occur through nucleophilic substitution [26, 27], which however, does not explain why mechanical properties of zinc phosphates films differ substantially between valleys and peaks [12,13,14]. Thus, the final aspect of film formation should be film strengthening, which is the topic of this work.

In the absence of asperity contact, it was shown that the rate-limiting process for ZnP formation is the ZDDP decomposition [27], which would be consistent with the observation that elevated temperature [28, 29] and non-contact shear stresses [6] speed up the film formation. It seems plausible that the decomposition remains the rate-limiting process in contact-driven film growth. However, another process may become rate limiting, if contact stresses further reduce the energy barrier to dissociation.

It is certainly established that thermal films made up of ZDDP decomposition products at temperatures of 150 °C and 200 °C are quite soft [7] and do not significantly protect surfaces from wearing [2]. This, besides the observation of stiff films on top of asperities and soft ones in valleys, supports the idea that the properties of zinc phosphate networks arise as pressure-hysteresis effect, in which large contact pressures induce a greater connectivity and thus enhanced stiffness of the ZnPs [15, 30]. This conjecture was based on simulations [30] that had revealed an irreversible densification of initially disconnected ZnPs near 6 GPa, accompanying an irreversible coordination change on zinc atoms. In this initial work, a reversible hybridization change on zinc atoms was observed at pressures near 17 GPa alongside with irreversible development of connectivity. We find this estimate substantially reduced in this work in agreement with experimental findings on similar compounds [31].

In a series of high-pressure experiments, Shakhvorostov et al. [31,32,33,34] did not only provide additional evidence of but also important further refinements to the pressure-assisted ZnP-network-formation (PANF) conjecture. Compression of zinc \(\alpha\)-orthophosphates in a diamond-anvil cell confirmed the prediction that zinc atoms abandon their tetrahedral coordination at hydrostatic pressures at around 6 GPa [32]. More importantly, Raman spectra revealed irreversible coordination changes on zinc atoms after a flat-punch uniaxial compression of ZDDP decomposition model compounds having been placed on a copper foil, which underwent severe plastic deformation in the process. A comparison [31] of X-ray diffraction (XRD) and infrared (IR) absorbance spectra of various model substances, similar to those used here and in Ref. [30], revealed similarity of thermal films to uncompressed model compounds and to be of similar stiffness as tribo-films located in the films between asperities, while films on top of asperities were found to correlate with ZnPs having been decompressed from large pressures. Last but not least, hydroxylation of ZnPs was found to be reversible: stiff tribofilms (90 GPa indentation modulus) had softened to 30 GPa indentation modulus after year-long exposure to humidity [34]. This softening was completely reversed through nanoindentation with an interfacial force microscope [34], thereby providing evidence that stresses may not only change the hybridization of Zn atoms in ZnP films but also their stoichiometry.

Calcium phosphates do not show the same complexity as zinc phosphates [32, 33], probably due to the absence of directed bonds allowing a competition of and switching between different hybridizations [6]. This would explain why calcium phosphates are not used to inhibit wear. Likewise, undecomposed ZDDP has no anti-wear functionality [16,17,18]. The absence of interesting pressure hysteresis in ZDDP [35] is therefore in agreement with the PANF conjecture, see also Ref. [36], but yet occasionally held as evidence against it [6].

Atomic force microscope (AFM) experiments provided the arguably most direct support for the PANF conjecture: the growth rate of ZnP films increased exponentially not only with temperature but also with the applied compressive stress [9]. In addition, ZnP films formed on aluminium [37], once it was sufficiently work hardened to accommodate the contact stresses needed to promote the film growth. This observation is in line with the stiffening of ZnP powder deposited on copper, which only occurred on top of severely plastically deformed spots during a large-scale/flat-punch indentation [32]. Further support for the emerging idea that stiffening enhances the growth of resilient ZnP films comes from new AFM studies revealing that films grow on hard MgAl alloys while no films form under usual rubbing conditions on elemental magnesium or aluminium crystals [38].

A drawback of most theoretical studies on the pressure-assisted network formation of ZnPs is that pressure is applied isotropically, while that of uniaxial indentation experiments is anisotropic. In addition the degree of stress anisotropy during the film formation is unknown and impractical to determine, since at best three out of six stress-tensor elements can be (crudely) estimated from experiment while rubbing, i.e., those carrying at least one times the z index This is why, confirmations of the PANF conjecture with [9, 37] or without [32] (intentionally added) shear stress must be taken with a grain of salt.

The effect that stress anisotropy has on structural/chemical changes in highly compressed matter is typically investigated for crystals [39,40,41,42,43,44,45,46]. Part of the reason for this may be that anisotropic stress in disordered media (i.e., glasses, which can be seen as complex liquids with extremely large viscosity) cannot be described in terms of a linear-response equilibrium theory relating stress and strain. A clear time-scale separation of structural relaxation and shear/compression would be required to achieve that. In this sense it is only possible to apply a true equilibrium stress anisotropy on crystals but not on glasses.

The overall trend that stress anisotropy has on phase transformations in crystals is that the hydrostatic pressures, at which the transformation is triggered upon compression, is reduced when the stress-tensor eigenvalues are not all equal [40,41,42,43,44,45,46], i.e., in the presence of shear. This effect was observed for the bcc to hcp transition in Fe [42, 43] and cubic diamond to \(\beta\)-tin structure in silicon [44, 45]. Gao et. al. [46] found a reduction of transition pressure from graphite to hexagonal diamond and to nanocrystalline cubic diamond phases by up to a factor of 100 in the presence of shear. The pressure at which \(\alpha\)-SiO\(_2\) becomes amorphous also turns out to be reduced by shear [40, 41]. A similar observation was made for \(\alpha\)-AlPO\(_4\) [47]. These observations therefore are valid not only for structural phase transitions between two crystalline phases, but also for amorphization transformations. The presence of shear during compression can even trigger the generation of metastable phases that are inaccessible without shear [46, 48, 49]. To make things more complicated, in some cases, such as for the transition from hexagonal to wurtzitic BN [50] and the transition between B1 and B2 phases of NaCl [51], non-hydrostatic compressions do not change the hydrostatic transition pressure to a measurable extent.

In this work we want to explore how stress anisotropy, i.e., the presence of shear stress, affects the hydrostatic pressures needed to promote the formation of stiff, simple ZnP networks under the assumption that the ZDDP decomposition had already taken place and an initial soft film was formed. This includes an analysis of how shear stress affects the structure and elastic properties of the stress-modified “films”. For this purpose, we perform density-functional theory (DFT) based molecular dynamics simulations using different deformation modes giving raise to varying stress anisotropies during compression. Model and methods are presented in Sect. 2, results in Sect. 3, and conclusions drawn in Sect. 4.

2 Model and Methods

To model the formation of a stiff ZDDP-derived network, we chose to start from ZnP-based decomposition products [16,17,18] for reasons of computational convenience. (Decomposition products that do not get incorporated inside the tribofilm would not disappear in a periodically repeated simulation cell.) To do so, we follow earlier theoretical studies on pressure-induced ZnP network formation using appropriate model molecules [30, 52]. Our periodically repeated, originally cubic simulation cell contains two triphosphate (\(\hbox {P}_3\hbox {O}_{10}\hbox {H}_5\)) and two zinc- phosphate molecules (\(\hbox {Zn}[\hbox {PO}_4\hbox {H}_2]_2\)), which is the same stoichiometry but twice the number of atoms compared to the reference simulations by Mosey et al. [30, 52]. Such molecules are meant to represent intermediate chain lengths similar to those observed experimentally for ZDDP [53]. We also studied compounds with higher Zn content, including crystalline orthophosphates to ensure that the trends reported in this work are robust. Since this work is mainly concerned with the question of how stress affects elastic properties of a film during growth, we do not address its binding to the substrate or other substrate-related effects. In addition, we decided to keep our focus on relatively simple systems not containing sulphur, as they usual allow the most general conclusions to be drawn, which may explain why Shakhvorostov et al. [31] succeeded in making a full-circle comparison of XRD and Raman spectra of model substances (DFT and experiment) and real tribofilms.

Four types of simulations were conducted in this study: energy minimizations at (i) constant stress and (ii) predefined strain tensor as well as finite-temperature simulations, which were conducted either with (iii) time-dependent pressure p during an isotropic compression of the cubic simulation cell, or (iv) at a predefined box-geometry or strain-tensor, which could change linearly with time . In a constant-p simulation, the volume of the simulation is allowed to vary but not its shape, while the shape is also treated as being dynamic at constant \(\sigma\).

The temperature was set to \(T = 600\) K in all finite-temperature simulations. This is a little more than 200 K above the operating temperature of engine oils. We chose this slightly increased value to speed-up chemical reactions, in an attempt to reduce the gap between our effective compression rates and those that occur during asperity collisions. At the same time, we remain well below estimates for flash temperatures [54], whose correctness we do not dare to judge.

In the following sections, we separate the description of further details on the simulation method itself and on the protocols used to generate the initial structure and to impose the deformation as well as on the observables, some of which are not frequently reported.

2.1 Simulation Details

All calculations were based on the density-functional theory (DFT) [55, 56] using Gaussian Plane Waves (GPW) [57] method as implemented in the CP2K package [58]. The Perdew–Burke–Ernzerhof exchange-correlation functional [59] was employed in combination with the empirical van der Waals corrections by Grimme [60]. We used Gaussian basis sets of double-\(\zeta\) quality [61] for all atoms in our system (H, O, P and Zn) in combination with Goedecker-Teter-Hutter pseudopotentials [62, 63]. Since the bandgap in ZnPs tends to be of order 5 eV and thus rather large, only the Gamma point needed to be sampled. The energy cutoff was increased from 120 Ry in the reference simulations to 400 Ry in our DFT-based molecular dynamics. It was further increased to 600 Ry in all static simulations including the energy minimizations.

The “canonical-sampling-through-velocity-rescaling” (CSVR) thermostat [64] was applied to atoms and the temperature was set to \(T=600\) K. If finite-temperature simulations were run at constant stress (or pressure), Nose-Hoover chains [65] were used as barostats.

In NpT-simulations, pressure was changed in steps of 1 GPa, and the system was given 10 ps to equilibrate at each pressure leading to an effective pressure rate of 0.1 GPa/ps. In strain-controlled simulations, the strain was changed in quanta of 0.02, which was followed by equilibration periods of 6 ps leading to an effective strain rate of approximately 3.3 GHz. This is roughly 12 orders of magnitude faster than strain rates in diamond-anvil-cell experiments, however, probably not too far away (on a logarithmic scale) from those that arise at 1 m/s sliding speeds between nanometer-sized asperities. The discrepancy between our simulations and the relatively mild elastohydrodynamic lubrication conditions in tribometer experiments [6] of 1 MHz is already only three decades. In general, we are not too concerned about the gap in time scales in the context of stress-induced transformation as long as our systems do not heat up, because comparisons between super-fast molecular simulations and super-slow diamond-anvil cell experiments are routinely done quite successfully, see, e.g., Refs. [32, 33, 66,67,68]. More details on strain-controlled simulations are given further below in Sect. 2.3.

2.2 Initial Configurations

The molecules were first placed at \(T = 600\) K in a relatively large cubic simulation cell with 30 Å long edges and equilibrated for 5 ps. Then, the system was gradually compressed by changing the linear size of the simulation cell in steps of 2 Å down to 12 Å. Each step involved another 5 ps equilibration. After this initial compression, an external pressure of 0.5 GPa was applied, and a relaxation of 20 ps was performed at constant pressure. One initial configuration is depicted in Fig. 1.

In the last 10 ps of the final relaxation run, configurations were dumped out each 2.5 ps yielding five different small-pressure configurations. Snapshots of such produced configurations reveal four separate molecules, interacting predominantly through hydrogen bonds, whose topology clearly changed between two subsequent configurations. Yet, the elastic tensors elements determined on those structures were very similar so that we decided to only keep the last configuration for the compression analysis.

Fig. 1
figure 1

Representative snapshots of the simulation box a at \(p=0.5\) GPa (initial structure), b at \(p=5\) GPa (isotropic compression), c at deformation \(\varepsilon =-0.24\) (uniaxial compression along z-axis), and d at deformation \(\varepsilon =-0.30\) (anisotropic, density-conserving deformation). Zn, P, O, and H atoms are drawn as gray, purple, red, and white balls, respectively. The colors of the polyhedra are consistent with that of their central atom. \(\varepsilon _{\alpha \alpha } < 0\) means that the deformation along \(\alpha\)-axis is compressive. \(\rho = {\text{const}}\) means that the particle density is conserved (see text for details)

Full size image

2.3 Imposing Deformation

The original structures were compressed: (i) isotropically, (ii) uniaxially at fixed areal density in the normal direction, i.e., the area of the simulation cell normal to the compression direction was kept constant, and (iii) uniaxially at fixed mass density. The elements of the strain tensor in the three compression modes obey

$$\begin{aligned} \epsilon _{\alpha \beta }^{\text{(i)}} = {\left\{ \begin{array}{ll} \varepsilon &{} {\text{ for }} \alpha = \beta \\ 0 &{} {\text{ else }} \end{array}\right. } \end{aligned}$$
(1a)

in the case of the isotropic compression,

$$\begin{aligned} \epsilon _{\alpha \beta }^{\text{(ii)}} = {\left\{ \begin{array}{ll} \varepsilon &{} {\text{ for }} \alpha = \beta = 3\\ 0 &{} {\text{ else }} \end{array}\right. } \end{aligned}$$
(1b)

for compression mode (ii), which will be referred to as simple uniaxial compression in the following, and

$$\begin{aligned} \epsilon _{\alpha \beta }^{\text{(iii)}} = {\left\{ \begin{array}{ll} \varepsilon &{} {\text{ for }} \alpha = \beta = 3\\ (1+\varepsilon )^{-1/2}-1 &{} {\text{ for }} \alpha = \beta \ne 3\\ 0 &{} {\text{ else }} \end{array}\right. } \end{aligned}$$
(1c)

for compression mode (iii), which will be referred to as density-conserving compression in the following. The scalar \(\varepsilon\) is always chosen non-positive, i.e., the system is always compressed parallel to the “3”-direction. For compression modes (ii) and (iii) we assume the “3”-direction to be parallel to each of the three unit cell vectors in subsequent runs. This allows us to lift the bias of choice for the ”unique” direction.

The two uniaxial compression modes are meant to roughly mimic the situation that occurs in a ZnP film, which is indented by a counter asperity so that the stress is largest in the direction normal to the tribological interface. Such a film will attempt to expand within the plane, but, assuming the Poisson’s ratio of ZnP films to be (slightly) positive, it will do that to a lesser attempt than if density were conserved. Thus, we see compression modes (ii) and (iii) to sandwich the real situation. In addition, in a real-laboratory friction experiment, there will be a non-diagonal stress-tensor element, whose magnitude can be non-negligible when the friction coefficient is of order unity. While it might have been interesting to consider this in-plane non-isotropy of the stress tensor explicitly, we argue that the natural fluctuations of the stress tensor, which are induced by the finite size of our system should suffice to implicitly induce such a non-diagonal stress. This happens, because the degeneracy of the two smaller stress tensor eigenvalues, which will be introduced next, is lifted as a consequence of the finite size.

2.4 Observables

2.4.1 Energy of Reaction

To address the question whether a stress-induced reaction is exothermic or endothermic, we computed the energy of reaction as

$$\varDelta E = E({\text{product}}) - E({\text{reactant}}).$$
(2)

This is done by completely relaxing the “reactants” (ZnPs before deformation) and the “products” (ZnPs after deformation) to their closest energy minima through a conjugate-gradient minimization as provided in the software package CP2K [58]. The enthalpy minimization is done at constant external (isotropic) stress without constraints on the shape of the simulation cell.

Our estimate does not include thermal effects and thus neglect the correction \(\int dT\, \varDelta c_p(T)\), where \(\varDelta c_p\) is the difference of specific heats between product and reactant and T is the temperature. Assuming that \(c_p\) of both product and reactant do not deviate substantially from the rule of Dulong Petit (at least in classical treatments of the nuclei), \(\varDelta c_p(T)\) is the difference between two small differences from that rule. This is why we believe that \(\varDelta E\) is a reasonably accurate measure for the (experimentally relevant) free energy of reaction with room temperature and ambient pressure being close to the relevant thermodynamic reference state. In other words, we believe entropic and anharmonicity effects on the free energy difference between products and reactants to be minor.

2.4.2 Stress-Tensor Invariants

Before introducing stress-tensor invariants, it is in place to clarify that reported stresses are compressive stresses, in which case the hydrostatic pressure p, please welcome the first stress-tensor invariant, is nothing but the mean of the diagonal elements of the stress tensor, i.e., \(p = \sigma _{\alpha \alpha }/3\) using summation convention over identical indices. Since the stress tensor is symmetric, it has D real eigenvalues in D spatial dimensions, which fully defines the stress state of an originally isotropic system. Thus, the stress tensor is fully characterized by those D eigenvalues and the orientation of the coordinate system in which the stress tensor is diagonal with respect to the laboratory.

Rather than stating the stress states of an originally isotropic system in terms of the stress-tensor eigenvalues \(\sigma _i^{\text{E}}\), it is frequently more meaningful to state invariants that can be constructed from the various \(\sigma _i^{\text{E}}\). One such invariant is the so-called von Mises stress, \(\sigma _{\text{vm}} \equiv \sqrt{D\,J_2}\) with \(J_2 \equiv s_{\alpha \beta } s_{\beta \alpha }\), which turns out proportional to the standard deviation of the stress-tensor eigenvalues. Here, we introduced the so-called stress deviator tensor through \(s_{\alpha \beta } \equiv \sigma _{\alpha \beta } - \delta _{\alpha \beta } \sigma _{\gamma \gamma }\). In a two-dimensional system, p and \(J_2\) are the only invariants needed to specify the stress state and \(J_2\) is nothing but the shear stress in the coordinate system for which \(\sigma _{11} = \sigma _{22}\).

In three spatial dimensions, the von Mises stress remains a measure for the stress anisotropy and thus shear stress, but an additional measure for the loading type can be made with the so-called Lode angle \(\theta _{\text{L}}\) [69]. It allows one to ascertain the position of the middle eigenvalue, let’s say, \(\sigma _2^{\text{E}}\), relative to the minimum and maximum eigenvalues, which we denote by \(\sigma _1^{\text{E}}\) and \(\sigma _3^{\text{E}}\), respectively. To define the Lode angle, we first introduce the third invariant in three spatial dimensions, \(J_3\), which is simply the determinant of the just-introduced tensor s. Finally, the Lode angle is given as

$$\cos (3\theta _{\text{L}}) = \frac{J_3}{2} \, \left( \frac{3}{J_2} \right) ^{3/2}.$$
(3)

Readers may or may not want to convince themselves that the smallest possible Lode angle of \(\theta _{\text{L}}^{\text{min}} = 0\) is taken when the middle eigenvalue is equal to the smallest eigenvalue \(\sigma _3^{\text{E}} > \sigma _2^{\text{E}} = \sigma _1^{\text{E}}\), while the maximum Lode angle, \(\theta _{\text{L}}^{\text{max}} = \pi /3\), occurs when the middle eigenvalue equals the largest eigenvalue, in which case \(\sigma _3^{\text{E}} = \sigma _2^{\text{E}} > \sigma _1^{\text{E}}\). If the middle eigenvalue \(\sigma _2^{\text{E}}\) is the mean of \(\sigma _{1,3}^{\text{E}}\), then \(\theta _{\text{L}} = \pi /6\).

2.4.3 Elastic Properties

The bulk modulus B of a material specifies how resistant that material is to compression. It can be defined through the volume derivative of pressure as

$$B(p) = -V_0\left. \frac{dp}{dV}\right| _{V=V_0},$$
(4)

where \(V_0\) is the volume of the system at pressure p under the condition that the material adjusts its shape to minimize its enthalpy or rather its Gibbs free energy. We used this definition to determine the (zero temperature) bulk modulus numerically and performed the derivative dp/dV by minimizing the enthalpy at \(p = \pm 0.1\) GPa and by measuring the volume (changes) needed to minimize the enthalpy.

We also determined individual elastic tensor elements, which are defined by

$$C_{\alpha \beta \gamma \delta } = - \frac{ d \sigma _{\alpha \beta } }{d \varepsilon _{\gamma \delta }}$$
(5)

when \(\sigma\) indicates compressive stress. Finite differences were taken by setting individual strain tensor elements to \(\pm 0.001\). In the following, we will leave the tensor notation for elastic tensors and use the Voigt notation instead, in which pairs of indices are reduced to one index, i.e, \(11\rightarrow 1\) through \(33 \rightarrow 3\) and \(23 \rightarrow 4\) through \(12 \rightarrow 6\).

For the determination of elastic properties, all (reference) configurations entering the analysis were first relaxed to zero temperature and zero stress, thus allowing the simulation shape to deform, before the strain or volume changes were imposed for the measurement of the various elastic constants. After imposing the strains, the energies were relaxed again but not the simulation cell shape except for the determination of bulk moduli.

The decompression of the networks formed was done in steps of 0.5 GPa. At each pressure, the enthalpy was minimized by relaxing the atomic position and the volume using a conjugate-gradient without preconditioner so that no rate can be associated with the decompression.

To measure anisotropic stiffness changes induced during the network stiffening, we determined mean values of \(C_{33}\) and mean values of \(C_{11}\) and \(C_{22}\). The latter could be symmetrized for uniaxial compression, as the “3” axis is the symmetry axis. In a finite cell, there are yet maximum and minimum in-plane eigenvalues of the \(C_{\alpha \beta }\) tensor, where both indices are less than 3. The orientation of the coordinate system, in which this sub-tensor is diagonal, is not necessarily oriented with the simulation cell. Thus, we also computed off-diagonal elements like \(C_{12}\). More details on elastic tensor rotations are presented in Sec. III of the Supplementary Material.

If our system sizes had been very large, the results for the elastic tensor would have obeyed the symmetry relations for isotropic solids, at least, for the initial structure and those obtained by isotropic compression. Due to the finite size, we observe non-negligible deviations from that symmetry in each individual configuration. However, in all tested cases, it appeared as if the elastic tensor was similar to that of an orthorhombic crystal, in which tensor elements \(C_{ij}\) were rather small if one of the two Voigt indices was \(\le 3\) while the other was \(\ge 4\). Also, while fully relaxed configurations always happened to be triclinic, the maximum deviation of the angle between any two unit cell vectors from 90\(^\circ\) always turned out less than 3\(^\circ\). This obviously led to corrections to the calculated (smallest and largest) in-plane values of \(C_{11}\) and \(C_{22}\) of order 1%. As we are interested in trends rather than in precise numbers, we neglected these contributions for the sake of simplicity.

2.4.4 The Mixed Radial, Angular, Three-Body Distribution Function

The most frequently studied function from which information on local structure is deduced from molecular simulations is the radial distribution function (RDF) g(r). It states the probability density to find an atom (of a given type) at a distance r in units of the mean density of that atom type. It plays an important role because it can be directly linked to the structure factor and thus to diffraction patterns. Unfortunately, bond angles are difficult to deduce from g(r), in particular in non-elemental and/or disordered systems. While bond-angle distribution functions (ADFs) are, as their name suggests, sensitive to bond angles, they cannot “see” past the nearest neighbors.

The recently proposed mixed radial, angular, three-body distribution function (RADF) [70], \(g(r_{ik},\cos \vartheta _{ijk})\), contains implicitly most information from RDFs and ADFs. However, it allows additional insight to be gained in a single graphical representation, such as, typical angles between nearest-neighbor and next-nearest-neighbor “bonds”. Specifically, \(g(r_{ik},\cos \vartheta _{ijk})\) can be defined as the probability density of finding an atom k at a distance \(r_{ik}\) from atom i when the angles between vectors \({\mathbf{r}}_{ik}\) and \({\mathbf{r}}_{ij}\) take the value \(\vartheta _{jik}\) under the condition that the atoms i and j are nearest neighbors [70]. While there is some ambiguity to the precise choice beyond what maximally allowed bond length \(d_{ij}^{\text{max}}\) two atoms i and j cease to be considered neighbors, the precise value for \(d_{ij}^{\text{max}}\) usually does not play a significant role when it is chosen with moderate care. Only Peierls or Jahn-Teller distorted systems prove difficult to treat. For more (mathematical) details on this distribution function, we refer to the original literature [70].

In this study, we measure \(g_{\text{ZnOO}}(r_{ik},\cos \vartheta _{ijk})\) in which case the atom i must be a Zn and the two remaining atoms j and k oxygens. As the Zn-O RDF shows well separated first-neighbor and second-neighbor peaks, the precise choice of \(d_{\text{ZnO}}^{\text{max}}\) is uncritical. In order to lift the remaining, small ambiguity, a mean bond length \(d_{\text{ZnO}}^{\text{mean}}\) was deduced from a skewed-normal distribution (SND) analysis, as described in Ref. [71], and the standard deviation of the bond length was added to this number to yield typical values of \(2.2\pm 0.05\) Å for \(a_{\text{ZnO}}^{\text{max}}\).

3 Results

As described in the method section, we exposed an initially cubic simulation cell containing two triphosphate (P\(_3\)O\(_{10}\)H\(_5\)) and two zinc phosphate molecules (Zn[PO\(_4\)H\(_2\)]\(_2\)) to various deformation modes, specifically (i) an isotropic compression, (ii) a simple uniaxial compression, in which the strain in the plane normal to the compression axis was set to zero, and (iii) a density-conserving uniaxial compression. In real experiments, deformations induced by an indenting tip should fall roughly between modes (ii) and (iii). Due to our systems being relatively small, stress-tensor elements can deviate from those that would be expected macroscopically from the choice of stress-tensor elements. This puts us into a position to argue that the eigenvalues of the stress tensor also occasionally took values that would be characteristic for a ZnP film below a sliding tip. The deformations were imposed until the measured energy changed substantially in a quasi-discontinuous fashion. Systems were then uncompressed. Heat of reactions as well as elastic properties reported in the following below refer to the uncompressed samples.

In the initial configurations, zinc is tetrahedrally coordinated as shown in Fig. 1a. For completeness and later discussion, we report the bulk modulus of the initial sample to be \(B \approx 30\) GPa and its Poisson’s ratio to be close to 1/3. This translates to a Young’s modulus \(E \approx B\) and an indentation modulus \(E^*=E/(1-\nu ^2) \approx 34\) GPa.

All systems obtained after a full compression/decompression cycle reveal that the coordination of one of the two zinc atoms changed to a seesaw geometry while the other adopted a square pyramidal geometry, as shown in panels (b)-(d). At the same time, the shape of the simulation cells shows a hysteresis, which is indicative of the previously imposed deformation, i.e., the simulation cell that was uncompressed from a density-conserved deformation looks the most flattened while the one arising from the isotropic deformation has the least modified shape.

Although the changes induced in the coordination of zinc atoms is similar for all three compression runs, central properties differ between them. The deformation-induced reaction is most exothermic for the isotropic compression, see Table 1 for precise numbers.

Table 1 Various properties in the vicinity of the deformation-induced structural change in a zinc phosphate model system for different deformation modes: energy of formation \(\varDelta E\), bulk modulus B of the uncompressed samples, and three stress-tensor invariants (hydrostatic pressure p, von Mises stress \(\sigma _{\text{vm}}\), and \(J_3\)) as well as the Lode angle just before and after the zinc atoms changed their coordination in the compressed state
Full size table

The general trend of low-energy structures of a given stoichiometry being stiffer than high-energy structures is also followed by the investigated ZnPs, as revealed by the bulk moduli B listed in Table 1: the lower the energy of the structure, the stiffer it is. Not only energy of reaction and stiffness of uncompressed systems differed between the compression modes but also the hydrostatic pressure, p, at which the irreversible deformation occurred, as well as the other stress-tensor invariants. We will later come back to their discussion further below.

3.1 Structural Properties

Although we computed various radial distribution functions in detail, especially those associated with Zn–O and Zn–P, we did not find their analysis particularly beneficial for the detection and characterization of the structural changes in the ZnPs. Most changes in the RDFs were subtle shifts of peak positions and intensities or the enhancement of shoulders. However, a quantitative analysis of the Zn–O RDF in terms of a skewed-normal distribution analysis of the nearest-neighbor peak [71] revealed a quasi-discontinuous change in the mean coordination number Z, which is depicted in Fig. 2 (the way how Z is obtained from the RDFs is shown in Fig. S1–S3 of the Supplementary Material). Specifically, Z increases abruptly from \(Z \approx 4\) to \(Z \approx 4.5\) during the deformation and remains close to 4.5 for larger deformations, until it increases discontinuously again at much larger compressions, which we do not discuss here.

Fig. 2
figure 2

Mean coordination number Z as function of reduced strain \(\varepsilon /\varepsilon ^*\) for the three compression modes, \(\varepsilon ^*\) being the value of the strain at which the reaction is observed. Values for \(\varepsilon ^*\) are (i) isotropic compression \(\varepsilon ^* = -0.08\), (ii) simple uniaxial compression \(\varepsilon ^* = -0.24\), and (iii) density-conserved compression \(\varepsilon ^* = -0.30\)

Full size image

To gain further insight on the structure, we computed the mixed radial, angular, three-body distribution function [70] introduced in Sect. 2.4.4. Results on RADFs are presented in Fig. 3 for the original structure in panel (a) and those that were obtained right after the hybridization change on zinc atoms had occurred, see panels (b–d). For the latter cases, the RADFs remain unchanged to the eye by a decompression to zero external stress. Fig. 3 reveals bond angles close to the ideal tetrahedral bond angle for the initially disconnected ZnPs, as reflected by a broadened peak at \(\cos \vartheta = -1/3\) at typical nearest-neighbor distances. This finding goes in line with the representative snapshot of the simulation box in Fig. 1a.

Fig. 3
figure 3

The mixed radial, orientational correlation function of triplets O–Zn–O of a the reference structure at pressure of 0.5 GPa, b after isotropically compressing the reference structure to a pressure of 5.0 GPa, c after a uniaxial compression along z-axis, and d after a density-conserving deformation. All data are taken under temperature of 600 K. The black circle in subfigure a shows the position of the peak for a perfect tetrahedral structure. The black pluses and crosses in subfigures b, c and d indicate positions of the peaks for perfect seesaw and perfect square pyramidal geometries, respectively. Distances are normalized by the mean Zn–O bond length

Full size image

For all three investigated compression modes, the bond angles take values near 90\(^\circ\) and 180\(^\circ\) after the stress-induced reactions occurred. Unfortunately, the analysis of RADFs does not allow seesaw and square pyramidal geometries to be distinguished from each other, as both have the same relative number of 90\(^\circ\) and 180\(^\circ\) bond angles. For a sample as small as ours, it is then easiest to make that distinction by visual inspection, which revealed for each compression mode one zinc atom to adopt a seesaw geometry and the other square pyramidal. Qualitative differences between the various structures can at best be detected by intensities arising from oxygens at distances r in the range \(1.5 \lesssim r/d^{\text{max}}_{\text{ZnO}} \lesssim 2\), i.e., in the second neighbor shell of zinc atoms.

3.2 (Critical) Stress Tensor (Invariants)

The critical stresses, or rather, the critical stress-tensor invariants, where the hybridization changes on zinc atoms occurred varied quite distinctly between the compression modes, e.g., at a critical hydrostatic pressure of \(p_{\text{c}} = 4\) GPa for the isotropic deformation and at \(p_{\text{c}} = 1.23\) GPa for the density-conserved compression. At the same time, the critical von Mises stress was noticeably larger for the density-conserved than for the isotropic compression, i.e., \(\sigma ^*_{\text{vm}} = 2.09\) GPa versus \(\sigma ^*_{\text{vm}} = 1.23\) GPa, see Table 1 for more details. Here we report the last available deformation before the transition happens, at which the crystallographic positions correspond to the thermal equilibrium positions. Post-reaction stress-tensor invariants for deformation-controlled simulations were obtained by letting the newly formed structure adopt the last cell shape, for which the reactant had still been stable.

While increasing off-diagonal stresses clearly reduces the hydrostatic pressure at which the transition occurs, there is no substantial reduction of the shear stress after the transition occurred in isotropic or simple uniaxial compression. We would therefore argue that the transition is driven by the hydrostatic pressure but assisted by the shear stress, as increasing \(\sigma _{\text{vm}}\) clearly decreases the transition pressure. However, for the density-conserving uniaxial compression, \(\sigma _{\text{vm}}\) drops distinctly more after the transition than the hydrostatic pressure. It might thus be in place to call the network stiffening to be shear driven and (potentially) pressure assisted in that particular case.

The origin of non-negligible (critical) von Mises stresses arising in response to an isotropic deformation can be rationalized as a system-size effect. It would disappear if we started the simulations from stochastically independent initial configurations and determined the expectation values of stress-tensor elements before deducing the von Mises stress. We could also make it disappear by symmetrizing the stress tensor itself with allowed symmetry operations, i.e., by relabelling the names of x, y, and z axis and/or by changing the handedness of the used coordinate system. However, an average of the von Mises stresses of individual configurations would always lead to a finite value, which however, disappears with the inverse square root of the system size, according to the law of large numbers. Yet, microscopic stress fluctuations arise at small scales also for macroscopic systems and we found it useful to investigate how this local stress affects chemical changes in ZnPs.

3.3 Elastic Tensor Anisotropy

Since the stress anisotropy breaks the (expected/average) symmetry of the system, the elastic properties of the “glassy” ZnPs obtained after a full compression/decompression (c/d) cycle may turn out anisotropic, even if they remain disordered. However, for the initial system and the one obtained after a full isotropic c/d cycle, violations of elastic isotropy conditions, such as \(C_{11} = C_{22}\), \(C_{12} = C_{23}\), \(C_{16} = 0\), or \(C_{66} = (C_{11}+C_{22})/2\), arise as finite-size effects. However, they tend to be relatively small, i.e., typically \(<0.05~B\) for our initial system and \(<0.1~B\) for the isotropic c/d cycle.

The elastic anisotropy is distinctly enhanced after a density-conserving compression and even more so after the simple uniaxial compression, as depicted schematically in Fig. 4. The “soft direction” turns out to be the one in which the ZnPs had been most compressed. We rationalize this observation as follows: Atoms are squeezed deeply into repulsion during the compression, and the structural relaxation attempts to reduce the amount of most extreme repulsion in the compressed state, i.e., repulsive forces acting in the ’3’-direction. After the transition and after the removal of the external stress, atoms relax most (out of the repulsion) in the direction of the originally highest stress. Recompression along that direction is consequently easy to achieve, which explains the relatively small values of \(C_{33}\).

Fig. 4
figure 4

Elastic tensor elements of different ZnP structures of a the original structure and those of those obtained after decompressing from a b isotropic, c simple uniaxial, and d uniaxial, density-conserving compression

Full size image

We speculate that if the critical stress tensor had a Lode angle approaching 60\(^\circ\), we would have obtained two stiff and one soft direction. However, this loading condition would be atypical in a tribological situation, which is why we did not consider it in this study.

4 Conclusions

In this work we studied how a system built of ZDDP model decomposition products—two triphosphate molecules and two zinc phosphate molecules—reacts to different deformations, which included one isotropic and two uniaxial compressions, one of which conserved density, while the other kept all strain-tensor elements constant except for one diagonal tensor element. In all deformation modes, we observed that one zinc changed its initial tetrahedral coordination to a seesaw geometry while the other converted to a square pyramidal structure. Although both system size and relative number of zinc atoms were small, we believe our observations to be characteristic for ZnPs: in crystalline \(\alpha\)-\(\hbox {Zn}_3(\hbox {PO}_4)_2\), two thirds of Zn atoms, which are all coordinated tetrahedrally initially, change their local environment from tetrahedral to seesaw, while the remaining ones adopt a local coordination of five as in a square pyramidal geometry. Also the hydrostatic pressures, where the changes occur in the absence of significant shear stresses, differ by a factor of two, i.e., \(p_{\text{{c}}} = 5\) GPa for amorphous ZnPs and 9 GPa for the crystal. The latter is reported in Sect. II of the Supplementary Materials.

Our simulations corroborate the conjecture originally proposed by Mosey et. al. [30] that mechanical stress is the decisive factor to promote a hybridization change on zinc atoms, which is needed to activate the anti-wear functionality of ZnP films. However, in addition to the previous correction of the overestimation of the needed hydrostatic pressure \(p_{\text{c}}\) to induce irreversible coordination changes on zinc atoms from originally 18 GPa [30] to 5–7 GPa [31,32,33], we support the refinement of the theory in which shear stress is argued to reduce \(p_{\text{c}}\). In addition, we propose that the elastic properties of the ZnPs depend sensitively not only on \(p_{\text{c}}\) but also on the values of other stress-tensor invariants, most notably the von Mises stress at the point where the films undergo structural changes. Films generated predominantly by shear stress turn out comparatively soft, potentially too soft, like thermal ZnPs films, to protect surfaces from wearing [72]. This speculation is certainly consistent with the lower durability of ZDDP tribofilms formed in the absence of asperity contact [27].

Films generated under one of the two uniaxial compressions turn out noticeably softer in that direction than in the two remaining directions. Peak stresses in tribological contacts tend to be largest in the direction normal to the interface, at least as long as the friction coefficient remains less than unity. Thus, we expect films to be softest in the direction normal to the interface. We expect this elastic anisotropy to allow the films to be sacrificial under large shear stresses, even if it is significantly less than that of true layered compounds like graphite or molybdenum disulfide.

An interesting test of our stiffening hypothesis would be to repeatedly indent a thermal film or a film that had been produced through non-contact shear. The first indentation should reveal a relatively soft film, which we would expect to remain soft, until the applied load becomes sufficiently large to induce the coordination changes discussed in this work. On retraction and in subsequent indentation, we would expect the measured modulus to have increased and the film to be more wear resistant. Ideally, such an experiment would be done in conjunction with a chemical characterization so that the stiffening can be ruled out to have occurred in response to stress-induced stoicheometric changes, as we believe it to have been the case in films having been exposed to humidity for several years [34].

In fact, a process that our simulations do not capture but which we believe to be very important in real systems is that our periodically repeated simulations cell does not allow individual atoms to disappear (automatically). However, in the vicinity of the crosslinking stress state, a significant rearrangement of hydrogen atoms occurs. We expect some of them to break lose in reality and to drift in a direction opposite to the pressure gradients, whereby dangling bonds in the remaining ZnP film would need to be saturated by other dangling bonds rather than by hydrogens, or by other (small) radicals that are produced through large local stresses. This would obviously enhance the connectivity and thus the stiffening of the network. Support for this idea comes from the observation that some hydrogen atoms exhibit large bond-length fluctuations after decompression, which is indicative of a reduced chemical stability of these bonds.