Modeling the Contact Mechanics of Hydrogels

A computationally lean model for the coarse-grained description of contact mechanics 1 of hydrogels is proposed and characterized. It consists of a simple bead-spring model for the 2 interaction within a chain, potentials describing the interaction between monomers and mold or 3 confining walls, and a coarse-grained potential reflecting the solvent-mediated effective repulsion 4 between non-bonded monomers. Moreover, crosslinking only takes place after the polymers have 5 equilibrated in their mold. As such, the model is able to reflect the density, solvent quality, and the 6 mold hydrophobicity that existed during the crosslinking of the polymers. Finally, the such produced 7 hydrogels are exposed to sinusoidal indenters. The simulations reveal a wavevector-dependent 8 effective moduls E∗(q) with the following properties: (i) stiffening under mechanical pressure, 9 and a sensitivity of E∗(q) on (ii) the degree of crosslinking at large wavelengths, (iii) the solvent 10 quality, and (iv) the hydrophobicity of the mold in which the polymers were crosslinked. Finally, 11 the simulations provide evidence that the elastic heterogeneity inherent to hydrogels can suffice to 12 pin a compressed hydrogel to a microscopically frictionless wall that is ondulated at a mesoscopic 13 length scale. Although model and simulations of this feasibility study are only two-dimensional, its 14 generalization to three dimensions can be achieved in a straightforward fashion. 15

In the following, different aspects of the models are described in separate sections. Some of them  A simple model or representation of polymers is the freely jointed chain, in which the polymer is partitioned into P K Kuhn segments of fixed length b, which have no angle correlation within the chain. In the Kuhn model, the mean-square end-to-end radius is given by It applies to so-called θ-solvents, in which R 2 EE scales linearly with the real degree of polymerization 74 P real for sufficiently long polymers. Note that a polymer is its own θ-solvent. How P K and b are 75 conventionally chosen is discussed further below.

76
The distribution function Pr(R EE ) ressembles that of a Gaussian -except for very large values of R EE . Thus, the coarsest possible description of a polymer, in which (small) mechanical forces only apply to chain ends -as it happens when crosslinkers attach to them -is achieved by modeling the entire polymer with a single thermalized spring having a zero equilibrium length and a stiffness of where k B T is the thermal energy and D the spatial dimension. Thermal fluctuations need to be reflected, 77 because they encapsulate the effect of the polymer's internal degrees of freedom. Alternatively, the 78 chain could be partitioned into P G thermal springs, each of which having a stiffness of P G k EE to yield 79 another good approximation to the true R EE distribution function. This constitutes a Gaussian chain, 80 whose configurations can be interpreted as realizations of random paths satisfying -on average -the 81 diffusion equation. Partitioning the Gaussian chain into many small segments allows the (distribution 82 of) polymer configurations to be represented when external potentials are present, e.g., in the form of 83 an inpenetrable wall. 84 As an incidental remark, it is noted that the freely-jointed chain and the Gaussian chain show the 85 same linear response regarding their force-extention characteristics. Differences in the two models 86 only become apparent when the polymers are streched by a significant fraction of R EE . The Gaussian 87 chain does not change its stiffness, while the freely-jointed chain becomes infinitely stiff when the two 88 ends are pulled P K b apart. Of course, real force-extention curves lie in between those limits, see also 89 Sect. 2.2.3.

