Nanocontacts and Gaussian Filters

The inherent surface roughness at atomic scale contacts has profound effects on their local stresses, contact areas, and so on, which yield deviations from the predictions of continuum models. While these deviations can be interpreted as a breakdown of the models, attempts have been made to reconcile the two via different paths, e.g., proposing new methods for estimating the contact area, or filtering the contact stress fields. The applicability of the latter to randomly rough nanocontacts, however, has not yet been explored hitherto. Here, we show that this new method should be employed with care. Investigating two smooth and rough contacts, we find that the filtering method can make the stresses resemble each other by varying the filter’s resolution. Our results demonstrate that the filtering method improves comparisons between atomistic contacts and Hertz theory, while accounting for roughness in the rough contact results in more accurate solutions.


Introduction
From an exhaustive investigation of nanoscale contacts by means of atomistic simulations, Luan and Robbins concluded the breakdown of continuum models in describing mechanical contacts, due to the inevitable surface roughness at the atomic scale [1,2]. As a remedy, Mo et al. argued that a correct definition of contact area should be used for the post-processing of atomistic simulations, so that the results would be comparable with contact mechanics models [3,4]. They proposed that instead of measuring contact area A c by enclosing the edges of contact with a circle or a convex hull, A c should be estimated via T c A a , with T c and A a being the number of contacting atoms and the contact area of an individual atom, respectively. In order to utilize this method, the value of A a should be known. Moreover, a criterion is required for identifying the contacting atoms. This criterion may need to be defined by assuming the validity of continuum models for non-adhesive contacts [5,6], while for the post-processing of adhesive ones, such an assumption can be avoided [7,8]. The readers are referred to the literature, e.g. [5,9,10], for further details on the topic.
Very recently, Müser proposed that the results of smooth and rough nanocontacts can be compared more meaningfully to each other, if the stress fields of both are processed through the same low-pass filter [11], than if the result of a simulation with small-scale details (such as atoms, steps, etc.) is directly compared to a continuum solution of a model lacking these details. He employed two Hertzian tips (smooth and stepped), and studied their stress fields to show that applying a Gaussian filter with a resolution close to the broadest terrace of the quantized tip makes the stress distributions of both contacts indistinguishable. Hence, he concluded that a correct post-processing of atomistic and continuum solutions using appropriate filters, e.g., Gaussian ones, is required for a meaningful comparison.
This study focuses on estimating contact areas of smooth and rough nanocontacts using their filtered stress distributions. To this end, the atomistic contacts between two (smooth and rough) Hertzian rigid tips and a flat deformable substrate are examined using molecular statics. Then, the idea of smoothening the stress distributions using Gaussian filters and their effects on the final estimated interaction areas is investigated.

Model Setup
The contact is investigated by bringing a rigid indenter in contact with a deformable flat solid. The deformable solid is built from individual atoms of fcc crystal structure, with a lattice parameter of a 0 = 5.5884 Å corresponding to calcium [12]. The simulation box has lateral lengths of L x ≈ 248.4 Å and L y ≈ 253.5 Å along the crystalline directions of [ 1 10 ] and [ 11 2 ], respectively. Moreover, the originally flat solid has a thickness of L z ≈ 7.1 nm aligned with the [111] direction. The last atomic layer at the bottom of the solid is fixed, providing the needed support for the deformable solid, while periodic boundary conditions are applied along the lateral directions.
Two different rigid indenters are employed: smooth and rough. For the smooth indenter, a solid fcc structure with the lattice parameter of calcium, and with the same directions as the deformable body is generated; the thickness of the structure is selected to be three atomic layers, i.e., h = 3a 0 � √ 3 ≈ 9.68 Å . Then, the atoms of this structure are shifted along the z axis accordingly to follow the geometry of an elliptic paraboloid with a radius of curvature R = 100 Å . Finally, the height of the indenter is set to h, and the rest of the atoms are removed. The rough indenter is generated differently. First, an amorphous structure is generated using Atomsk [13], which is then subjected to a twostep cutting process: in the first step, this is truncated by a layer with thickness h, and in the second step, the positions of the atoms are compared with a Hertzian tip ( R = 100 Å ); the ones located outside of the paraboloid are removed. The rough tip is not subjected to a relaxation process as it could affect its surface roughness [14] and radius of curvature. Figure 1 shows the atomic structure of these two indenters.

