Temperature Activates Contact Aging in Silica Nanocontacts

Matthias Vorholzer, J. G. Vilhena, Ruben Perez, Enrico Gnecco, Dirk Dietzel, and André Schirmeisen Institute of Applied Physics, Justus-Liebig-Universität Giessen, 35392 Giessen, Germany Department of Physics, University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma de Madrid, E-28049 Madrid, Spain Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain Otto Schott Institute of Materials Research, Friedrich Schiller University Jena, 07742 Jena, Germany Center for Materials Research, Justus-Liebig-Universität Giessen, 35392 Giessen, Germany


I. INTRODUCTION
Around 350 A.D. the Greek philosopher Themistius reported that [1] "One may more easily further the motion of a moving slider than to move a body at rest." In spite of the ubiquity of this observation, understanding the origins of static friction is a challenging task and therefore the subject of vibrant research [2][3][4][5][6][7][8][9]. The complexity of the process can be best understood starting from the idea that virtually any macroscopic surface exhibits a significant roughness, leading to a certain number of true contact points at the interface [10,11]. The dynamics of the junctions is then governed by a broad range of complex physical-chemical processes including bond formation [12][13][14][15][16], plastic flow [17], capillarity [18,19], and atomic attrition [20]. Through the time-dependent nature of these processes the contact strength is expected to increase logarithmically with time. Indeed, it has been shown that the static friction increases logarithmically with the hold time [3,21], which is usually described in the rate-and-state models. However, an atomistic understanding of this phenomenological "law" is still missing.
In order to understand the atomic origins of aging, recent work has focused on silica contacts, which are relevant in a huge variety of technical and environmental systems [22][23][24]. It was shown that aging can be categorized into two different phenomena: quantitative and qualitative aging. The former, an established and widely accepted mechanism, describes the time-dependent increase of the real contact area due to plastic deformation and material creep [21,25]. More recently, researchers have focused on the so-called qualitative aging, i.e., the change of the contact strength per unit area usually controlled by the interfacial chemistry. This has spurred experiments based on friction force microscopy (FFM), where a single asperity contact geometry ensures a fixed contact area and any sign of aging should be dominated by qualitative changes of the interface. In fact, strong aging effects have been discovered with FFM, in particular for hydrophilic silica surfaces under ambient conditions [3]. Ab initio simulations have confirmed that covalent bond formation, i.e., the formation of siloxane bridges, can lead to the observed logarithmic increase of contact strength with time [12].
If such a multibond formation process is relevant for contact aging, it should play an important role in any contact where chemical bond formation is feasible. In particular, this should hold for atomically clean surfaces, where bond formation is not inhibited by passivating layers. Contact aging effects were indeed observed on a clean metal surface in ultrahigh vacuum (UHV) at cryogenic temperatures [8]; however, in metals bond formation and material creep are so intimately connected that is almost impossible to disentangle those two mechanisms. Logarithmic contact aging has likewise been observed for metal nanoparticles sliding on graphite [26], but in that case atomic bond formation between Sb and graphite is impossible and aging is due to atomic creep mechanisms. FFM experiments on NaCl surfaces also resulted in frictional strengthening, but here also due to creep caused by atomic attrition mechanisms [20].
Another aspect of aging due to multiatom bond formation is that those processes are inherently thermally activated. Temperature should thus strongly influence the aging dynamics and, indeed, theory predicts that the aging rate scales linearly with temperature [12], which, however, lacks experimental confirmation. Also, according to Eyring's law, the energy barriers responsible for bond formation can be significantly reduced by shear stress, leading to a load dependence of aging. Recent experiments for silica contacts in ambient conditions could not confirm this, but instead found "that the average pressure, and thus the average bond formation rate, was load independent within the accessible load range" [14]. This leaves open the question of if and how strongly chemical bond formation influences contact aging.
The most direct route to validate chemical-bond-based contact aging is thus to measure the temperature dependence of aging rates of atomically clean single asperity model interfaces. Indeed the friction of clean silica-silica contacts in UHV have been successfully described by a multibond formation process and might provide an excellent sample system [15]. However, recent slide-hold-slide experiments of such a system completely failed to observe contact aging at room temperature [27]. One reason may be that atomic bond formation for silica bonds proceeds too quickly in this case, as suggested by previous friction experiments [15]. Only at low, cryogenic temperatures bond formation rates could slow down such that they fall into experimentally accessible timescales.
At room temperature, current research indicates that qualitative contact aging for silica scales with humidity and is irrelevant in dry environments. For rock friction this would mean that in the absence of water, which is expected to be the case in seismic faults, the shear stress per real unit area remains constant, with corresponding consequence for earthquake simulations in the frame of rate-and-state models. Equally, under these assumptions, for wafer bonding processes humid and ambient environments must be considered crucial to facilitate the gradual bond formation across the interface, while using dry or vacuum conditions might be regarded as a viable route to prevent aging and wear in silicon-based microelectromechanical systems (MEMS).
Here we directly challenge the broadly accepted assumption on the absence of aging in dry SiO 2 nanocontacts through friction force microscopy experiments of dry silica nanocontacts in ultrahigh vacuum conditions. We demonstrate that the conventional FFM-based static friction measurement protocol is inadequate to observe contact aging, further rationalizing previous failure to observe this effect. Instead, we show that the time evolution of the contact stiffness gives clear evidence of qualitative contact aging in agreement with the thermally activated atom bond formation picture, including the predicted scaling of the aging rate. Moreover, we find an Eyring-type shift of the energy barriers due to shear stress. Complementary allatom MD simulations of silica nanocontacts corroborate the experimental results.