90
A Kuhn segment may contain several repeat units that cannot be approximated as being freely 91 jointed due to implicit or explicit bond angle potentials. If such potentials exist, the thermal expectation 92 value cos ϑ can differ from zero, where cos ϑ n = a n · a n+1 /a 2 0 . Here a n is the bond -or rather segment 93 -vector n, which is numerated from 1 to P r , while a 0 refers to the fixed segment length and P r to the 94 number of repeat units into which a Kuhn segment is discretized. Note that the product P r P K may but 95 does not have to coincide with P real . In addition, the term "segment" refers to a fraction of the polymer, 96 which may represent just one but, in most cases, is meant to encompass several (real) monomer units.
Assuming that the probability density for bond or segment angles satisfies Pr(ϑ) = Pr(−ϑ) and that thermal fluctuations of different bond angles are independent from each other, it follows that exp i n 1 +∆n ∑ n=n 1 ϑ n = cos ϑ ∆n (3) This in turn implies that the correlation of the chain's orientation decays according to cos ϑ ∆n = e −∆n a 0 /λ p (4) with the so-called persistance length being If the chain length is much less than λ p , it behaves like a rigid rod, while chains are floppy when their 98 lengths clearly exceed λ p .
(3) also allows the end-to-end radius in the more refined model to be determined to with c = cos ϑ . In the limit c P r → 1, the chain ressembles a stiff rod, in which case R EE = P r a, as can be deduced from L'Hôspital's rule. In the limit c P r 1, the chain ressembles a Gaussian chain satisfying R EE ∝ √ P r . The Kuhn length b is commonly chosen twice the persistence length, so that for long chains the number of Kuhn segments needs to be set to

101
When modeling the mechanical properties of a hydrogel it is desirable to dispose of a description 102 in which the representation of the polymer can be set arbitrarily somewhere between fine and 103 coarse-grained scales. We attempt to achieve this with the following energy expression where P is the number of segments into which the "polymer path" is discretized and r n−1,n = r n − r n−1 , jointed chain. In this limit, it would be most efficient to implement constraints rather than to use small 115 time steps. In the opposite limite, i.e., when making k 1 finite and setting a 0 and k 2 to zero, a Gaussian 116 chain is obtained, which might be useful for the modeling of a hydrogel in which extremely long 117 polymers are crosslinked at their ends. Chosing a positive value of k 2 introduces a bending stiffness to 118 the polymer, because this creates a force pushing a monomer towards the center-of-mass coordinates of 119 its two neighbors, whereby the chain straightens out. In the following, the generic properties resulting 120 from the energy expression in Eq. (8)  has a small effect on the root-means-square (rms) segment length a rms , which however is incorrect 124 unless k 2 is small compared to min(k 1 , k B T/a 2 0 ). The effect of k 2 on a rms is most easily quantified in a 125 one-dimensional model, because this constitutes a harmonic potential, which is amenable to analytical 126 treatments. Towards this end, we explore the continuum approximation that results from expressing 127 the displacements of atoms u n ≡ x n − x eq n,0 from the their equilibrium positions x eq n,0 = n a 0 , in which 128 case u n → u(x) and (u n − u n−1 ) → a 0 u (x). As such, we may write The continuum approximation to the model in one spatial dimension reads In a Fourier representation, this expression allows a wavenumber dependent stiffness to be defined as Due to Parseval's theorem, mean squares of fields can be either added up in real space or in the Fourier 130 representation. This is why we may write using equipartition withq = q a 0 and k eff = k 1 π 2 k 2 /(4 k 1 ) atan π 2 k 2 /(4 k 1 ) .
Assuming continuum dispersion relationsships up to q max = π/a 0 is certainly simplifying and 132 using the correct ones would be feasible, at least when integrating the pertinent expression for k eff 133 numerically. It would lead to a slightly reduced value of k eff , because the continuum approximation 134 overestimates the stiffness at the border of the Brillouin zone. However, it is felt that the current 135 treatment is sufficiently refined at this explorative stage and that having simple analytical expressions 136 is also advantageous.

137
An accurate analytical treatment of the model in two spatial dimensions is substantially more 138 elaborate than in one dimension, because the potential is no longer a harmonic function of the 139 displacements. The best we are currently able to do is to use k eff as an effective local bond stiffness. The 140 thermal average of the second moment of the instantaneous segment length ra 2 rms being nothing but r 2 averaged in two dimensions rather than in one -then becomes (β = 1/k B T, tildes indicate 142 that a variable is expressed in units of √ k B T/k eff ): 143ã 2 rms,2D = d 2 r r 2 e −β k eff (r−a 0 ) 2 /2 d 2 r e −βk eff (r−a 0 ) 2 /2 (14) = 2 +ã 2 0 e −ã 2 0 /2 + π 2ã 0 3 +ã 2 0 erf ã 0 so that he two asymptotic regimes satisfy ifã 0 1 (as for Gaussian chains) .
When k 2 a 2 0 is large compared to k B T, it is clear that the chain is locally quasi-one dimensional and that the one-dimensional formula for a 2 rms , see Eq. (12), is the appropriate choice, while in the opposite case, the two-dimensional treatment resulting in Eq. (15) is to be prefered. A switching function s is needed to embed both limits in one formula, e.g., through At this point, we are only in a position to argue heuristically. Assuming the switching function to be exponential in k 2 , the only functional form it may have from a dimensional analysis point of view is It was found that values of g less to but in the order of unity provide satisfactory cross-over functions.