Simulation Setup
Interactions between atoms in the deformable body are described by the embedded atom method (EAM) potential [15] with the database developed by Sheng et al. [12]. Moreover, non-adhesive contact was simulated by manipulating a 12-6 Lennard-Jones (LJ) potential [16], by truncating it at its minimum and shifting it to the zero energy level to ensure that the force and energy go smoothly to zero at the cutoff. The LJ parameters were selected as LJ = 0.2145 eV and LJ = 3.5927 Å [17]. The molecular statics (MS) method is performed with LAMMPS [18], using the Polak-Ribière version of the conjugate gradient algorithm [19]. The minimization technique iteratively adjusts the position of the atoms in a direction with the steepest potential energy gradient, until the system reaches its local potential energy. The reason for choosing MS over molecular dynamics is to have better comparisons between the molecular simulations and athermal continuum models, i.e., Hertz contact mechanics in this work.

Simulation and Data Collection
The contact is simulated by moving the rigid tips toward the deformable body with a step size of 0.1 Å for a total displacement of 5.8 nm and 6.8 nm for the smooth and rough tips, respectively. The total displacement is selected so as to be sure that the irreversible plastic deformation has started.
During the minimization process, the atoms' positions of the tip are fixed for the selected configuration, while the atoms of the counterpart are free to interact with each other and with the tip as before. Following the minimization, the acting forces upon each atom are measured at two stages, where the interaction between the bodies is (1) on and (2) off. The normal interacting force of each atom is calculated via F z,interact = F z,1 − F z,2 , where z refers to the normal component of force, and 1 and 2 correspond to the abovementioned states, respectively.

Binning Procedure
Once F z,interact is obtained, the xy plane of the simulation box is discretized into N 2 rectangular bins, with N being the number of bins along both the x and y axes. The area of each bin is A N = x y , with x and y being the lateral lengths of a bin along the x and y axes, respectively. The applied normal stress from an atom m on a corresponding bin j is obtained from mj = F z,m ∕A N . Then, the portion of each surrounding grid point of the bin j, i.e., j , is calculated based on the relative atom's position on the xy plane: the lateral distances between the selected atom and a grid point i of the corresponding bin Fig. 1 The smooth (left) and rough (right) tips of this investigation are calculated as x,mi and y,mi . Then, a weighting parameter is introduced in the form of w mi = b x,mi b y,mi , with b x,mi = 1 − x,mi ∕ x and b y,mi = 1 − y,mi ∕ y . Finally, the normal stress acting on each grid point is estimated from The values of A N can be compared with the approximated contact area of an individual atom A a,Ca ≈ 12.2 Å 2 [5] by defining a ratio (N) = A N ∕A a,Ca . It can be found that with N = 72 (for the current study), this ratio would be ≈ 1 . Various values of N ranging between 2 and 128 are used in this investigation, which are summarized in Table 1 along with their corresponding values of . One may argue that N > 72 results in < 1 , which is physically meaningless; accepting this argument, the possible use of gridding with < 1 is discussed in Sect. 3.2.

