How Thermal Fluctuations Affect Hard-Wall Repulsion and Thereby Hertzian Contact Mechanics

Contact problems as they occur in tribology and colloid science are often solved with the assumption of hard-wall and hard-disk repulsion between locally smooth surfaces. This approximation is certainly meaningful at sufficiently coarse scales. However, at small scales, thermal fluctuations can become relevant. In this study, we address the question how they render non-overlap constraints into finite-range repulsion. To this end, we derive a closed-form analytical expression for the potential of mean force between a hard wall and a thermally fluctuating, linearly elastic counterface. Theoretical results are validated with numerical simulations based on the Green's function molecular dynamics technique, which is generalized to include thermal noise while allowing for hard-wall interactions. Applications consist of the validation of our method for flat surfaces and the generalization of the Hertzian contact to finite temperature. In both cases, similar force-distance relationships are produced with effective potentials as with fully thermostatted simulations. Analytical expressions are identified that allow the thermal corrections to the Hertzian load-displacement relation to be accurately estimated. While these corrections are not necessarily small, they turn out surprisingly insensitive to the applied load.


INTRODUCTION
One of several drawbacks when applying continuum theory to small-scale contact problems, as they occur, addition, large standard deviations of experimentally measured depinning forces of atomic-force microscope 23 tips have been observed, which were accompanied by unexpectedly large reductions of the depinning force 24 with increasing temperature [8]. It is possible that thermal surface fluctuations, which were not included 25 in the modeling of temperature effects on tip depinning, are responsible for a significant reduction of 26 effective surface energy and thereby for a reduction of the depinning force. In fact, it has been shown that 27 thermal fluctuations limit the adhesive strength of compliant solids [9]. Finally, in the context of colloid 28 science, it may well be that thermal corrections have a non-negligible effect on the surprisingly complex 29 vector, and q its magnitude. 66ũ (q) = 1 A d 2 r e −iq·r u(r) (2) denotes the Fourier transform of u(r). The short-hand notation u 0 =ũ(q = 0) will be used for the 67 center-of-mass coordinate. 68 For flat indenters, only u 0 will be used to denote the mean distance, or gap, between indenter and the 69 solid surface. Here, we define the displacement d as a function of temperature and load according to where d 0 denotes the displacement for an ideal, athermal Hertzian indenter at a given load. In the current 78 work, we compute d T through the following approximation where r X is the point that is the most distant from the center of the Hertzian indenter. We found that the 80 first three to four digits are accurate in this estimate if the athermal Hertzian contact radius is less than 81 one quarter of the simulation cell's linear dimension. This is because the (true) surface displacement fields 82 converge quite quickly to their asymptotic 1/r form outside the (original) contact radius in the case of 83 short-ranged potentials and because the finite-size corrections to the true surface displacements are not 84 very sensitive to temperature. 85 The interaction with a counterface is modeled within the Derjaguin approximation [18] so that the surface 86 energy density depends only on the local interfacial separation, or, gap, g(r) = u(r) − h(r), between the 87 surfaces, i.e., the interaction potential is obtained via an integration over the surface energy density γ(g) via In the full microscopic treatment, hard-wall repulsion is assumed, i.e., 89 γ(g) = ∞ if g < 0 0 else .
Finally, the probability of a certain configuration to occur is taken to be proportional to the Boltzmann 90 factor, i.e., where β = 1/k B T is the inverse thermal energy. 92 One central "observable" in this work is the distance dependence of the mean force per atom, f (u 0 ), for 93 flat surfaces and finite temperatures. One might want to interpret this function as a cohesive-zone model, or, 94 in the given context better as a repulsive-zone model. Because of the so-called equivalence of ensembles, 95 which is valid for sufficiently large, systems, it does not matter if the separation is fixed and the force 96 measured, or, vice versa.

97
Note that we will go back and forth between continuous and discrete descriptions of displacement fields. 98 For the discrete description, the elastic solid is partitioned into atoms, which are arranged on a square 99 lattice with the lattice constant ∆a. This was done for reasons of simplicity, even if other discretizations 100 are possible, e.g., into a triangular lattice [15]. Transitions between discrete and continuous representations 101 in real space can be achieved with the substitutions while transitions between summations and integrals in the wavevector domain can be achieved with To simplify the analytical evaluation of integrals, the square Brillouin zone (BZ) of the surface will be 104 approximated with a circular domain. In this case, the upper cutoff for q is chosen to be q max = √ 4π/∆a whereF (q, t) is the Fourier transform of all external forces acting on the surface atoms. The terms m q and 114 η q represent inertia and damping coefficients of different surface modes, which may depend on the wave 115 vector. For isotropic systems, these terms only depend on q but not on the direction of q.