144
This is demonstrated in Fig. 2, which considers one model in which k 1 = k 2 is chosen and another one 145 with k 1 being fixed at a value clearly exceeding k B T/a 2 0 .
(a)  The scaling dependence of a rms differs in the two graphs shown in Fig. 2, because the chain 147 becomes similar but not identical to an ideal, two-dimensional Gaussian chain in the limit k 1 = k 2 → 0. compared to a 0 -when k 1 is fixed such that it clearly exceeds k B T/a 2 0 . At large values of k 2 both models 150 behave asymptotically as effictively one-dimensional chains. The small errors in Fig. 2 Obtaining analytical estimates for R 2 EE requires approximations for cos ϑ to be known in addition 155 to a rms . Towards this end, the potential energy is rewritten in terms of the segment vectors a n as The model resulting from the approximation in Eq. (20) -in the limit k 1 a 2 0 /(k B T) → ∞ is also known withk 2 = β k 2 a 2 rms /4 and I n (...) being the modified Bessel function of the first kind.

159
Rather than reporting cos ϑ , Fig. 3 shows the dependence of R EE on k 2 for the cases k 1 = k 2 160 and k 1 a 2 0 /(k B T) 1. Analytical results are obtained by using our approximations for a rms , as 161 shown in Fig. 2, in Eq. (6). The resulting predictions turn out quite satisfactory except for the case 162 k 1 = k 2 k B T/a 2 0 . However, when coarse graining polymers in a systematic fashion, k 2 will 163 always disappear more quickly than k 1 , since the k 2 term scales with the inverse fourth power of the 164 wavelength, while the k 1 term only scales with the inverse square, as can be deduced, for example, 165 from Eq. (11). This would make our model approach the limit of an ideal Gaussian chain quite quickly 166 so that analytical expressions for R EE would be quasi exact. Remember that R EE exceeds a 0 times the 167 number of segments in that limit.

Monomer-wall interactions 179
In coarse-scale descriptions, the interactions between monomers and confining walls is commonly represented in terms of hard-wall potentials. While implementing such constraints is possible within molecular dynamics, we chose another route, namely, a finite but short-range repulsion between monomers and wall. We make the interaction as stiff as possible without having to reduce the time step by setting the curvature of the monomer-wall interaction potential similar in magnitude as the chain stiffness, specifically, where Θ(..) is the Heavyside step function, while r p is the (signed) distance of a monomer from the 180 wall. A negative sign of r p indicates an penetration into the wall. The rms-slope of the used wall 181 profiles is always sufficiently small to have a unique value for r p .
When mimicking walls attracting polymers more strongly than solvent, adhesion is added. To make it computationally tractable, it is chosen as Here γ is the surface energy and r max the range of interaction. between monomers belonging to different polymers are lumped together as intermolecular interactions.

188
The leading-order approach in a density-functional or self-consistent field theory (SCFT) to these 189 interactions is to make an extra monomer "pay" an energy density of χ · ρ(r) if positioned at a site

195
A similar approach, as those pursued in SCFT is adopted here, however, with two approximations to make the approach computationally efficient. First, the density of a mononmer is smeared out so that an estimate for the coarse-grained density is obtained where ∆a has the unit length and r n is the coordinate of monomer n. ∆a was set such that it was similar 196 to a 2 rms .