The Method for Identification of Local Stress Maxima and the Line of Maximum Gradient of Stress Distribution
The filtered stress distributions obtained from binning with N = 128 are post-processed with in-house codes written in MATLAB 9.3 particularly for this purpose. Firstly, in order to increase the accuracy of the procedure, the data are mapped onto a 1024 × 1024 rectangular gridded plane. Then, the local stress maxima are identified using function "imregionalmax" available via the Image Processing Toolbox 10.1. While the point of maximum stress can be easily obtained for the smooth contact, maxima in the stress distribution of the rough one have to be inspected as it is possible that various local maxima exist. Afterwards, the gradient of the stress is calculated, and the line of maximum gradient around each local stress is obtained through the following procedure. The Cartesian coordinate is transformed into a polar one, by setting its pole at the position of the selected maximum. Then, the angular coordinate is divided into 36 equivalent sectors, so that a single point of the maximum gradient line can be identified within each one, resulting in a total of 36 data points around the selected maximum. To do so, the stress values within each sector are sorted as a function of their radius, and the local maximum with the smallest radial coordinate is selected as one point on the line of maximum gradient. After processing all 36 sectors, the coordinate system is transformed back into a Cartesian one. Once the procedure is applied to all of the identified local maxima, the contact area is estimated for the traced maximum gradient lines using function "boundary." Figure 2 shows the contact's normal force F z as a function of indentation depth for both contacts. In order to convert the displacement to indentation depth, ind = 0 is selected from the position of the indenter at one step prior to the initiation of interaction between the bodies, i.e., where F z ≠ 0.

Results and Discussion
In both contacts, the normal force increases to a maximum and is followed by a sudden drop, indicating the initiation of irreversible plastic deformation within the deformable body. By fitting a power law equation of the form F z = c e ind to the force-indentation depth curves, the values of c and e are found to be 3.5 and 1.76 for the smooth tip, and 1.61 and 2.14 for the rough one, respectively. Comparing the normal forces for the two contacts prior to the force drops, it is obvious that, at any given indentation depth, the contact force for the smooth indenter is higher than for the rough one. This is due to the higher number of interacting atoms for the case of a smooth indenter. Furthermore, it can be noted that the contact force for the rough indenter features a sudden change in its slope that first occurs at ind = 4.4 Å , which corresponds to F z ≈ 37.4 nN . The reason for this change is the lateral rearrangement of the interacting atoms of the deformable body. Figure 3 shows the displacement vectors of the substrate's surface atoms calculated from the atoms' positions at two indentation depths of 4.3 Å and 4.4 Å .
Lateral rearrangement of the substrate's atoms

Hertz Solution
The simulated contact deviates from the Hertz contact theory in different ways: (1) the surfaces are not continuous, (2) although the contact is investigated within the elastic deformation regime, elasticity is not isotropic, (3) the contact is not frictionless due to the inevitable roughness at the atomic scale, and (4) a finite-range repulsion, as opposed to a hardwall interaction, describes the interaction between the bodies. Nevertheless, it can be argued that the contacts of the atomistic simulations and continuum theories could be compared, if their stress distributions would go through the same filtering procedure [11]. The Hertz stress distribution has a form of (r) = 0 (1 − (r∕r c ) 2 ) 0.5 , where r c is the contact radius, and 0 is the maximum compressive stress. Correspondingly, the contact radius and maximum contact stress are respectively. Figure 4 shows the maximum stress of the Hertz theory for the investigated contacts. Moreover, following Hertz contact theory, the elastic contact modulus E * can be found as follows: The contact elastic modulus can be estimated via E * = E∕(1 − 2 ) = 28.57 GPa , with E = 26 GPa and = 0.30 for the simulated material [12]; however, it has already been discussed that, assuming the applicability of Hertz theory at the atomic scale, E * is not a constant, but varies with the indenter's roughness, radius of curvature, indentation depth, and the governing interacting energy between the tip and the deformable body [8]. Figure 5 shows the calculated E * from Eq. 3 for both the investigated contacts. As is expected, the results show that the estimated values are not constant: both start at a low value, and increase with F z , albeit with different slopes. The E * of the smooth indenter increases for the whole range of the elastic contact, The elastic contact modulus E * as a function of F z , estimated based on Hertz contact theory. The reference bulk value of E * = 28.57 GPa is found from the governing potential for the simulated material while the rough indenter's remains close to the bulk value. Appendix 1 shows that the samples' thickness values were large enough to not affect the estimated values of E * .

