Mechanochemical Ionization – Differentiating Pressure-, Shear-, and Temperature-Induced Reactions in a Model Phosphate

Using density-functional theory based molecular dynamics simulations, we study stress and temperature-induced chemical reactions in bulk systems containing triphosphoric acid and zinc-phosphate molecules. The nature of the reactants depends sensitively on the imposed conditions, e.g., isotropic and even more so shear stress create ionic reactants. Instead, zwitterions and water molecules emerge from thermal cycles. Hydrostatic stresses required for reactions to occur lie well below those typical for tribological micro-contacts and are further reduced by shear. During the moment of impact, many protons are mobile, even if they bond to groups that formally remain inert in a compression/decompression cycle. Our ﬁndings are consistent with a variety of observations on the triboelectric charging of insulators but contradictory to the hot-spot model.


Introduction
Mechanochemistry is primarily concerned with the interplay of mechanical forces and chemical bonds [1][2][3] . This includes the analysis of the breaking or unfolding of individual molecules under mostly tensile forces, which is of great scientific interest 1,4 , but also, stress-induced chemical changes in bulk systems and/or at interfaces 2,5 . These latter processes, which are central to tribochemistry and its independent subfield ball milling 6 , are arguably of greater applied relevance. Two questions in mechanochemistry are surprisingly unexplored. First, how does the precise form of the stress tensor affect stress-induced chemical changes? Most works, including influential reviews 1, 2 , treat stress as being small or large. However, mathematically speaking, stress is a three-by-three matrix and not simply a scalar. Thus, for certain reactions to be preferentially induced in an originally isotropic system, the three eigenvalues of the stress tensor may have to approach three independent target numbers closely. In other words, generating an optimum condition for a certain reaction to occur could necessitate the simultaneous tuning of shear and hydrostatic stress plus one additional condition, and even more conditions for ordered systems. Second, how do stress and its anisotropy affect the generation of ions and/or their mobility? The answer to this question could be central to correctly explain various observations made on the friction-induced charge transfer between two rubbing, large band-gap insulators, which is one of the earliest studied natural phenomena 7 . In fact, as McCarty and Whitesides 8 discuss very thoroughly, tribo-charging is difficult to explain when not only electrons but also ions appear to be immobile. Thus, a central question to be addressed is how (contact) stress affects number and mobility of ions and not only the electronic band structure.
There is an arsenal of bulk systems for which the precise form of the stress tensor could affect the chemical structure after the external stress is released, irrespective of whether or not it involves the separation of an interface. We decided to first focus on phosphate-based systems, because their properties appear to be particularly sensitive to external stresses: the formation and functionality of zinc-phosphate (ZnPs) based anti-wear films on rubbing surfaces 9, 10 as well as their patchiness with stiff films occurring on the highly loaded asperities and soft films in the valleys 11 was interpreted as a stress hysteresis 12 . Experiments supported the conjecture that the fast formation of stiff films benefits from hard substrates in that mechanical properties of the substrates were found to be more critical than their stoichiometry 13,14 . Metal phosphates are also used for many other applications, e.g., as electrolytes in Galvanic 15 and fuel cells 16,17 . The question arises to what extent the conductivity of (metal) phosphates [16][17][18][19][20][21] can be increased through the application of mechanical stress, given that computer simulations of triphosphoric acid molecules indicated a proton transfer reaction to occur between originally identical and neutral molecules at a hydrostatic pressure of p ≈ 3.5 GPa, whereby free cations and anions were produced 22 . Related stressinduced processes may explain why the conductivity of a 1:1 mixture of CsHSO 4 and CsH 2 PO 4 was substantially enhanced through ball milling 20 . The interplay of mechanical stress and phosphate chemistry even appears in organic systems, where mechanically stressing integrin receptors was observed to enhance the tyrosine phosphorylation of cytoskeletally anchored proteins 23 . These few examples are certainly only the tip of the iceberg of cases, in which the stress and its anisotropy affect the chemistry of phosphates and the tip of an entire mountain range if other group-5-element containing (molecular) solids were included.
We have recently studied the response of a bulk system containing triphosphoric acid and zinc phosphate molecules to externally imposed deformations 24 . All compression/decompression cycles lead to exothermic chemical changes, however, the reaction energy, i.e., the energy difference between the fully relaxed initial and final structures, generally turned out largest in magnitude under isotropic deformation and smallest when the deformation was most anisotropic. The bulk modulus of the decompressed structures revealed similar trends. Unfortunately, our previous analysis of the structural or chemical changes in terms of two-and three-body correlation functions did not allow us to rationalize our results. The original motivation of this work was to achieve this with a more detailed analysis of the chemical bonding. We discovered, much to our surprise, that the Lewis structures of the configurations, which were produced while compressing three-dimensional, orthorhombic simulation cells, could be represented in terms of two-dimensional Lewis structures. This feature allowed us to interpret the stress-induced changes in zinc phosphates in rather simple chemical terms and to reveal a process, which might be dubbed stress-induced (zwitter-) ionization. Such a process is likely to have immediate consequences for rubbing or contact induced static charge, also known as triboelectricity. While we analyze to a large extent the structures produced during our previous study, we add many additional simulations in which the initial dense structure is heated and cooled again or in which the stability of decomposition products is analyzed. In addition, we run compression/decompression simulations for nylon, in order to address the question why it tends to pick up positive tribo charges.