197
The second approximation is that the force acting on monomer n is not taken to result from the derivative of χ d 2 rρ 2 cg (r)/2 with respect to r n but this term is approximated with the gradient at a monomer's center-of-mass position, i.e., through Strictly speaking, this part of the model constitutes a force field rather than a potential. However, we 198 expect differences between an "exact" treatment and the pursuied approach to be small. Moreover, the 199 main goal of the model is to obtain a scaling similar to that predicted by Flory, i.e., R EE ∝ P 3/(D+2) 200 in the limit of large P when a single polymer is immersed in a good solvent. This scaling (or more 201 precisely the one with the exact rather than the Flory exponents [12,30]) will most certainly be achieved 202 in our model, as it penalizes energetically a recurrent random walk. When parameterizing the model 203 to a specific polymer-solvent system χ needs to be set such the correct prefactor is obtained.

204
Finally, the smearing out of density proposed in Eq. (24) was done in an approximate fashion.

205
First, a grid was used to which the effective monomer density ρ(r) = ∑ n δ(r − r n ) was binned 206 such that the center of mass remained unchanged. Rather than using a fast Fourier transform, the 207 density at a given grid point was smeared out such that one ninth of it was reassigned to the original 208 grid point as well as to its nearest, and next-nearest neighbors on a square lattice. The procedure 209 is repeated four times, resulting in a smearing out that starts to look Gaussian. Density gradients  Each bead can connect to a maximum of one end-monomer of a polymer, as depicted schematically in 234 Fig. 1(b). A connection is modeled with a harmonic spring of zero equilibrium length and a stiffness 235 k 1 , where k 1 pertains to the value used for the polymers.

236
As the crosslinkers were short molecules, no coarse-grained density was defined for them, 237 neither were they coupled to the polymer's coarse-grained density. Connections between crosslinkers 238 and polymer ends were established irreversibly after equilibration had taken place whenever a free 239 crosslinking monomer approached an unconnected polymer end within a 0 /3.

240
In our model the formation of crosslinks and crosslink topology evolves such that polymer ends 241 may end up dangling. There is also the possibility of clusters to be formed that do not belong to the 242 percolating cluster. As these constituents may also exist in reality, it was decided to not eliminate 243 dangling ends or disconnected polymers. In real systems they certainly also affect the viscoelastic 244 response functions until they are squeezed out of the hydrogel.

Model parameters for indentation simulations 246
For the remaining part of this paper, the following parameters were selected to be: P = 24, 247 a 0 = 1.5, k 1 = 16, k 2 = 4, and χ = 0.25. This leads to a persistance length λ p of roughly 16 segments 248 so that crosslinks were separated by polymers being approximately 1.5 times λ p . The crosslinkers 249 were modeled with P = 4, a 0 = 1.5, k 1 = 10, k 2 = 10, and χ = 0. Finally, the temperature was set to 250 T = 0.25.

251
The total number density of the monomers was set to unity, which resulted in a relatively dense 252 mesh, which allowed us to test if the network would swell after releasing it from the mold, through

Contact mechanics 266
The usual approaches to contact mechanics are based on the assumption that knowing the (normal) surface displacement u(r) of an elastic body, which is imposed through a frictionless contact with a counterbody, is sufficient to deduce the deformed body's elastic energy. In other words, U ela is a functional of u(r). For flat, linearly elastic, semi-infinite bodies, this functional is most easily expressed in a Fourier representation where E * is the contact modulus and q a wavevector, or, in our case of a 1+1 dimensional hydrogel, a 267 wavenumber. The validity of Eq. (26) furthermore assumes isotropy, and the surfaces to be extended, 268 or, more precisely, to be periodically repeated. More so, conventional linear elasticity is a scale-free 269 theory.

277
The first generalization is the use of a scale-dependent elastic modulus, in which case the following substitution needs to be made It allows two things to be achieved: First, the scale-dependent stiffness of (bulk) elastomers can be

283
In addition, for deformations that are not frictionless at the primary interace, the surface displacement 284 must be represented as a vector and E * (q) as a tensor of rank two.