116
The effect of thermal fluctuations can be cast as random forces, which have to satisfy the fluctuation-117 dissipation theorem (FDT) [19]. In the given formalism, random forces must have a zero mean, while their 118 second moments must satisfy, Good numbers for the exponent n and the dimensionless hard-wall stiffness κ o need to be chosen. In 156 order for the effective hard-wall potential to have a minimal effect on ∆t, the (non-negative) exponent n explained at length in any better text book on molecular dynamics [11,12]. While these arguments can be somewhat academic when the discontinuities are small, we are going to send κ o to large numbers resulting in significant force discontinuities. Thus, n must be chosen greater equal two. This appears to make n = 2 161 the optimal choice.

162
The next question to be answered is: Given a time step ∆t and an exponent of n = 2, what is a good value 163 for κ o ? Here, it is useful to keep in mind that we do not need very accurate dynamics in the "forbidden" 164 overlap zone. The main purpose of the stiff harmonic potential is to eliminate overlap as quickly as possible, 165 i.e., to effectively realize a collision of the particles with the position of the (original) hard wall. However, 166 the stiffness should remain (well) below a critical value above which energy conservation is violated in 167 the absence of a thermostat even when a symplectic integrator, such as the Verlet algorithm, is used. For

168
Verlet, the critical time step for a harmonic oscillator is ∆t c = T /π, where T is the oscillator period, i.e., 169 for ∆t < ∆t c , the trajectory may be inaccurate, but the energy is conserved (except for round-off errors).

170
This can be achieved by setting the overlap stiffness to where k s = ∆u 2 /(k B T ), while m is the inertia of the considered degree of freedom. ν o is a numerical 172 factor, which must be chosen less than unity. At and above the critical value of ν o = 1, energy conservation 173 would be no longer obeyed in the absence of a thermostat. At the same time, dynamics but also static 174 distribution functions are very inaccurate, even if a thermostat prevents the system from blowing up.

175
The optimum value for k o certainly depends on the specific investigated problem. However, the analysis 176 of simple models can provide useful preliminary estimates. This will be done in Sect. 2.3.3.

Numerical case studies 197
To explore the relative merit of the two proposed hard-wall methods, we investigate the following 198 single-particle problem: an originally free, harmonic oscillator with a (thermal) variance of ∆u 2 . This 199 harmonic oscillator is then constrained to have no negative deflections from its mechanical equilibrium 200 site. The analytical solution to this problem stating the force F needed to realize a given constraint is 201 contained in the mean-field approximation to the full elastic problem, which is presented in Sect. 3.2.2.

202
The given constraint of the spring sitting exactly on the hard wall corresponds to a value, where u 0 203 crosses over from its short-range to its long-range asymptotic behavior. Therefore, we see this case as 204 being representative for both scaling regimes.

205
In essence, the problem we investigate corresponds to the choice where k B T , k, and m are used to define  At a given value of ∆t, the approximate collision rules clearly outperform the approximate hard-wall 213 interactions. However, u 0 has leading-order corrections of order ∆t in both approaches. With the choice 214 ν o = 0.1, the asymptotic result for the parabolic, effective hard-wall potential has an accuracy of better run at two different values of ∆t, say e.g., at ∆t = 0.25 and ∆t = 0.15 in order to perform a meaningful 217 ∆t → 0 extrapolation. In a full contact-mechanics simulation, the number of required Fourier transforms 218 doubles when using the approximate collision rules, which in turn leads to increased stochastic errors given 219 a fixed computing time contingent. For this reason, but also because approximate collision rules require 220 significantly more coding-in particular when averaging wall-surface forces from collisions when using 221 wavevector dependent inertia-we decided to use the harmonic, effective hard-wall potential for the full 222 contact-mechanics simulations.

THEORY
The main purpose of this section is to identify an analytical expression for the thermal expectation value Minor errors in the treatment presented below appear in numerical coefficients that result, for example, 229 by having approximated the Brillouin zone of a square with a sphere, or, by having replaced a discrete set 230 of wave vectors (finite system) with a continuous set (infinitely large system). However, these and related 231 approximations are controlled, because errors resulting from them can be estimated and they could even be 232 corrected systematically. Since the free surface is the reference state, we start with its discussion. An important quantity, in 235 particular in a mean-field approach, is the variance of atomic displacements due to thermal noise. For a 236 fixed center-of-mass coordinate, it is defined as the following thermal expectation value: It can be evaluated in its wavevector representation in a straightforward manner. Specifically, where we made use of equipartition for harmonic modes, see also Eq. (29).

239
Of course, up to the prefactor of 2/ √ π ≈ 1.1284, which is very close to unity, Eq. (19) follows directly 240 from dimensional analysis. However, in a quantitative theory, we wish to know and perhaps to understand 241 its precise value. A numerical summation over a square BZ assuming a square real-space domain with N atoms reveals that ∆u 2 can be described by In 255 reality, i.e., in less than infinite dimensions, there is always a correlation of thermal height fluctuations.