Results and Discussion
Starting point of our simulations are two Zn[PO 4 H 2 ] 2 and two P 3 O 10 H 5 molecules. The response of these molecules to isotropic stress was investigated before 12 to rationalize the stress-induced formation of anti-wear films. This is why we (rightfully) expected interesting mechanochemistry to happen. The molecular Lewis structures are shown in Fig. 1(a) and (b), respectively. The four molecules are placed in a cubic simulation cell with a volume of V = (3 nm) 3 , and then compressed to V = (1.2 nm) 3 within 50 ps at temperature of T = 600 K. Next, an external pressure p = 0.5 GPa was applied and the system equilibrated for 20 ps at the same temperature. During this time, all covalent, ionic, and mixed bonds remain intact. However, due to the rotation of terminal groups and other molecular rearrangements, the hydrogen-bond network changes over a-few-picoseconds timescales. Although the simulation cell is kept cubic at all times and the final, dense configuration is space-filling, see Fig. 1(d), it is possible to represent the structure as a single molecule, formally lying in a two-dimensional plane, i.e., all pairs of atoms connected through either covalent or strong ionic bonds can be drawn as neighbors, which is done in Fig. 1(c). Dashed arrows indicate ionic bonds. Zinc, phosphorus and oxygen are presented in gray, orange and red colors, respectively. Formal positive charges are always assigned to the phosphor atom in a phosphate group rather than to hydrogen atoms. Two colors were used for hydrogen atoms when drawing Lewis structures: red for those that are bonded to the same oxygen atom before and after the compression and green otherwise.