285
E * (q) of a crosslinked network is expected to have no significant wavelength dependence on q, 286 once 1/q cleary exceeds the distances on which the hydrogel is heterogeneous, say, a few times the 287 typical distance between crosslinks. More generally speaking, 1/q needs to excced any characteristic 288 length in the system, such as λ p , however, it must not exceed the hydrogel's thickness, in which case 289 E * (q) is sensitive to the boundary condition -e.g., constant stress versus constant displacement -on 290 the opposite surface. In a non-crosslinked melt E * (q → 0) → 0, since the (static) contact modulus of a 291 melt is zero in the continuum limit. house-written C++ code, which was specifically designed for this study. In our MD simulations, 295 the end points of the segments are treated as dynamical variables rather than the segments themselves.

301
In addition to the conservative forces acting on the monomers, whose enumeration we start 302 at n = 0 rather than at n = 1 as for the bonds or segments. Constant temperature needs to be imposed, which was achieved with thermostats. An ideal thermostat emulates the (dynamical) effects 304 of those degrees of freedom (DOFs) that have been coarse-grained out on the remaining DOFs. Here, 305 however, we are predominantly interested in (static) thermodynamic properties, which is why a quick 306 equilibration is desired. Towards this end, we use a Langevin thermostat, which consists of a damping 307 term −m γ v n , where v n is the velocity of monomer n and a random force with mean zero, and a 308 second moment of 2 m γ k B T/∆t, where ∆t is the time step. The damping constant is best set such that 309 the slowest internal mode in a polymer, namely R EE , is close to being underdamped, which can be 310 reasonably well achieved by chosing γ = k eff /(P m). consequence of the structural and thus elastic heterogeneities that exist at small scales.

319
As can be appreciated in Fig. 7, the studied in-silico hydrogel also undergoes non-affine  although this argument only applies to a manifold of fixed height; fixed shape would require more 335 elaborate considerations. Thus, for large contact sizes but thin slabs, the elastic manifold would be 336 pinned on astronomically large times, even though the walls are perfectly smooth microscopically.

337
From the discussion in the last paragraph, it might be tempting to conclude that the relation 338σ (q) = q E * (q)ũ(q)/2, is of rather limited benefit with respect to modeling the contact mechanics of  Essentially all solids stiffen when being compressed. The reason for simple, i.e., densly packed solids is that atoms are pushed more deeply into the repulsive parts of their potentials, where the potentials' curvature -and thus effective spring constants -increase. In the case of hydrogels, compression leads to an increase in the crosslink density, which is argued to linearly affect a hydrogel's shear stiffness at small crosslink densities [12]. To explore these effects in the pursued model, the mold is replaced with an adhesionless indenter having a sinusoidal shape at fixed thickness or height h 0 , i.e., h(x) = h 0 +h q cos(q x).
The value ofh q is set such that qh q = 1/4. By chosing qh q constant, it is ensured that the contact line -353 as it would be determined in a continuum approximation -is constant, which facilitates the deduction 354 of meaningful E * (q) relations. A relatively large value of qh q guarantees relatively good signal to 355 noise ratios.

356
In Fig. 8, the stress profiles emerging in response to a sinusoidal indenter with fixed surface 357 shape -specifically an indenter having a wavelength of 100 as the one depicted in Fig. 6 -is shown