256
To deduce an estimate for the distance over which height fluctuations are correlated, we calculate the 257 thermal displacement autocorrelation function (ACF) C uu (r). It can be defined and evaluated to obey: where J 0 (x) is the Bessel function of the first kind and 1 F 2 (...) is a generalized hypergeometric function.  A quite reasonable approximation or rather generalization of C uu (r) to a continuous function can be made 266 by constructing the simplest expression with the correct asymptotic behaviors summarized in Eq. (26): As can be seen in Fig. 2, this asymptotic approximation is quite reasonable already at a nearest-neighbor 268 spacing of r = ∆a and has errors of less than 5% (in the limit of large N ) for larger values of r. While The asymptotic ACF has decayed to approximately 30% of its maximum value at the nearest-neighbor 273 distance. This means that the displacements of adjacent lattice sites are essentially uncorrelated.

274
The last property of interest used in the subsequent treatment is the partition function of a free surface 275 (fs): using path-integral techniques [24]) and to assign m q such that it reproduces the correct zero-point variance 285 of the mode in question.

286
In the mean-field (Einstein solid) approximation, the partition function simplifies to with ∆u having been introduced in Eq. (19) and λ mf being a mean-field de Broglie wavelength.  The arguably simplest analytical approach to the contact problem is an adaptation of the so-called Einstein 296 solid, which was already alluded to in Sect. 3.1, to surface atoms. We first do it such that a degree of 297 freedom is a hybrid of an atom in real space and a delocalized, ideal sine wave. Specifically, we first assume 298 that elastic energy of an individual atom reads In order to maintain a zero expectation value of u, it is furthermore assumed that the interaction energy 300 with a counterface placed at a distance u 0 from the atom's mean position is given by This means, an oscillation of an atom entails an undulation. With this assumption, u 0 automatically 302 corresponds to the atom's mean position.

303
The excess free energy per particle ∆F/N for a fixed center-of-mass position satisfies where the term "excess" refers to the change of the free energy relative to that of a free surface. For 305 hard-wall interactions, the integral in Eq. (33) can be evaluated to be Hence, For reasons of completeness, the force predicted from this first mean-field approximation is stated as: In the limit of u 0 → 0, repulsion diverges proportionally with 1/u 0 , while it decays slightly quicker than where f needs to be chosen such that u = u 0 so that the lattice position of the particle u eq is situated 316 at u eq = u 0 + βf ∆u 2 . At u eq , there is no elastic restoring force in the spring. The requirement u = u 0 317 automatically leads to the following self-consistent equation for f : This line of attack leads to similar results for the f (u 0 ) at small u 0 as the first mean-field approach.

319
However, for large u 0 the predicted force turns out half that of the first mean-field approximation. In fact, 320 the second mean-field theory turns out to be a quite reasonable approximation to the numerical data for any value of ∆u, see the results and discussion presented in Sect. 4.

Probabilistic approach 323
The exact expression for the excess free energy of an elastic body in front of a hard wall can be defined 324 by a path integral, where D[u(r)] denotes an integral over all possible displacement realizations and In the case of hard-wall repulsion, the r.h.s. of Eq. (40) is easy to interpret: It represents the relative number 327 of configurations that are produced with the thermal equilibrium distribution of a free surface (fs), whose 328 maximum displacement is less than u 0 , i.e., where µ gev is the mode of the Gumbel distribution, i.e., the most likely value for u max to occur, and β gev a 338 parameter determining the shape of the distribution. For a normal Gaussian distribution Φ G (u/∆u), they 339 are given by in the limit of large N . Here erf −1 (...) stands for the inverse function of the error function [25].

341
In fact, Fig. 3 shows that the distribution of u max as produced with GFMD and by taking the maximum Pr(u max < u 0 ) = erf and therefore This result turns out to apply to large separations, that is, to u 0 /∆u 1. The functional form of F(u) is 352 identical to the one obtained in the first mean-field variant, except for the prefactor, which is reduced by a 353 little more than a factor of two.
where the contact modulus E * and the contact radius R c were effectively used to define the units of pressure where F pa ≡ F/N denotes the hard-wall, free-energy normalized to the atom. The approximation in where the last approximation is only valid at small temperatures. Taylor expanding this last expression 376 leads to

High-temperature approximation 379
At very large temperatures, d T is in excess of d 0 so that deformations of the elastic solids are very small.

380
In a first-order perturbative approach, it then makes sense to assume the displacement field to be a constant, 381 i.e., to be d T . In that approximation, individual forces can be simply summed up with a mean gap of d T + r 2 n /(2R c ). Recasting the sum as an integral yields Eq. (58) can be solved for d T with the help of the Lambert W function W (x) ≈ ln x − ln ln x for x 1:

