Time-Dependent Electrical Contact Resistance at the Nanoscale

Conductive-atomic force microscopy (C-AFM) and molecular dynamics (MD) simulations are used to investigate time-dependent electrical contact resistance (ECR) at the nanoscale. ECR is shown to decrease over time as measured using C-AFM and estimated using two approaches from MD simulations, although the experiments and simulations explore different time scales. The simulations show that time dependence of ECR is attributable to an increase in real contact area due to atoms diffusing into the contact. This diffusion-based aging is found to be a thermally activated process that depends on the local contact pressure. The results demonstrate that contact aging, previously identified as an important mechanism for friction, can significantly affect electrical conduction at the nanoscale.


Introduction
Electrical contact resistance (ECR) plays a key role in the function and performance of a variety of electromechanical components used in a wide range of applications, including switches/relays [1][2][3][4][5][6][7] and electronic connectors [8,9]. Beyond the conventional engineering scale, ECR is particularly important for electromechanical components at the micrometer scale, due to their large surface-to-volume ratio. At such small length scales, it has been found that ECR is affected by multiple parameters and conditions, including contact pressure, temperature, surface properties, and the environment [1,3]. Applications that rely on ECR for proper function must ensure that these parameters are tightly controlled during operation.
Another factor that affects the reliability of electrical connections, but cannot be controlled, is time. Time-dependent ECR at the microscale has been extensively studied using radio frequency microelectromechanical system (RF-MEMS) switches [2, [5][6][7]. These devices have very short contact times (i.e., very high switching frequency), so the 1 3 50 Page 2 of 9 change in ECR with time is an important parameter affecting overall performance [10].
There are multiple mechanisms that potentially affect time dependence, and some cause an increase in ECR with time, while others cause a decrease. The chemical character of the contact may gradually degrade due to oxidation or fretting corrosion, resulting in higher ECR [8,11,12]. Alternatively, adsorbed chemical species can be removed from the contact by stress [13], which would cause ECR to decrease over time. Also, an increase in contact area due to asperity deformation (i.e., creep) can cause a decrease of ECR with time. Creep can be coupled with Joule heating that increases ECR [6], but, at room temperature, the increasing contact area is dominant such that decreasing ECR over time is typically observed.
Creep is described by a power law relationship where resistance R decreases with time t rapidly initially and then approaches a final limiting value. To capture this behavior, time-dependent ECR at the microscale can be expressed as R(t) = At − + B [2, 6,14]. Physical interpretations of each term in this relationship were proposed based on single asperity contact and strain rate formulas that assume creep with a power law dependence near room temperature [2]. In this context, A is a parameter that is related to interfacial geometry and roughness, is determined by the creep behavior of the asperity material, and B is a theoretical limiting resistance at t → ∞ . Despite the frequent use of this model and its ability to generally capture trends observed in RF-MEMS experiments, the physical mechanisms underlying resistance change over time are not fully understood.
At the microscale, where multiple asperities are involved in contact, it is often difficult to correlate changes that occur locally at each asperity with time-dependent ECR that is measured globally. Alternatively, single asperity contacts at the nanoscale can be studied to isolate and understand the physical mechanisms underlying the change of ECR with time. At the nanoscale, the contact can evolve over time through processes collectively referred to as qualitative and quantitative contact aging [15,16]. Contact quality is the contact strength per unit area and is subject to time evolution at nanoscale by formation or breaking of bonds at the interface. Contact quantity is the size of the contact. Contact size can increase with time due to creep, i.e., time-dependent plastic deformation under a constant applied load where the growth rate depends on temperature and the mechanical properties of the contact [15,16].
Contact aging has recently been investigated by assessing friction variations for silica-silica [16] and Si-highly oriented pyrolytic graphite (HOPG) contacts [17] as well as metallic nanoparticles on HOPG [18]. However, the same physical mechanisms that lead to friction variation can also affect ECR, since resistance decreases when there are more conduction pathways across the contact. The number of conduction pathways across a given contact can increase due to an increase in either the real or apparent contact size (or both). Apparent (also called nominal) contact area is defined by the global dimensions of the bodies in contact, while the real (also called actual or true) contact area is defined by the region where mechanical contact is achieved via mating asperities of the surfaces [19,20].
Despite extensive studies of contact aging at the nanoscale and time-dependent ECR at the microscale, time-dependent ECR has remained relatively unexplored at the nanoscale. Conduction at single asperity contacts can be measured experimentally using conductive atomic force microscopy (C-AFM). C-AFM has been used to study ECR of various material systems including thin SiO 2 films grown on doped Si substrates [21], HOPG [22,23], carbon nanotubes [24,25], self-assembled monolayers of organic molecules on gold [26], and gold islands on HOPG [27]. However, such studies have not investigated the time dependence of ECR at the nanoscale.
While C-AFM can be used to characterize the nanoscale electrical properties of a sample, it cannot provide atomistic detail about the structure and chemistry of the "buried" contact between the nanoscale tip and the sample surface, which hinders interpretation of observed phenomena. Therefore, simulations of single asperities need to be used to model C-AFM. Ab initio calculations can capture electronic effects and, therefore, have been used to calculate current across nanoscale interfaces directly [28][29][30]. However, such simulations are typically limited to hundreds of atoms and so cannot capture contacts more than several atoms wide. Further, the inability to model a significant amount of material away from the contact means that such methods cannot incorporate material deformation in response to load.
The size scale issue associated with ab initio approaches can be addressed by molecular dynamics (MD) simulations that have been widely used to model single asperity contact [31,32]. Since classical MD does not explicitly model electrons and cannot directly calculate conduction, previous MD-based studies of single asperity contact conduction have assumed that current is proportional to contact area [33,34]. In this approach, contact area is obtained by calculating the number of atoms in contact and then conduction is calculated from a classical relationship between current and area. Alternatively, current can be approximated from reactive MD simulations using a method called EChemDID [35], that enables calculation of current from the electrochemical potential of atoms. This approach has been used to study conduction across nanoscale contacts previously [23,27,36,37]. Lastly, the relationship between atom-atom distance and conduction from first principles calculations can be fit to an empirical model which is then used with MD atom distances to estimate conduction [38]. However, none of these studies explored the change of ECR with time.
Here, we investigated time-dependent ECR at the single asperity level using C-AFM experiments and MD simulations with two approximations for calculating current. Both experiments and simulations showed that resistance decreased over time, albeit at different time scales. We hypothesized that this trend could be explained by a contact aging mechanism [16,18,39]. In C-AFM, contact area was correlated with measured resistance using expressions for ballistic electron transport. Both real and apparent contact area were calculated directly from atom positions in the simulations. The contact areas from experiments and simulations were fitted to an exponential equation for contact aging and time constants were extracted. Finally, the simulations were used to identify the atomic scale origins of contact conduction aging on short time scales that may also suggest mechanisms for longer time scales.

Experiments
The C-AFM measurements were performed on ZYBquality HOPG substrates (TedPella). For each experiment, an HOPG substrate was freshly cleaved via scotch tape and immediately placed inside the sample chamber of a commercial AFM instrument (Asylum Research, Cypher VRS). A constant flow of 99.999% purity nitrogen (N 2 ) gas through the sample chamber was maintained overnight to purge ambient gases out. This resulted in an inert N 2 environment inside the sample chamber and limited adsorption of contaminants onto the HOPG surface during the experiments. Under the N 2 environment, the HOPG substrate was annealed at 100 °C for 60 min to desorb adsorbates from the sample surface. Following the annealing procedure, the HOPG substrate was cooled down to a temperature of 38 °C, which was maintained stable (± 0.1 °C) throughout the measurements via PID control. This temperature fluctuation might lead to some level of lateral drift between tip and sample. However, as our experiments were conducted on atomically flat terraces of the HOPG sample, away from step edges and other defects, effects from thermal drift in the experiments should be negligible. Either conductive diamond tips ( CDT − CONTR, Nanosensors ) or Pt tips (HQ:NSC18∕Pt, masch ) were used for the C-AFM measurements. The radii of the diamond and Pt tips used were 120 and 30 nm (per manufacturer specifications), respectively. Prior to electrical characterization, the adhesion force between tip and sample was measured via force-distance (F-d) spectroscopy.
Subsequently, current-voltage (I-V) spectroscopy was performed, whereby the applied bias voltage ranged from − 100 to + 100 mV. During I-V spectroscopy, either no normal load (for conductive diamond tips) or a small normal load (5 nN, for Pt tips) was applied in addition to the adhesion force between the tip and the sample (which was in the range of 20 to 30 nN for conductive diamond tips and around 10 nN for Pt tips). I-V spectroscopy measurements were performed over varying periods of time (ranging from a few minutes to a few tens of minutes). The resistance values were extracted from the linear region of each I-V curve around zero bias. Note that the frequency at which I-V curves were taken did not affect the general trend in the results. Additionally, I-V curves are separated by durations of ~ 2 s or more, which should be long enough for most transient electromagnetic effects to dissipate. Following I-V spectroscopy experiments, the adhesion force between tip and sample was measured once again to verify that there were no significant tip changes during measurements. Finally, the area around the I-V spectroscopy site on the HOPG surface was scanned to confirm that there was no damage on the HOPG surface.

Simulations
MD simulations were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS) package [40] and atomistic visualizations were obtained using the OVITO package [41].
The model system was designed based on previous simulations of contact aging that focused on the phenomena for SiO 2 [15,16,42]. Here, we model a hemispherical amorphous SiO 2 tip apex with 9 nm radius contacting a four-layer HOPG substrate, as shown in Fig. 1. The amorphous SiO 2 structure was generated following the approach described in Ref. [16]. The atoms in the topmost few angstroms of the tip were treated as a rigid body, while the bottommost layer of graphene in HOPG was held fixed. A Langevin thermostat was applied to all the non-constrained atoms to control temperature throughout the simulation. The interatomic interactions within the HOPG substrate and the SiO 2 tip were described using the AIREBO [43] and Tersoff [44] potentials, respectively. The Lennard-Jones 12-6 potential model was used to describe the interactiona between tip and substrate with the parameters taken from Ref. [45]. Dynamic simulations were run with a 1 fs time step. The Fig. 1 Schematic of C-AFM on HOPG. Inset is a snapshot of the atomistic simulation tip was brought into contact with the HOPG substrate and a normal load of 80 nN was applied. Simulations were run at six different temperatures ranging from 350 to 650 K for a duration of at least 10 ns.
To calculate the contact area, we first identified the atoms that were "in contact" using a distance criterion. Specifically, if the distance between a tip atom and its nearest substrate atom was less than 0.4 nm [36,46], then this tip atom was considered to be in contact with the substrate [19]. Then, the real contact area was calculated from the number of contact atoms multiplied by an approximate atomic area of 0.021 nm 2 . This value was estimated as the area of a circle with diameter 0.165 nm, the distance to the first peak of the radial distribution function for SiO 2 . Lastly, the apparent contact area was calculated by measuring the area of a polygonal region (a convex hull) enclosing the contact atoms [47].