Filtering the Stress Distributions
In order to compare the stress distributions of the investigated contacts calculated for different numbers of bins N (further discussed in Sect. 2.4), i.e., (N) , the distributions are smoothened with a low-pass filter. Herein, following the suggestion in [11], we use a Gaussian filter of the form , with x m and y m being the coordinates of the center of the xy plane, which are set to zero in this investigation, is the width of the filter, and A = 1∕2 2 ensures that ∫ G( , r)dr = 1 . Hereafter, we denote this filter as G( ) . The filtered distributions (N, ) are calculated via where F and F −1 denote the Fourier and its inverse transform, respectively. In this investigation, we define as a function of the nearest neighbor distance, i.e., we employ (p) = pd NN , where p is a dimensionless prefactor; we apply different values of p between 1 and 6 with an interval of 0.5. The nearest neighbor distance can be obtained as d NN = a 0 √ 2∕2 ≈ 3.9516 Å for the simulated material. It should be noted that the force-displacement curves of the two indentation processes have no intersection in the elastic deformation regime (see Fig. 2), meaning that there is no state in both processes where their corresponding Hertz's stress distributions can be directly compared with each other; therefore, we present the results of this section for the state of the contacts prior to the first force drop at F z ≈ 55 nN , to illustrate how the filtering method reduces the deviations between the results of the atomistic simulation and the Hertz solution. Figure 6 shows the filtered stress distributions for (N ≤ 32) and p = 1 for both contacts. The results show that the filtered distributions of both contacts are quite similar for these small values of N.
Choosing a value for N is rather arbitrary, although it seems to be natural to define N by enforcing = 1 , i.e., the area of each bin is equal to the projected area of an individual atom; however, the filtered stress distribution would not change for (N) < 1 ; see Appendix 2. Therefore, it is wise to select an arbitrarily large value of N, e.g., 128 as in this investigation.
In order to compare the theoretical stress distributions with the ones obtained from the atomistic simulations, they should all be processed through the same filter.
compares the filtered stress distributions of both contacts and their corresponding Hertzian ones for different values of . The results show that a filter with a resolution of = 5d NN ≈ 2 nm smears out the stress distributions such that the contacts become indistinguishable.

Interacting Area
Prior to the force drop, the projections of the contact stress distributions (indicating the interacting area) for N ≤ 32 are found to be almost identical for both contacts. The results are reported in Appendix 3 and show that at N = 32 the estimated interacting area becomes localized and circular. Once a value of N = 72 is employed, at which (N) = 1 , the differences between the interacting areas of the two systems can be easily noticed, as shown in Fig. 8. The interacting bins of the smooth indenter are neighbors, and follow the crystallography pattern of the contacting bodies, while the rough tip results in randomly distributed interacting bins. Figure 9 shows a comparison between the interacting area as a function of F z for selected values of N. Moreover, it reports the results of two other methods: prediction of Hertz contact mechanics ( A Hertz = R ind ) and estimation of the real interacting area. For the latter, the interacting atoms at each indentation depth are counted and denoted T int , and the real interacting area is estimated via A int,real = T int A a,Ca [9]. The results show that, with increasing N and consequently decreasing bins' sizes, the estimated interacting areas become smaller. Moreover, A N=72 is found to be close to A int,real at low values of F z . The deviations between these two values could be related to the differences between the real shape of the contact and its projection [9]: for the smooth indenter, the deviation could be simply due to the curved shape of the contact, while severe deviations could be expected for the rough tip. Therefore, as soon as the real shape of the interacting area starts to deviate from its projection on the 2D binned space, these methods result in different estimations.
Comparing the results of A N=72 with A Hertz shows that the interacting area for the smooth indenter is close to A Hertz at the beginning of the indentation process, i.e., F z ≤ 17 nN ( ind ≤ 2.5 Å ), and the Hertz prediction gradually underestimates the contact area. For the rough contact, however, Hertz theory overestimates the interacting area A N=72 for the entire process.