II. RESULTS AND DISCUSSION
In our experiments we used temperature-dependent friction force microscopy [28] to measure the lateral forces at the interface between an AFM cantilever covered by SiO 2 and p-type (111) Si wafers with either a native oxide layer or a 300-nm oxide layer. Since humidity is known to have a strong influence on microscopic and macroscopic contact aging [12,13,29], all experiments were conducted under ultrahigh vacuum conditions.
The measurement protocol was inspired by the slidehold-slide procedure introduced by Li et al. in their seminal work on this topic [3]. Essentially, the tip is elastically pulled laterally along the substrate, then halted for a welldefined time t hold , after which the pulling process is restarted [ Fig. 1(a)]. The sliding friction F sliding can be derived from the average value of the force profile in the first sliding phase [corresponding to the green region in Fig. 1(b)]. After a sudden drop, which is related to cantilever creep during the hold time, the lateral force increases almost linearly when the cantilever is pulled again along the surface. In this hold phase the tip is effectively pinned to the sample, and once sliding is reinitiated, a force constant can be calculated from the slope of the force profile which corresponds to the effective spring constant k of the pinned contact region combined with the AFM cantilever under torsion. The cantilever creep is related to the close match between cantilever torsion and friction force but is effectively limited to the first few milliseconds of the hold time (see the Appendix B) and is thus not significant for the general interpretation of the slide-holdslide process where the contact can still be assumed as static. The instant at which the contact is broken corresponds to a maximum value defining the static friction F static . At this point, the tip motion is resumed and the lateral force starts its decay toward the initial sliding friction value, which occurs within a certain transient time.
Altogether, the force variations were recorded along a path of 300 nm and the measurements were repeated between 160 and 350 K.
As a first result, we observe that both the static friction and the sliding friction decrease with temperature T [see Figs. 1(c) and 1(d)] in line with the well-known picture of thermally activated slip originally proposed by Prandtl for atomic-scale sliding friction [30][31][32]. In this model thermal fluctuations lead to a premature slip of the contact, which would otherwise occur when the theoretical shear strength of the interface at zero temperature is reached. The higher the temperature, the larger the amplitude of these fluctuations, which consequently results in a decrease of the friction force, which was demonstrated experimentally [32]. Our analysis shows that the static friction increases with contact time for all temperatures, in line with the contact aging models, but the increase is more pronounced than a simple logarithmic behavior. Even more puzzling is that the growth rate of the static friction decreases with increasing temperatures. This is in direct contradiction with the principles of thermally activated bond formation at the interface, which should be activated by temperature.
The reason for this apparent paradox is related to the dynamic nature of the static friction peak. So far, this peak in slide-hold-slide experiments has been assumed to be proportional to the number of interfacial bonds [3,12,14,27]. This is problematic, since any static friction measurement is a complex process including contact rupture dynamics on top of the aging process itself, both of which are characterized by distinct temperature dependencies. The result is a highly complex dynamic process where static friction is not a direct measure of the number of interfacial bonds. This is particularly relevant for nanoscale contacts, where bond rupturing can be significantly thermally activated [15,26,33], and it becomes especially apparent if experiments at different temperatures are compared, as is our case. This also explains why recent attempts to measure static friction of dry silica contacts have proven to be not so conclusive [27]. At room temperature, bond rupturing is strongly thermally activated [15,20,26], basically burying the dynamics of the bond formation process and making it difficult to observe qualitative aging.
A more direct measure of the number of atomic bonds is the lateral contact stiffness k. Figure 2 shows the dependence of k on the hold time t hold , extracted from the same experiments. The minor slope decrease preceding the slip is well known from the theory of atomic-scale and is caused by the slow onset of the tip slip [34]. For this reason, k is better estimated from the increase of the lateral force built up immediately after the hold phase, when the tip starts to be pulled again. If we analyze only forces well below the rupture force, the interface remains in a quasistatic state and the initial contact stiffness is a good measure for the interface condition after t hold . In an atomic bond picture, the aging of silica contacts is described by the time evolution of the number of atomic bonds at the resting interface [12,13]. Assuming that each interfacial bond corresponds to a nanospring with average spring constant k 1 along the pull direction, a parallel arrangement of bonds will simply result in an effective stiffness Nk 1 ; i.e., it will be proportional to N. A logarithmic increase of the number of interfacial bonds, Nðt; TÞ ¼ N 0 ðTÞ þ αðTÞ lnðt=t 0 Þ, should then result in an increase of the effective contact stiffness, kðt; TÞ ¼ Nðt; TÞk 1 . Indeed, a very regular logarithmic behavior of k with t hold is observed on the native oxide layer [ Fig. 2(a)], as well as on the 300-nm thin silica film [ Fig. 2(b)]. In both cases we find that higher temperature results in a faster logarithmic increase of k with the contact time [Figs. 2(c) and 2(d)]. In particular, the aging rate αk 1 scales linearly with temperature. Interestingly, the observed logarithmic behavior is also expected [35] for any two-level process, i.e., for bond formation or rupture. More specifically, a simplified analytical model proposed by Liu and Szlufarska [12], based on thermally activated bond formation, predicts that α is directly proportional to temperature with a constant prefactor, i.e., αðTÞ ¼ α 0 T, in agreement with our findings. At last, since both native and thick oxide surfaces behave similarly, this provides further evidence that only the upper atomic layers are involved in the aging process, as expected for the thermally activated bond formation aging mechanism and in line with previous simulations [12,13,36]. Note that the absolute values of contact stiffness in Figs. 2(a) and 2(b) cannot be compared directly, as two different tips were used in the two sets of experiments. Despite Figs. 1 and 2, to our knowledge constituting the first experimental evidence for thermally activated aging in dry silica contacts, there is a wealth of experiments [37][38][39] and simulation [12,13,36,37,[40][41][42] works focused on stress-induced bond formation in a-SiO 2 , indirectly supporting our findings. Namely, several MD simulations [12,13,36] predict that the indentation process results in siloxane bond (i.e., Si─O─Si) formation between tip and surface, which was also confirmed by our MD simulations; see Fig. 3(a) cyan inset. Interestingly, aberration-corrected transmission electron microscopy [37], pair distribution function analysis of 2D x ray [38], and NMR spectroscopy [43] revealed the soft nature of this bond. In fact, these works conclusively [38,[40][41][42][43] demonstrate the importance of bond angle change (of Si─O─Si or O─Si─O) versus bond stretching for strain accommodation in amorphous silica. Therefore, in our experiments, once the siloxane bond between tip and surface has been formed, the deformation process resulting from the slide-hold-slide protocol will first cause and angle distortion of the SiO 4 tetrahedron, as schematically represented in Fig. 3(a) gray inset. It would be interesting to compare this soft deformation, with a torsional stiffness k θ in the order of 10 −19 J=rad 2 [44], corresponding to a lateral stiffness k 1 ¼ k θ =R 2 ∼ 1 N=m (where R ¼ 0.35 nm is the length of an unstretched siloxane bridge), to the values of k presented in Fig. 2. Although the initial values at high temperature (300 K) are in the same order of k 1 , a direct comparison is not possible without taking into account the complex network of chemical bonds (here represented as springs) of the contacting objects in the proximity of the contact area. In particular, it would be wrong to estimate the number of contact points in the experiment simply as N ¼ k=k 1 . It would mean that the contacting springs are rigidly arranged in parallel, but this brutal assumption underestimates N. A detailed analysis of this problem is left for future work.
Building on the success of MD in describing the mechanical properties of a-SiO 2 [37,38,[40][41][42][43] while providing the enlightening atomic-scale picture of such processes, here we follow these works in order to obtain the atomic detail of the thermal activated contact aging process. We considered the indentation process of a realistically large a-SiO 2 tip or surface composed by 352137 atoms with an experimentally matched roughness [45] and let them evolve under a constant loading force F N for long periods of time (see Appendixes C, D, and E and Fig. 3 for details). Although the well-known limitations of fully atomistic methods applied to tribological processes [12,13,[46][47][48]  . The contact stiffness shows logarithmic aging, and fits (solid black lines) for hold times up to t hold ¼ 0.6 s reveal temperature-dependent slopes, as shown for the native oxide layer (c) and the 300nm-thick oxide layer (d). In both cases, the temperature dependence of the slopes can be well approximated by a linear function with almost identical slope. (e),(f) Initial contact stiffness at t ¼ 10 ms as a function of the temperature. In both cases, the initial contact stiffness decreases with temperature. restrict us to much smaller timescales than experimental ones, these are still long enough to appreciate the logarithmic time dependence of the bond formation process expected from the experimental results [12,13]. The results in Fig. 3(d) show the number N of bonds formed between tip and surface as a function of the contact time t. In the cyan boxed inset of Fig. 3(a) we show one representative siloxane bond occurring at the tip-surface interface. For all temperatures T considered, Nðt; TÞ increases logarithmically with time as N 0 ðTÞ þ αðTÞ lnðt=t 0 Þ, where α represents the aging rate and N 0 ðTÞ corresponds to the atomic bonds formed before the logarithmic aging is observed [see Fig. 3(d) and Figs. 3(b) and 3(c) for spacial distribution of these bonds at different temperatures]. More importantly, as shown in Fig. 3(e), the slope α is found to increase linearly with T in agreement with our experimental observations of the contact stiffness. Therefore, MD simulations not only support that contact aging occurring in dry silica nanocontacts emerges from the gradual formation of siloxane at the contact, in agreement with previous works [12,13,36], but also that an increasing temperature results in a faster aging (bond formation) rate, as shown in Figs. 3(b) and 3(c). Furthermore, quite surprisingly the complex bond dynamics occurring at the interface of such large contact is well described by the simple analytical model by Liu and coworkers [12,13] based on drastic assumptions, such as a homogeneous distribution of energy barriers in a welldefined energy interval ranging from 0 to E 0 . Indeed, in the intervals between 20 and 80 ns that are here considered, bonds are repeatedly formed and broken at different sites in the contact region, with a complex interdependent energy barrier distribution. Therefore, both our simulations and experiments provide support for the nontrivial linear temperature dependence of α proposed by Liu and Szlufarska [12] and validate the conclusions of this simplified model well beyond its restrictive assumptions. Furthermore, the very similar rates obtained at different loads aging rate does not depend on F N [whereas the initial contact area or number of bonds is considerably larger when the load is increased, see Fig. 3(f)]. Finally, we address the observation that in FFM measurements the initial value of the contact stiffness decreases with the temperature, as shown in Figs. 2(e) and 2(f). According to the aforementioned atomic bond picture, this seems to imply that the higher the temperature, the lower the number of initial bonds in the contact. Given that bond formation is a thermally activated process, this result looks, at first, surprising. In fact, our MD simulations at F N ¼ 40 nN predict a small increase in the initial number of bonds with the temperature (∼7% in the range of temperatures used in the experiments). Nevertheless, it can be understood by considering how the contact is prepared experimentally. Since the tip is halted while sliding, the interface starts aging under finite shear stress conditions which are highly temperature dependent. As shown in Fig. 1(d), the sliding friction decreases with increasing temperature. As a result, the higher the temperature, the lower the shear stress in the contact. Comparing the MD simulation at different loads we observe a similar behavior: i.e., lower loads, or equivalently lower contact stress, correspond to a smaller initial number of bonds [see Fig. 3(f)] that result in a smaller initial contact stiffness, in agreement with the experiments. Most importantly, the simulations clearly show that the rate at which the number of bonds increases with the temperature, i.e., α, does not depend on the contact stress [see Fig. 3(e)] for the range of loads employed here. Moreover, these findings are in agreement with Eyring's law [49] according to which this stress shifts the energy barrier distribution for bond formation, a mechanism that was also confirmed to apply for nanocontacts in several previous studies [14,50]. In our case, the shear stress allows us to form a larger number of bonds during the first 10 ms of hold time, while the subsequent bond formation rate related to bonds of higher energy barriers remains unchanged by the shear stress.
Overall, our findings highlight not only the importance of qualitative aging (contact strength) with respect to quantitative aging (contact area), but also shed light on its atomic origins. Guided by its atomic origins, one may expect that our results go well beyond the case of simple silica contacts. In fact, it is sensible to assume that in any homocontact of covalently bonded materials (e.g., diamondlike carbon, SiC, etc.) the naturally present dangling bonds, defects, and step edges naturally present will provide a route to a chemically induced qualitative contact aging analogous to SiO 2 contacts. It is worth mentioning the recent experiments by Weber et al. [51] that revealed a surprising increase of the real contact area [52,53] (quantitative aging) during macroscopic slip events, which are accompanied by a reduction in friction. Then it follows that the smaller shear stress per unit area (qualitative aging)-directly related to the contact stiffness shown in Fig. 2-underpins the weakening of dynamic friction with respect to the static one. This suggests that a synergetic combination of single asperity slide-hold-slide experiments with real contact area measurements [51][52][53] allows us to further explore the thus far elusive qualitative aging from idealized to realistic conditions, thus providing new insight into our understanding of third-body chemistry and how this process governs aging processes at nanoscales and, consequently, macroscales.

III. CONCLUSION
In summary, the temperature-dependent measurements presented in this paper show that contact aging of silica interfaces in the absence of water is observed on the timescale of up to several seconds, and is a thermally activated process. This is demonstrated by first noticing that the lateral contact stiffness k is a more suitable marker of contact aging, as compared to the static friction. Being proportional to the number of interfacial bonds, k scales logarithmically with the contact time, and its rate scales linearly with temperature, in full agreement with the accompanying MD simulations of the contact evolution. Our results also shed light on the distinction between contact "quantity" and "quality," which has been recently proposed in the literature. The average number of interfacial bonds seems to be all that is needed to characterize the qualitative aging of a contact between two surfaces in dry conditions. Finally, our experiments also highlight that faster contact aging does not necessarily lead to a faster increase of static friction, an effect which can be attributed to the often overlooked role of rupture in static friction measurements.

APPENDIX A: EXPERIMENTAL METHODS
All low-temperature AFM experiments have been performed using a standard Omicron VT-AFM system at pressures of p ≤ 5 × 10 −10 mbar and Si cantilevers with a normal force constant of k ¼ 0.2 N=m (Nanosensors, PPP-LFMR). The lateral forces were calibrated using the calibration technique introduced by Bilas et al. [54]. To remove residual water and hydrocarbons, all tips and samples have been cleaned by thermal treatment at a temperature of about 250°C for 3 h under UHV conditions prior to the AFM measurements.
Our contact aging measurements are based on an experimental approach similar to the slide-hold-slide protocol introduced by Li et al. [3]. We recorded conventional friction loops of 300 nm width, where scanning was suspended for a defined time t hold in the middle between the two cantilever turning points to initiate contact aging. For both static and sliding friction, the required reference level of the lateral force signal is extracted from the full friction loops. A relatively high scan speed of 2.5 μm=s has been used throughout all of our experiments, primarily to minimize effects of gradual and thermally activated contact breaking. For measurements with systematic variation of the hold time at a constant interface temperature, one line of 300 nm width was repeatedly scanned with a nominal normal load of 1.2 nN while varying the hold time. Starting from t hold ¼ 10 ms the hold time was steadily increased up to 2.5 s and subsequently decreased back to 10 ms. No significant differences were found between these two ramps. After each cycle in hold time, the tip was moved to a new sample area and the procedure was repeated, resulting in 40 positions each measured for 20 hold times between 10 ms and 2.5 s.
In the slide-hold-slide experiments, the effective contact stiffness k between the tip and the sample is evaluated from the slope of the lateral force buildup that occurs directly when sliding is reinitiated. Only the first 30% of total height of any friction peak are used for contact stiffness determination to avoid including data points already within the onset of rupture.
Before and after the slide-hold-slide experiments, conventional AFM images larger than the sample area are taken to ascertain a clean surface area and to evaluate the occurrence of any alterations induced by the experiment. Temperatures below 150-200 K occasionally showed wear and were thus altogether excluded from the contact stiffness evaluation. To ensure that no systematic influences related to tip degradation affect the experiments, the temperature sequence has been randomized, and frequent reference point measurements at room temperature turned out to be highly reproducible.

APPENDIX B: CANTILEVER CREEP EFFECTS
In Fig. 1, a sharp decrease of the lateral force was observed at x ¼ 0, i.e., after the hold time. This observation is universal for all slide-hold-slide curves and can be attributed to a creeping relaxation movement of the tip during the hold time. This creep motion becomes understandable, when considering that the cantilever torsion measured as sliding friction is exactly the force to sustain cantilever sliding. This means that when the tip is halted, a lateral force equal to the sliding friction initially acts on the interface between tip and sample. This force, in combination with thermal activation, leads to a motion of the tip in scanning direction and accordingly to a reduction of the shear stress at the interface with time. From Eq. (B1), the time-dependent velocity of this creep motion can be estimated based on the lateral force values F min recorded directly after each hold time: where F min;i and F min;i−1 are the minimum lateral forces after the hold times t hold;i and t hold;i−1 , respectively. Additionally, k x denotes the effective lateral stiffness of the cantilever. The resulting creep velocities as a function of the hold time are shown in Fig. 4 for a typical set of slide-hold-slide experiments and realistic values of k x . Velocities down to the picometer per second range can be resolved by this approach. A finite velocity is observed at any given hold time, yet the velocity rapidly decreases with time. After 50 ms, the creep velocity has already dropped to approximately 1 nm=s, which is 3 orders of magnitude slower than the sliding velocity. For the example shown in Fig. 4, the total slip distances during the maximum hold time of 2 s vary between 0.7 nm (300 K) and 1.4 nm (200 K), with typically around 50% of the creep already occurring during the first 15 ms. We thus conclude that creep of the tip can affect the aging behavior only during the first milliseconds of the hold time, but can be neglected for extended hold times.

APPENDIX C: MD SIMULATION DETAILS
The MD simulations were carried out using LAMMPS [55]. The equations of motion are integrated with velocity-Verlet algorithm [56] with a time step of 0.001 ps, which correctly samples the dynamics of the lightest atoms, i.e., oxygen, at high temperatures or pressures [57,58]. The coordinates were saved every picosecond. The temperature and pressure are regulated via a Nose-Hoover [59][60][61][62] thermostat and barostat, respectively. In the preparatory annealing-quenching cycles, NPT (constant temperature and pressure) ensemble simulations were performed, whereas all other simulations were performed in the NVT (constant volume and temperature) ensemble.

APPENDIX D: SiO 2 ATOMIC MODELS (FORCE FIELDS) AND THEIR SUITABILITY FOR INDENTATION STUDIES
Si─Si, Si─O, O─O interactions were modeled via a bond-order three-body Tersoff [63,64] potential parametrized by Munetoh et al. [63] to describe a Si, O mixed system. Munetoh [63] parametrization is based on experimental data of α-quartz as well as ab initio calculations of small fragments. This parametrization provided results in excellent agreement with experimental data [63]: namely, the structural properties and energies of the various crystalline phases as well as the amorphous phase for silica [57,63], radial distribution functions, and phonon density of states of silica glass [63]. Furthermore, it has also been experimentally validated to describe atomic rearrangements in silica glass surfaces induced using aberration-corrected transmission electron microscopy [37], transformation mechanisms [65] and thermal conductivity [65,66] in porous silica, and also the Hugoniot curves [58,67] describing the dynamic shock-induced response of silica. Additionally, a good agreement with experimental results has also been found in other mechanical properties such as the deformation behavior of amorphous silica upon indentation [40,41].

APPENDIX E: SIMULATION PROTOCOL
The simulation protocol was composed by two stages: a preparatory stage, where we generate the initial coordinates of the tip and the surface, and a second stage, where we indent the surface with a constant normal load and temperature. In the first stage, we built an amorphous SiO 2 surface of 60 × 60 × 18 nm 3 composed by 3345408 atoms by using the atomic configurations of amorphous SiO 2 cubes of ð5Þ 3 nm 3 reported in Ref. [68]. Then, in order to remove the superperiodicity of 5 nm, we improved SiO 2 amorphization by following the standard melting-quenching procedure [2, 12,13,40,41,57,58,63,[65][66][67]69]. This protocol consisted of melting SiO 2 by heating it up to 5000 K in 250 ps at zero pressure followed by a quenching to 300 K with a cooling rate of 11 K=ps. Note that during the annealing stage the temperature was raised up to 5000 K, which is well above the melting point of SiO 2 obtained experimentally [70]. Although already at temperatures of 1900-2400 K SiO 2 started to melt, in agreement with previous works [57], annealing to higher temperature promotes a sufficiently fast diffusion of the atoms so that we could produce an equilibrated liquid in the short time span (as compared to timescales used in the experiments) of the simulation. Additionally, as in Ref.
[2], by using freeboundary conditions in the out-of-plane direction, the obtained surface displays a roughness on the top and bottom part of the surface. Moreover, the surface stress of the asprepared surface is allowed to relax during the ∼0.5 ns quenching stage. Then, we extract a slab of 25 × 25 nm 2 and a thickness of 8.5 nm. Similarly, we obtain the tip atomic coordinates by selecting all atoms within a sphere with a radius of 10.25 nm and then we kept only the bottommost atoms 8.7 nm. In order to avoid any synchronization on the roughness of the surface and the tip, the tip roughness is provided through the bottom part of the larger slab, while the roughness of the surface is provided through the top part of the larger slab. In total, this system composed by the tip and surface contains a total of 352137 atoms and the charge neutrality was assured throughout the whole process by imposing the ratio SiO 2 . The density of the amorphous silica thus created is of 2.22 g=cm 3 , in agreement with the experimental value [71] of 2.20 g=cm 3 and previous MD simulations [40,41,57,58,63,[65][66][67]69,71]. The amorphization of the sample has been confirmed by calculating the radial distribution functions of Si─Si, Si─O, and O─O pairs whose first peaks occur, respectively, at 3.13, 1.65 and 2.67 Å, in agreement with previous MD simulations as well as with experimental results [40,41,57,58,63,[65][66][67]69,71]. The root-mean-square roughness of the surface thus obtained is of 1.56 Å, which is similar to a previously reported MD simulation whose protocol was followed here [2] and close to the experimental value [45]. In the second stage, the position of the bottom atoms of the surface (bottom 0.7 nm, shown in green in Fig. 1) is kept fixed during the indentation process. Initially the tip is brought into contact with the surface while assuring no atomic bonds exist at t ¼ 0, and then the indentation starts by applying a constant normal force of 40 nN (with the exception of the 80-nN simulation) on the tip's top atoms (top 0.7 nm, shown in green in Fig. 1). All the other atoms are allowed to move freely during the indentation process. When the contact is formed we follow the number of bonds formed between the tip and the surface atoms. The criteria for stable bond formation at the interface correspond to Si─O or O─O pairs of atoms that are within the maximum equilibrium distance reached due to solely thermal fluctuations at room temperature, i.e., 1.7 and 1.55 Å, respectively. Additionally, these distances are close to the average distance defined in the force field [57,63] and also agree well with the experimental distances of these bonds [74]. By recording the coordinates of all atoms every picosecond, we then compute the average number of bonds between the tip and the surface every 50 ps. The results thus obtained are shown in Fig. 3.