386
In this section, we investigate to what extent the three approaches introduced in Sect. 3.1 reproduce mean-field approach appears to be asymptotically exact for small u 0 , while the approach based on the law 392 of large numbers seems to be asymptotically exact for large u 0 . In between these two regimes, there is a 393 smooth transition between them. This transition is reflected quite well by the second mean-field approach.

394
Unfortunately, we did not identify a closed-form analytical expression for it, which would nevertheless be 395 nice to have when implementing a potential of mean force into a simulation. However, as is demonstrated 396 in Fig. 4, simple switching functions introduced next allow one to approximate numerical data reasonably 397 well.

398
Since both force-distance asymptotic dependencies have the same functional form and since the transition 399 between them is quite continuous, it is relatively easy to come up with switching functions allowing 400 the numerically determined free energy to be approximated reasonably well. Defining F mf1 through the 401 free-energy expression in Eq. (35), this is done via with the weighting functions

413
The solution of the continuous displacement field has no dimensionless number if the contact radius a c is 414 taken to be the unit of length. However, a c /∆a starts to matter as soon as it is no longer large compared to 415 unity. Since discreteness effects are a different topic discussed elsewhere [26], a c /∆a is chosen sufficiently 416 large so that the discrete problem reflects the continuous Hertz contact reasonably well.

417
To test the applicability of the thermal repulsive-zone model in the realm of Hertzian contact mechanics, 418 the following parameters were chosen as useful defaults after some trial and error: R c = 256 ∆a and a 419 normal load of L = 131 E * ∆a 2 leading to a c ≈ 30 ∆a within regular Hertzian contact mechanics. In the 420 athermal Hertzian contact, the mean contact pressure turns for these parameters is p ≈ 0.049 E * . Results for the stress profile at a temperature of k B T = 0.2 E * ∆a 3 are shown in Fig. 5.  Figure 5. Left: Interfacial stress σ as a function of distance r from the symmetry axis in a Hertzian contact geometry. The (blue) circles reflect zero temperature data from the hard-wall overlap potential. The full (blue) line represents the analytical solution to the Hertz problem. The (red) open squares show finitetemperature data from full simulations, while the (red) dotted line shows zero-temperature simulations, in which, however, the effective potential was constructed to reflect thermal vibrations at the given temperature. The arrow marks the point of largest slope for the thermal indenter. Right: Displacement field u(r) as a function of distance r from the symmetry axis.
An interesting but perhaps also obvious outcome of the data presented in Fig. 5 is that there is no abrupt 423 transition from finite to zero contact stress, once thermal fluctuations are finite. This observations is of 424 relevance when discussing the concept of "true contact area". Since collisions in a hard-wall potential are 425 instantaneous, the probability of observing two (finite) surfaces to be in contact has a statistical measure of 426 zero, so that the instantaneous contact area could be argued to be (almost) always zero. Contact exists only 427 in the isolated moments of time at which collisions take place. However, during these isolated moments of 428 time, the forces between surfaces is infinitely large such that time averaging yields a distribution which 429 resembles the well-known Hertzian stress profile; the smaller the temperature the closer the stress profiles 430 between original and finite-temperature stress profiles.

431
The question of how to meaningfully define (repulsive) contact area when repulsion has a finite range and 432 adhesion is neglected arises naturally. In a recent paper [26], it was proposed to define the contact line (or 433 edge) to be located, where the gradient of the normal stress has a maximum slope. In the current example, 434 this leads to a reduction of the contact radius of order 1%, which is significantly less than the reduction of 435 approximately 30% of the normal displacement in the given case study.

436
In contrast to contact radii, force and displacement can be defined unambiguously. Thermal noise will and The master curve shown in Fig. 6 reveals asymptotic regimes at low and at high temperatures, respectively.

447
They can be approximated with power laws. However corrections logarithmic in temperature need to 448 be made at low temperature to obtain quantitative agreement over broad temperature ranges. We find 449 numerically that Inserting the low-temperature approximation of the master curve into Eq. (64) and reshuffling terms In other words, the elastomer surface is effectively shifted by a little less than 1.5 times the thermal standard

SUMMARY
In this work, we analyzed the effect that thermal fluctuations can have on contact mechanics in the case of hard-wall interactions. To this end, we first demonstrated that thermal surface fluctuations are dominated 505 by short wavelengths undulations. They smear out the originally infinitesimally short-range repulsion to a 506 finite range of ∆u ≈ k B T /(E * ∆a). The functional form of the repulsive force was derived analytically 507 and shown to diverge inversely proportionally with the interfacial separation u 0 at small u 0 but to decay 508 slightly more quickly than exponentially in −u 2 0 at separations clearly exceeding ∆u.

509
To come to these results, the Green's function molecular dynamics (GFMD) technique was generalized 510 to include thermal noise. Particular emphasis was placed on the question how to handle (the original)