How Does Filtering Affect the Estimated Contact Area?
Once the stress distributions are filtered, no sharp border can be simply identified as the contact line; however, in order to use the filtered distributions for estimation of the contact area, it is necessary to define contact. Considering that the filtered distributions are differentiable, let us assume that the contact line is along their maximum gradient. Figure 10 shows the contact lines for both contacts prior to the force drop at F z ≈ 55 nN . The results show that the estimated contact lines for most of the contacts are more or less circular. Note that, for some cases (insets e and f), a number of local stress maxima were identified; however, it can be noticed that not all of the local maxima could be identified using the employed method (Sect. 2.5). This post-processing is performed on the entire indentation process of both atomistic simulations and their corresponding Hertz solutions. Then, the projected contact areas ( A pc ) are estimated as the enclosed area inside the identified contact lines. Figure 11 shows the estimated values of A pc for the atomistic simulations at different contact forces as functions of p. The first noticeable finding of this figure is some large initial  contact areas obtained from the filtered stress distributions. This is due to our method of defining contact area using the filtered data: the estimated contact area of an individual atom becomes increasingly larger as a wider filter is applied for smoothing the stress distribution. Moreover, the results show two general trends: (1) for low contact forces, the estimated A pc increases with p, and (2) for high contact forces, the estimated A pc decreases to a minimum followed by an increase, by increasing p. The reason for the decrease in A pc is that the small values of p preserve the non-uniformity of the stress distribution, which is a consequence of the inevitable surface roughness at the atomic scale. Increasing p smears out the effects of the surface roughness; thus, the identified contact area becomes more circular, and results in a smaller A pc . Further increasing p, however, makes the stress distribution spread out more, resulting in larger A pc . In order to compare the estimated values of A pc with the Hertz solution, Fig. 12 shows A pc as functions of contact force. The results show that, for either indenter, a value of p can be found at any given contact force resulting in a projected   Fig. 13. It can be seen that a unique value of p cannot be found to accurately apply to the whole range of the simulated process. See Appendix 4 for further discussions on applying a single value of p. Figure 14 shows the values of A pc estimated using two sets of p: (1) those at which the corresponding estimated contact area has its minimum value, i.e., the black lines in Fig. 11, and (2) those reported in Fig. 13. The first finding of this figure is that the filters diminish the deviations between the estimated projected contact areas obtained from the atomistic simulation A pc,sim and Hertz theory A pc,H . Moreover, the results show that the estimated contact behavior depends on the applied filter. Using the first set of p, the values of A pc for the smooth tip obtained from both the atomistic simulation Smooth Rough Fig. 13 The identified values of p at which the estimated A pc were closest to the Hertz solution and Hertz theory are lower than the original estimations for the entire indentation process. For the rough tip, however, A pc,sim is lower than A pc,H before the initiation of the lateral rearrangement of the substrate's atoms ( F z ≈ 37.4 nN ) and becomes higher than A pc,H at larger indentation forces. Applying the second set of p, the values of A pc,sim are either identical to or higher than A pc,H for the entire process. Moreover, the deviations between A pc,sim and A pc,H become significant only for the rough contact for indentation forces larger than F z ≈ 37.4 nN . Furthermore, as is expected, the values of A pc are in good agreement with the Hertz solutions; however, in comparison with A N=72 , the filtering method results in under-and overestimations of A pc for the smooth and rough indenters, respectively. The results reveal an important possible feature of the filtering method. Differentiating between the interacting and contacting atoms leads to defining a corresponding contact area A c smaller than the interacting one A int [5,8]; the validity of this approach has been shown for studying atomistic rough surfaces [6]. For the current studied rough contact, however, the behavior is partly (p set 1) or entirely (p set 2) reversed, i.e., A c is found to be larger than A int . This physically unjustifiable artefact is a consequence of the applied filter, and, perhaps, our definition of a contact line.
Furthermore, it should be noted that, while the first set of p can be identified by varying the width of the filter, the preparation of the second set requires a priori knowledge of a reference, e.g., the Hertz solution as in the current study, which can be obtained using indentation test measurements. Note that, although A N=72 can be used as the reference as well, this solution is not directly accessible from an experimental test. This means that, in practice, this solution is not available for calibrating the filter's width. The presented results of the Hertz model in Fig. 14, both the original and the filtered ones, are based on the Hertz stress field established from r c and 0 estimated via Eqs. 1 and 2, respectively; therefore, they require the values of R and ind in advance. Considering the goal of this work, i.e., investigating the applicability of the filtering method, the results which were dependent on some already known data serve solely for validating purposes.