Conduction Estimation
Classical molecular dynamics simulations do not model electrons, but electrical conduction can be approximated based on the positions of atoms in the model. Here, we used two approaches: EChemDID [37] and an analytical model for conduction as a function of atom distance based on density functional theory (DFT) data [38]. Both calculations were performed using atomic configurations that were extracted from different time instances of the contact simulations described above. Note that these methods estimate conduction based only on atom positions and the material itself, in this case SiO 2 , does not affect the result. So, although SiO 2 is not conductive, the conduction results can be generalized to any amorphous material.
The first approach to calculating current was EChemDID. Since EChemDID requires reactive potentials to describe the current flowing through the system, the atomic interaction models were replaced with ReaxFF with parameters taken from Ref. [37]. Next, the topmost atoms of the tip and the bottommost layer of HOPG were coupled to an external voltage bias and the relative current was obtained from the combination of Ohm's law and the continuity equation under the assumption of diffusive transport and no Joule heating [35].
The second approach to estimating conduction relied on empirical models for the relationship between atom-atom conduction and the distance between atoms [38]. Briefly, conduction channels were identified based on the distance between C atoms in the topmost graphene layer and Si atoms in the SiO 2 tip. The conduction for each channel was then calculated using empirical models developed by fitting DFT data for conduction between individual Pt atoms to exponential decay functions. This approach captures both constant direct metallic conduction and distance-dependent tunneling conduction. The approach also models the effect of O atoms to decrease conduction, as expected to be the case for any metallic AFM probe containing a native oxide [48]. Finally, conduction across the contact was calculated as the sum of the conduction through individual atom-atom channels. Note that the empirical functions used in this calculation were developed for Pt atoms, so the results can only be used to provide general trends for the model system being studied here.