3/8
After thermal equilibration, the system was compressed in three different modes: (a) isotropic compression, (b) areaconserving, uniaxial compression, and (c) volume-conserving uniaxial compression. During (b), the area normal to the direction of compression remained fixed, while in (c) it was scaled such that the volume of the simulation cell was constant. Typical deformation modes in a high-pressure diamond-anvil-cell experiments are close to mode (a). Ball milling and tribological contacts would frequently lie somewhere between (b) and (c) whenever normal contact forces clearly exceed lateral forces, a situation which we refer to as "implicit or no explicit shear". The precise location of such a stress state would depend on the Poisson's ratio ν of the compressed material, e.g., it would be close to compression mode (b) for ν = 0 and to mode (c) for ν = 1/2. Compressions were terminated when the enthalpy suddenly decreased quasi-discontinuously, which is indicative of a major structural rearrangement and/or a chemical reaction. Then, all degrees of freedom were allowed to fully relax to zero stress and minimum energy, that is, every atomic coordinate as well as the six degrees of freedom specifying the simulation cell, which are three lengths and three angles.
The Lewis structures of the decompressed system are shown in the left column of Fig. 2 (a)-(c) reflecting, in this order, deformation modes (a)-(c), while the three-dimensional structures are shown in the right column. Some general observations, valid for each of the investigated compression/decompression cycle, can be made. One zinc atom remains four coordinated with oxygen while the other forms an additional bond with another oxygen atom. Distorted tetrahedral and distorted square pyramidal molecular geometry are adopted around four and five-coordinated zinc atoms, respectively. An analysis of the electron density between nearby zinc and oxygen atoms reveals that each zinc has mixed bonds with two oxygen atoms (with most bond lengths in between 1.85 and 2.05 Å) while remaining bonds are purely ionic (predominantly between 1.95 and 2.15 Å). In each structure, at least one net proton transfers from a phosphate group to another one occurred. A significant fraction of hydrogen atoms changed the oxygen atom that they had been bonded to, even if the phosphate moiety that the mobile hydrogen atoms were bonded to did not participate in the stress-induced reaction. This means that during the brief moment of mechanical impact, i.e., of being in a stressed state, many protons were essentially free to move around. Immediately after the final coordination on the zinc cation was achieved, only a few protons remained mobile in the compressed state, namely those associated with charged groups. Even those hydrogen atoms became immobile after the systems were relaxed to the next available, zero-stress energy minimum, and then re-equilibrated at 450 K. During the initial compression, hydrogen atoms started to become mobile, even those not formally participating in the reaction, when the deformation exceeded roughly 50% of the critical deformation, at which enthalpy suddenly decreased. Given that hydrostatic pressures p and shear stresses τ (to be precise, the second deviatoric stress invariants), at which the transitions are expected to occur quasi-instantaneously are approximately 24  Systematic differences between different deformation modes could also be observed. The degree of structural rearrangement increases, not surprisingly, with increasing deformation anisotropy, i.e., with increasing (implicit) shear. The more anisotropic the deformation, the more linear becomes the final Lewis structure. The spatial separation of cations and anions increases with the degree of anisotropy, which is consistent with the lower energy of reaction for anisotropic deformation reported previously 24 . For mode (a) and (b) only one net or charging proton transfer occurs, while mode (c) leads to two proton transfers. These statistics are produced irrespective of the choice of the compression axis, although the precise nature of the reactants differed each time and it appeared to be random which of the two zinc atoms became penta-coordinate. A similar stress-induced coordination change in Zinc-orthophosphate [α-Zn 3 (PO 4 ) 2 ] was observed to occur near a hydrostatic pressure of 9 GPa in a combined experimental, theoretical study and attributed to the ability of Zn to grasp on to additional oxygen atoms by hybridizing a d orbital with a p orbital of a bridging oxygen into a σ bond 25 . The insensitivity of our main results (number of produced ions, energy of reaction, final bulk modulus, number of mobile atoms during high compression) on the direction of compression in a sample as small as ours could result from the relatively high symmetry of the tetrahedral bonding on zinc and phosphorous atoms and the relatively large number of terminating OH groups making them originally point in quasi-random directions.
A few additional details are worth discussing. All deformation-induced reactions generally lead to the formation of free ions, although zwitterions are also produced occasionally. In case of the shown reactant for mode (a), the cation is formally a protonated phosphoric acid molecule, which we occasionally also observed in modes (b) and (c), but which would probably not form in gas-phase or solution chemistry. We note that we find an isolated, neutral P(OH) 4 molecule to spontaneously decompose into a hydrogen radical and a phosphoric acid molecule, while an isolated, singly charged P(OH) 4 cation remains stable.
To compare stress with temperature-induced reactions, we heated the system to 5,000 K, which are often believed to occur at small scales in tribological contacts, while avoiding evaporation by setting the pressure to p = 0.5 GPa. Temperature was then brought back down in discrete steps. At T = 4, 000 K, 2,000 K, 1,000 K, 600 K, 400 K, and 250 K, the system was allowed to equilibrate for 5 ps, which was long enough for the enthalpy to reach a regime where signs of relaxation 4/8 were hidden in thermal noise. This protocol made the cooling process last approximately as long as a decompression. Yet, the thermal cycle was endothermic with an energy of reaction of ∆E = 0.507 eV per Zn atom, while all stress cycles were exothermic with ∆E ranging from −0.077 eV to −0.195 eV. More importantly, the reactants, which are shown in Fig. 3, differ substantially between thermal and stress cycles, for example, in that water molecules were produced through heating but not through stressing the system.
In addition, the thermal cycle promotes only zwitterioniziation while the stress treatments produce predominantly free ions. Another substantial difference between our stress and thermal cycles is that no oxygen atom became unbound from a phosphor through stress at 600 K, but all P-O bonds broke through temperature, since the system turned into an ionic liquid at 5,000 K.
Given our results, the occasionally made assumption that mechanochemistry results predominantly from local heating, as advocated in the hot-spot model 26 , is certainly false for our system of interest, and, as we believe, for most mechanochemical systems that we are aware of. In fact, we expect solids that are mechanically produced to be softest in the direction of largest compressive stress after the stress is released 24 , while temperature-induced reactions of originally isotropic system should not have preferred stiff directions after cooling, unless anisotropic single crystals happened to be produced. In addition, we are somewhat skeptical of typical flash-temperature estimates as they tend to assign the dissipated energy to a small interfacial zone, although important dissipative processes contributing substantially to friction, such as viscoelastic or plastic deformation, also take place far away from it. We believe our results on stress-induced ionization to matter in the context of triboelectcricity. Although the required local stresses of order 1 GPa may appear large from a continuum perspective, linear-elasticity theory predicts them to be of order root-mean-square height gradient times the smaller Young's modulus of the contacting materials 27 . They can exceed the macroscopic hardness, because plasticity-inducing defects are not energetically favorable at small scales and/or because of inertial confinement in a short-lasting asperity collision. Thus, even contacts formed by amber and keratin (wool, fur, epidermis of skin, etc.), which both have relatively large Young's modulus and form the basis of the classical amber against wool or cat-fur triboelectricity couple 7 , must be expected to reach local values exceeding 1 GPa. Nernst 28 speculated a long time ago that tribocharging between insulators is due to ion transfer. He proposed various equations describing how ions account for contact-induced charging in response to gradients or interfacial differences in chemical composition, ion concentration, temperature, pressure, and the like. Central parameters are ionic mobilities and number densities, which are utterly sensitive to stress in our model phosphate at typical contact stresses. It may appear daring to discuss triboelectricity w.r.t. our phosphates, in light of its (current) absence in popular triboelectric series 29,30 , in which material A is placed above material B if A acquires a positive charge through rubbing against B. Nonetheless, dry phosphate ores seem prone to triboelectric charging, because they can be beneficiated with triboelectric belt separators 31,32 .
A careful review of the literature 33 should prevent one from making strong claims of one mechanism, let alone one material property, to be the sole key to understand tribo-induced charging between insulators, even if it is enticing to reduce triboelectrical series to those materials allowing them to be ordered according to the dielectric constant 34 , the pK a 35 , the Seeman coefficient 36 , or the expected flexoelectric potential differences between two materials in contact 37 . As is the case for the origin of friction, where strikingly different microscopic mechanisms can lead to identical macroscopic friction laws 38 , triboelectricity being one of them 39 , the question to be answered is not what single mechanism explains everything but which mechanism dominates under what circumstances and how can we figure out, theoretically or experimentally, which one it is? In fact, the existence of a cyclic triboelectric series 8 , where material A is more tribopositive than B than C than A-somewhat reminiscent of the stairs in the famous lithography Ascending and Descending by the Dutch artist Maurits Cornelis Escher (1960)-appears to indicate that there is not a single dominant tribo-charging mechanism for the AB, BC, and CA interfaces. Likewise, the reversal of charging direction with increasing rubbing time 34,40 can be seen as an indicator for the existence of 5/8 competing tribo-charge carriers potentially having different response times.
While regular ion-transfer charging can already be at the root of various triboelectric observations, such as the charging between dielectrics of identical chemistry but different size 41,42 or mosaics of positive and negative domains 43 , additional phenomena might have their root in stress-enhanced ionization. This concerns in particular how strain alters the position of some materials in the triboelectric series 44,45 and the correlation of the latter with their pK a 35 , which appears natural, since a material with a large propensity to donate protons to an electrolyte will also release them under stress. Such an effect could even radically change the positioning of a material in a triboelectric series, e.g., when it has stress-releasable protons in addition to nucleobase moieties. In well-defined, single-asperity contacts, free ions would only be expected to occur above a certain threshold force and explicit stress would enhance but not be required for tribocharges to be generated. In randomly rough surfaces, the repositioning of a material in a triboelectric series due to increased surface roughness 46 could simply result from an increased rms-height gradients and thus increased contact stresses. We certainly do not mean to imply that the observation of every discussed phenomenon automatically implies stress-induced ionization to be responsible for it. However, given the ease with which ions can be apparently created under stress, our bias has certainly shifted toward the "ion-transfer school", despite us having made large efforts in the past toward the design of empirical electron-transfer potentials describing tribocharging between dielectrics 47 .
Of course materials acquiring their tribocharge through ions transfer while taking a prominent role in the triboelectric series do not necessarily release ions under stress. To substantiate this claim, we studied the response of nylon, which always appears at the very positive end of the triboelectric series. When compressing nylon 6,6 and nylon 7,7 first in the direction parallel to the molecular backbones and next in the direction orthogonal to the molecular backbone, stresses beyond 10 GPa were reached, which distinctly exceeds the Young's modulus of nylon. However, no single bond was broken even when heating the systems to 400 K. A natural explanation for nylon's prominent role in the triboelectric series is the presence of terminal amine groups.
To summarize, our simulations on phosphates revealed a significant stress-induced (zwitter) ionization at stresses that appear unavoidable in tribological contacts, where similar molecules are used as anti-wear and anti-oxidant additives. The final molecular structures and the number of stress-induced ions turned out to depend sensitively on the precise shape of the stress tensor during compression and differed much from products obtained through thermal activation at moderate hydrostatic pressure. This observation provides theoretical support but also refinement of the experimentally acquired picture that tribocharging can originate from the creation of radicals through mechanochemistry 48,49 . The large number of ions produced in our simulations-protons were mobilized long before the quasi-discontinuous changes in enthalpy occured, e.g., at p ≈ 0.6 GPa and τ = 2 GPa-should certainly suffice to create a net imbalance of one elementary charge per 100,000 surface atoms, which is the value needed to yield tribovoltages of several thousand volts on insulators 45 . The targeted exploration of the stressinduced ionization via simulations might benefit the search for triboelectrically active materials. Stress-released ions could lead to large tribo-voltages on a first stroke and thus be desired. However, they could also be detrimental should it be difficult to chemically rejuvinate the material after exploitation of the tribovoltage. In any event, it seems clear that gas-phase chemical considerations for the functionality of materials in contact situations might have to be augmented with analysis similar to ours, in order to ascertain if or how a certain substance or material fulfills its requirements for the expected contact conditions.

Model and Methods
Density-functional theory 50,51 simulations are done using the same CP2K open-source package 52 and similar methodology as our precedent paper 24 , including the Perdew-Burke-Ernzerhof exchange-correlation functional 53 with Grimme empirical corrections to the dispersion interactions 54 and Goedecker-Teter-Hutter pseudopotentials 55, 56 with the double-ζ Gaussian basis sets 57 . Energy cutoffs were 400 Ry in molecular dynamics simulations and 600 Ry in static calculations. Temperature (default value 600 K) was controlled through "canonical-sampling-through-velocity-rescaling" thermostat 58 . Constant stress and constant pressure were imposed through Nose-Hoover chain barostats 59 .
In N pT -simulations, pressure was changed in steps of 1 GPa. At each pressure, the system was given 10 ps to thermally equilibrate under isotropic volume changes 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. In minimum energy calculations, box shape and atomic coordinates were allowed to relax.

Data availability statement
Data supporting the findings of this study are available through the link, as well as from the corresponding author upon a reasonable request.

6/8 5 Competing interests
The authors declare no competing interests.