Accounting for Roughness
Another way of studying the rough contact is to consider the area of the local contact patches Ã pc (see Fig. 15), and estimating A pc = ∑Ã pc . Figure  be noted that, in the calculation of ∑Ã pc , the intersected areas are treated by merging them into one contact patch, if the distances between the corresponding local maxima are found to be smaller than closest approach between their identified maximum gradient lines.
Considering the surface roughness of the contact, the contact area is estimated using the following analytical expression for A pc of a Hertzian tip with small-scale roughness [20]: with where g is the root-mean-square gradient of the surface roughness estimated to be ∼ 0.7 for the studied rough contact after removing the global curvature of the indenter [21], and ≈ 1.6 [22]. Figure 16 shows the solutions of Eq. 6 (solid line), resulting in values of A pc close to A N=72 for the entire indentation process.

Conclusions
In this paper, the idea of applying Gaussian filters on atomistic contact stress fields for smearing out the effects of surface roughness is investigated. Through a binning procedure, the interacting areas for each contact are estimated; comparing the results with the corresponding solutions of Hertz theory shows that the model under-and overestimates the interacting areas for the smooth and rough contacts, respectively. Furthermore, the stress distributions of the atomistic simulations and analytical solutions are locally averaged by applying Gaussian filters with different widths ( ), and the contact line is defined as the line of maximum gradient of the stress distribution surrounding the contact patch.
The results of this study show that the filtering method improves the correlations between A pc of the atomistic simulation and Hertz theory within the range of elastic deformation for both contacts, albeit with different accuracies depending on the worked filter (Figs. 14, 20). Although the choice of is rather arbitrary (Sect. 3.2), the effect of varying its value on the estimated A pc shows that the applied filter should have a varying width; however, defining the correct range of is problematic. Our suggestion for identifying via finding the minimum of A pc as a function of the filter's width does not result in accurate estimations of A pc for both contacts and for the entire indentation process. On the other hand, identifying a set of p using a reference, e.g., the Hertz model, requires the assumption of applicability of that reference at the atomic scale. Our presented results show that either of these approaches for identifying the filter's width leads to deviations between the estimated values of A pc from the filtered stress distributions and those obtained directly from the atomistic simulations, without the possibility of determining a validity range. Therefore, while the filtering process reduces the deviations between the behavior of a non-adhesive atomistic contact and the corresponding Hertz solution, it does not validate the assumption of the applicability of the Hertz model for such a contact, especially for the rough case.
Furthermore, by taking surface roughness into consideration, we define A pc to be the sum of the contact area of local patches, i.e., ∑Ã pc . By comparing the results with the benchmark A N=72 , a filter with = 1.5d NN is identified for its accurate estimations of the projected contact area, prior to the lateral rearrangements of the surface atoms. Finally, our results with a model for rough Hertzian tips (Eq. 6) suggest that, instead of filtering the contact's pressure distribution, selecting an appropriate model or theory for the relevant physics may yield more accurate solutions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Eq. 6 [20] ( ) Fig. 16 The values of the interacting area A N=72 (open square) and A pc estimated via ∑Ã pc using different values of p (dashed lines) and Eq. 6 (solid line) show that enlarging the filter's width has two effects: (1) the values of A c,Sim and A c,H are close for larger indentation depths, and (2) the deviations between the estimated contact area from the two methods become smaller. The same general trend can be seen for the rough contact ( Figure 20b) for p = 2.5 and 5, although the deviations are more pronounced than the corresponding trends of the smooth contact. Moreover, the filter with p = 1 results in A c,Sim being smaller than the corresponding A c,H for the whole range of the indentation depth. The estimated contact area for the rough contact with p = 1 is smaller than both the estimated interacting area and the corresponding Hertz solution, showing no advantage over the calculated interacting areas for comparison with the continuum theory. Note that the results labeled as "original" are the data of Hertz and A N=72 before filtering (see Fig. 9) (Color figure online)