Results
A representative resistance measurement from the C-AFM experiments performed with a conductive diamond tip is shown in Fig. 2a. Over a period of about six minutes, the resistance decreased by approximately a factor of two. This decreasing trend was observed in other experiments (see Fig. S1), performed with two different tip materials (conductive diamond and platinum) and at three different temperatures (305 K, 311 K and 323 K).
Resistances estimated from simulations at 650 K using EChemDID and an empirical conduction vs. atom distance relationship are shown in Fig. 2b and c, respectively. Note that Fig. 2 Resistance as a function of time from a C-AFM experiments and MD simulations with conduction estimated using b EChemDID and c empirical models for conduction-atom distance relationships the unit of the resistance is arbitrary since it only reflects the change in atom positions over time, and does not depend on the involved materials. Comparing the two conduction estimation methods, the resistance estimated from the empirical conduction-atom distance relationship (Fig. 2c) decreases faster than that calculated using EChemDID (Fig. 2b). This difference is attributable to the fact that the calculation methods are based on intrinsically different approaches: EChemDID couples the system to a voltage bias and obtains the current flowing through the system by describing the equilibration of electrochemical potential throughout the system, while the other approach estimates the conduction between the tip and the substrate by identifying atomic channels through the contact and evaluates the conduction at each channel using a DFT-derived current-separation relationship.
Regardless, resistance estimated from the simulations using both methods decreases with time, which is qualitatively similar to the trend observed in experiments. However, there is a smaller change in ECR and the change occurs over a much shorter time scale in simulations than experiments. The difference between experiments and simulations may originate from the difference in size and time scales, as discussed later.
The decrease in resistance with time shown in Fig. 2 is probably not due to mechanisms that involve the presence (and subsequent stress-driven removal) of insulating species. The presence of insulating layers (in the form of adsorbed organic molecules from the environment and/or oxide layers) at the tip apex is expected to lead to non-Ohmic behavior at the tip-sample junction, which was not observed in the experiments. Additionally, the simulations, that were performed in perfect vacuum and therefore do not include insulating layers, also report a gradual drop in ECR. Therefore, it is reasonable to propose that the observed trend is due to a change in contact area. Specifically, an increase in contact area over time will increase the number of conduction channels and therefore decrease ECR.
In diffusion-based contact aging, when the tip is initially brought in contact with the substrate, the load is supported by the first atoms to come into contact. To achieve pressure equilibrium, additional atoms diffuse to the contact region such that the contact area increases and the local stress decreases. This process continues until the interface can support the stress without further migration of the atoms into the contact region. In this context, it has been proposed that the contact area A increases with time t and then saturates following [39]: where A 0 corresponds to initial contact area, is a constant that defines the maximum relative increase of A(t) and is a characteristic saturation time. (1) To analyze the experimental data in the framework of contact aging, the measured ECR needs to be converted to contact area. Since the contact radii formed between the tips and HOPG sample were small (2-5 nm, determined using the Derjaguin-Muller-Toporov (DMT) model [49]) compared to the electron mean free path in HOPG (200 nm [50]), ballistic transport of electrons can be assumed. Measured resistance R was converted to area using Sharvin's equation for ballistic electron conduction [51]: where is the resistivity and is the mean free path of the electrons. The values of and used to convert ECR values obtained with the conductive diamond tip to contact areas were 95 × 10 −6 m and 100.7 nm, respectively, whereas those in case of the Pt tip were 2 × 10 −5 m and 104.7 nm, respectively. These values were calculated from the arithmetic mean of the properties of the two materials in contact [52], either diamond and HOPG or Pt and HOPG.
Six independent data sets were fit to Eqs. 1 and 2, as shown in Fig. S1 with the resultant fit parameters reported in Table S1 in the Supplementary Information. The electrical contact areas obtained from this calculation range from 0.009 to 0.296 nm 2 . The variation from test to test reflects differences in tip size and shape. We also calculated apparent contact area using experimentally measured adhesion forces with the DMT model. The apparent contact areas were found to be on the order of a few tens of nm 2 , two to three orders of magnitude larger than the electrical contact areas extracted by fitting ECR data. The discrepancy between the mechanical and electrical contact areas is attributed to the fact that the portion of tip-sample contact that is electrically conducting is much smaller than the physical contact area [21,53].
The characteristic time scale of ECR aging was also obtained from this fitting, as reported in Table S1 in the Supplementary Information. The average characteristic time was found to be 564 s, although there is wide variation from test to test. This variation could be due to atomic mechanisms (e.g., diffusion to the tip apex) contributing to the contact aging. That is, the local atomic arrangement of the tip apex may be responsible for different time constants since some atomic structures can reach equilibrium faster than others.
The contact size in simulations was calculated directly from the positions of the atoms in the interface. Contact size was quantified as either the real or apparent area, where the former was obtained from the number of atoms in contact multiplied by a per-atom area and the latter was calculated from the area defined by the perimeter of the contact atoms. The real and apparent contact areas calculated from the simulation data are shown in Fig. 3. As expected, the apparent contact area was larger than the real contact area [19,20]. Also, the real contact area increased with time, whereas the apparent contact area remained constant (see Fig. S2 for atomistic visualization of the contact at 0 and 2.2 ns). The observed real contact area increase corresponded to the migration of almost 30 new atoms into the contact. The real contact area data were fit to Eq. 1 and the characteristic time constant was found to be 0.55±0.05 ns. The characteristic time constant in simulations is many orders of magnitude smaller than that in the experiments (564 s vs. 0.6 ns). However, the duration of the simulations is also about ten orders of magnitude shorter than the experiments (~ 1000 s vs. 10 ns). Also, the contact size in simulations is much smaller, such that there are far fewer atoms available to participate in a diffusion-based process. Therefore, we cannot directly compare experiments and simulations. Regardless, we assume that simulations are capturing the behavior of a conductive junction in the initial tens of nanoseconds of contact and use the simulations to investigate the atomistic origins of the aging behavior.
The results in Fig. 3 suggest that the contact aging in the simulations is primarily attributable to diffusion of atoms within the original perimeter of the contact. This is corroborated by the snapshots in Fig. S2 that show new contact atoms are observed most within the original contact perimeter. To confirm this, we analyzed the radial density of atoms in the contact. This was quantified as the change in local density between 0 and 10 ns at different radial positions. As shown in Fig. 4, the atom density increases with time in the middle of the contact (radial positions less than about 1.5 nm) but there is little change in the density near the perimeter of the contact.
In diffusion-based contact aging, atoms diffuse to the contact region such that the contact area increases and the local stresses decrease [39]. We have shown that the atom diffusion occurs primarily at the center of the contact (Fig. 4). This can be explained by the fact that the local contact pressure is higher in the center of the contact [19,54].
To confirm this in our simulations, we calculated the local pressure in the contact using the DMT contact equation. The local contact pressure is shown on the top x-axis of Fig. 4. At the perimeter of the contact the pressure is zero and the highest pressure is at the center of the contact. As shown in Fig. 4, the change in density increases with increasing pressure.
The results so far suggest that ECR aging is attributable to diffusion of atoms within the contact area over time. Diffusion is a thermally activated process, so the characteristic aging time should be described by an Arrhenius law as: In this expression, U is the attrition barrier, T is the temperature and k b is the Boltzmann constant [39]. To test the suggestion that ECR aging is thermally activated, simulations were run at temperatures ranging from 350 to 650 K. Note that simulations at lower temperatures were run longer than 10 ns to ensure the system reached steady state. For each case, the real contact area as a function of time was fit to Eq. 1. The data and fit lines for all cases are shown in Fig. S3, with fit parameters in Table S2 in the Supplementary Information. The best fit characteristic time scales are plotted as a function of inverse temperature in Fig. 5. This data were then fit to Eq. 3, as shown by the dashed line in Fig. 5. The fit to Eq. 3 resulted in an attrition barrier of U = 182 ± 41 meV . Attrition-based phenomena take place by activating (i.e., atom hopping) over an energy barrier that, for atomic friction, can be assumed to be on the same order as the energy needed for contact separation (and then slip) during stickslip sliding [39]. The barrier height calculated here is comparable to previously reported activation barriers for sliding on HOPG surfaces in vacuum (100-150 meV) that were determined using a Prandtl-Tomlinson model to interpret atomic stick-slip friction data [55,56]. The observation that the simulation trend is consistent with the Arrhenius model indicates that ECR aging is indeed a thermally activated process.

Conclusion
ECR at the nanoscale was shown to decrease with time using C-AFM experiments and MD simulations with two approximations for conduction. It was hypothesized that contact aging, previously shown to affect nanoscale friction, may also explain time-dependent ECR. To test this, contact area was calculated from both simulations and experiments and the results were found to be described well by an expression for contact aging. Then, the simulations showed that the increase in contact area (and decrease in ECR) was attributable to diffusion of atoms into the central part of the contact, i.e., due to an increase in real contact area. The greater diffusion into the center of the contact was correlated with higher local pressures predicted by continuum contact mechanics. Lastly, simulations run at different temperatures revealed ECR aging to be a thermally activated process.
It is important to note that similar trends were observed for all experiments and simulations regardless of the tip material, indicating that contact aging may play a role in ECR generally, although material-dependent parameters like attrition barrier and environmental parameters such as temperature will determine the magnitude of its effect. This finding is particularly relevant to small-scale devices for which ECR affects device performance. It also has scientific implications since it was shown that contact aging affects not only friction, as previously observed, but also electronic conduction. This suggests that aging is a more general phenomenon that may affect various processes that occur at nanoscale contacts.