366
The non-isotropic stress at an indenter-fluid interface revealed in Fig. 8(a) might appear 367 counterintuitive from a continuum perspective, because non-isotropic pressures imply a system's 368 ability to sustain non-zero shear. However, it needs to be kept in mind that any system having structure at a lengthscale λ in one way or another is not isotropic, which is why it can sustain moderate amount 370 of shear imposed at that or smaller lengthscales. In our specific case, stiff chains resist being squeezed 371 into an ondulation with a wavelength being less than the polymer's persistance length -whether or 372 not it belongs to a crosslinked network. A very small, non-isotropic stress would have even resulted 373 for a freely-jointed polymer model, in which case R EE would have been the characteristic lengthscale 374 below which it reveals noticeable "solid-like" behavior. The perhaps most simple argument to sustain 375 this latter claim is that long random paths cannot be squeezed into narrow structures or ondulations, 376 while short random paths would not object. 377 We will now consider the swelling of the hydrogel, i.e. its expansion upon reducing the 378 confinement by moving the ondulated indenter to a position allowing the hydrogel to double its 379 mean original width. Solvent is implicitly supplied in an unlimited fashion by keeping χ unchanged.  The effect of solvent quality on the stress profiles at mesoscopic length scales is studied next. Parameters remain similar as before with a focus on the crosslinked (ρ cl = 0.83) structure. However, this time, a hydroegel in Θ-solvent (χ = 0) is contrasted with that in a good solvent (χ = 0.25). Results are summarized in Fig. 9. Monomers are depleted near the walls for the Θ solvent compared to good solvent conditions, as shown before in Fig. 5(a,b) for two plane-parallel walls. Moreover, the better solvent increases the osmotic pressure with which the polymers push against the sinusoidal piston. It thus leads to a larger mean stress but also to a larger oscillatory response in the stress profile, which can be representad quite accurately as a truncated Fourier series via σ(x) = 2 ∑ n=0 σ n cos(n q x).
(29) with a qualitative analysis, which shows that the relevance of the first higher-harmonic σ 2 compared to 399 the ground harmonic σ 1 is almost twice as large for the Θ solvent than for the good solvent, i.e., 0.19 400 versus 0.11. Having a "softer" system show larger anharmonicity than a more rigid system at identical 401 displacement and wavelength can be seen as counterintuitive from more convential elasticity. the indenter's wavelengths is less than, say, λ p /2, which, in this study, also corresponds roughly to 413 the mean distance between crosslinks. Conversely, the disconnected or partially connected polymer 414 melts show a quickly decaying E * (q) with decreasing q, as to be expected for a melt. Connected symbols consider hydrogels, which had been crosslinked in a "hydrophilic" mold attracting the solvent but repelling the polymers, while open symbols show simulation results for hydrogels that were crosslinked in a "hydrophobic" mold leading to a weak adhesion to the polymers. All hydrogels were exposed to the same indenter at a given value q and the mean height was kept constant throughout all simulations. The hydrophobicity of mold and indenter were also kept identical. The default system, shown in Fig. 6, corresponds to an undulation with q = 2π/100 ≈ 0.063. the polymers during the crosslinking stage. The arguably strongest aspect of the model is that 420 its description becomes more reliable with increasing degree of polymerization P and increasing 421 distance between crosslinking points, because it is very much designed in the spirit of a self-consistent 422 field theoretical description of polymers. The latter become exact in the limit of large P. While 423 revealing many trends correctly, which are summarized further below, the current model certainly 424 needs refinement when applied to specific applications.

425
First, it needs to be generalized to three dimensions, which, however, is an exercise that can 426 be achieved in a rather straightforward manner. What makes simulations potentially simpler in 3D 427 compared to 2D is that short-range repulsion can be introduced whereby bond crossing is prevented 428 or at least made extemely unlikely without leading to artifacts as it would happen in 2D. Second, 429 parameters for k 1 , k 2 , and χ must be properly gauged, for example, from measurements of R EE in after averaging, at least within -but as we believe probably also beyond -the linear-response regime.

457
In future models, smaller bond stiffnesses and longer chains between crosslinks will be studied to 458 (better) separate the scales and effects due to crosslinking and persistance length.

459
The large (local) stress sensitivity of E * is certainly important to be taken into account in the 460 modeling of the contact mechanics of polymers, because it leads to an important feedback mechanism:

461
If a contact-mechanics description is conducted assuming linear elasticity, some zones in the interface 462 will carry more load than others. The hydrogel stiffens below these points, which means they carry 463 even more load, which leads to further stiffening. Once pressures are sufficiently large to induce 464 more than a 20% or 50% correction to the original E * (q), the linear-response assumption can no 465 longer be justfied. For randomly rough surfaces, this is expected to happen at the point where the require major modifications to them to make them reliable down to small scales. In contrast, we 469 foresee that Persson theory can be readily applied to the systems in question but require a reliable 470 determination of E * (p, q) across the scales for a study along these lines, which is planned for the future.