Daniele Loco

Étienne Polack

Stefano Caprasecca

Louis Lagardère

Filippo Lipparini
email: flippari@uni-mainz.de

Jean-Philip Piquemal

Benedetta Mennucci
email: benedetta.mennucci@unipi.it

A QM/MM approach using the AMOEBA polarizable embedding: from ground state energies to electronic excitations

A fully polarizable implementation of the hybrid Quantum Mechanics/Molecular Mechanics approach is presented, where the classical environment is described through the AMOEBA polarizable force field. A variational formalism, offering a self-consistent 1 relaxation of both the MM induced dipoles and the QM electronic density is used for ground state energies and extended to electronic excitations in the framework of Time-Dependent Density Functional Theory combined with a state specific response of the classical part. An application to the calculation of the solvatochromism of the pyridinium N-phenolate betaine dye used to define the solvent ET30 scale is presented.

The results show that the QM/AMOEBA model not only properly describes specific and bulk effects in the ground state but it also correctly responds to the large change in the solute electronic charge distribution upon excitation.

Introduction

The idea of studying an intrinsically quantum-mechanical (QM) process taking place in a complex environment by partitioning the whole system into a (smaller) subsystem (S) and the environment (E) has a long history in quantum chemistry. Within this framework, S is identified as the minimal subunit where the process occurs, and is described using an accurate albeit expensive level of theory, typically QM, whereas the environment is treated at a much lower detail, but sufficient to properly describe its effects on the process under study. Two alternative strategies are most commonly followed: one can either model the environment as a polarizable continuum 1,2 or employ an atomistic description of E, but introducing a Molecular Mechanics (MM) force field (FF). [3][4][5][6][7] The latter class of approaches is commonly referred to as QM/MM, and is widely used, particularly when electronic processes in complex environments are studied.

Since the first formulations, different QM/MM approaches have been proposed; in most of them the effective Hamiltonian defining the electronic properties of the QM subsystem in the presence of the MM system, is divided into a term describing the isolated S, and a term taking into account its electrostatic interaction with the classical environment through point charges centered on the MM atoms. This formulation, known as "electrostatic embedding" is nowadays the most common QM/MM formulation: it in fact allows one to include the effects of the classical subsystem in the determination of the QM electronic density and all its related properties. What, however, this approach still misses is the possibility for E to polarize in response to the charge density (and its eventual changes) of the S subsystem. To include such mutual polarization effect, an extension of the QM/MM formulation beyond the electrostatic embedding is necessary and different alternative strategies have been proposed so far. [8][9][10][11][12][13][14][15][16][17][18][19][20] In this work, we will illustrate the theoretical development and the computational implementation of a novel polarizable QM/MM approach, based on the AMOEBA polarizable FF, [21][22][23] which we will refer to as QM/AMOEBA. Within this framework, the environment polarization is achieved through the use of atomic point dipoles, which are induced as a response to the electric field generated by the S system as well as the same MM sites bearing distributed multipoles up to quadrupoles. An approach to couple AMOEBA to a QM package has been recently proposed 24 in the framework of a not self-consistent procedure. In this contribution, for the first time in the context of AMOEBA, the global relaxation of the mutual QM/MM polarization is achieved by solving the QM/AMOEBA equations in a self-consistent manner, without approximations.

The implementation is further extended to describe electronically excited states. As a matter of fact, non polarizable QM/MM formulations can be easily extended to excited states as the "new" operator has a one-electron nature. When a polarizable embedding like AMOEBA is used instead, a difficulty appears when a linear response (LR) approach as that commonly used in the Time-Dependent Density Functional Theory (TD-DFT) is used. In the LR formulation, in fact, the whole spectrum of the excitations of interest is determined in a single step calculation by solving for the poles of the proper response function. When the same problem is formulated within a polarizable embedding, an additional contribution has to be taken into account, namely the dynamic response of the E subsystem. In the standard formulation of polarizable models (both in their continuum and MM formulations) such a response is calculated through the transition densities corresponding to the different excitations: the oscillating transition density of the S subsystem induces an oscillating polar-ization in the E subsystem which creates an in-phase response acting back to the transition density. 25,26 This LR formulation has been shown to properly describe the dynamic environment effect in excitations involving bright states characterized by a large transition dipole moment. 27 The same formulation, however, lacks the capability of describing the relaxation of the environment in response to the changes in the QM density upon excitation: it is therefore not suitable to model excitations involving large changes such as Charge-Transfer (CT) like excitations. To overcome this shortcoming, various models have been introduced to recover a state-specific (SS) description of the response both within a continuum [28][29][30] or an MM formulation. 20 Here we adopt a perturbative correction which is exactly equivalent to the so-called corrected Linear Response (cLR) scheme originally developed for polarizable continuum models. 28 Within this framework, a relaxed density matrix is calculated for the excited state of interest and the corresponding excitation energy is corrected for the interaction with the proper induced dipole moments within the environment.

To show the potentialities of the method, we have selected a very well-known solvatochromic probe, namely the pyridinium N-phenolate betaine dye (or 'betaine-30'). Due to its large negative solvatochromism combined with a high solubility in many different solvents, it was proposed as a solvatochromic indicator for the determination of the solvent polarity, through an empirical scale called E T (30). 31,32 From a computational point of view this molecule is really challenging for at least two reasons: (i) it is zwitterionic in its ground state and (ii) the excitation is related to an intramolecular charge transfer (CT) between the pyrimidine and phenolate moiety which acts to reduce the zwitterionic character. For such a kind of excitation, the solvation model not only needs to be able to properly describe both specific and bulk effects in the ground state but it also has to correctly respond to the large change in the electronic charge distribution from the ground to the excited state. We will show that both these two requirements are fully satisfied by our implementation of the QM/AMOEBA.

Theory and Implementation

The AMOEBA Force Field

One of the main characteristics of AMOEBA is the improved description of the electrostatics through atomic multipoles, up to the quadrupoles, placed on each classical atom. Polarization effects are included by using an induced point-dipole model within the smeared Thole's damping interaction scheme. 33 This is achieved by also providing the classical atoms with atomic polarizabilities. In general, the polarizability α i is a symmetric matrix; in practice, the scalar isotropic polarizabilities α i ≡ 1 3 tr(α i ) are used instead. In a purely classical framework, at each polarizable site i, the static multipolar distribution generates an electric field E i that induces an electric point-dipole moment, µ i . The set

{ µ i } N Pol
i=1 is the unknown of the polarization problem and the minimizer of the functional 34

E Pol = - N Pol i=1 E α i µ α i + 1 2 N Pol i=1 (α -1 i ) αβ µ α i µ β i + 1 2 N Pol i=1 N Pol j =i T αβ ij µ α i µ β j , (1)
where indexes α and β generally indicate cartesian components, for which the Einstein's summation convention is assumed; E i and µ i are the electric field and the induced dipole, both at the atom site i. Note that we will generally assume that, while all the classical sites are characterized by multipole moments, not all of them will be polarizable, i.e., associated with a polarizability. The number of classical and polarizable sites are indicated as N MM and N Pol , respectively, with N MM ≥ N Pol . Further details can be found in previous works on the force field; 23,34 here, we only recall the expression for the (damped) dipole interaction tensor:

T αβ ij = - δ αβ | r i -r j | 3 λ 3 (u ij ) + 3 | r i -r j | α | r i -r j | β | r i -r j | 5 λ 5 (u ij ). (2)
Several damping schemes to avoid the well known polarization catastrophe can be found in the literature. 33,[35][36][37] The scheme employed in AMOEBA makes use of an exponential damping:

       λ 3 (u ij ) = 1 -e -au 3 ij λ 5 (u ij ) = 1 -(1 + e -au 3 ij )e -au 3 ij . (3)
where

u ij = r ij /( α i α j ) 1 6
is the effective distance between two polarizable sites i and j.

In a more compact notation,

E = E T 1 , E T 2 • • • E T N Pol T and µ = µ T 1 , µ T 2 • • • µ T N Pol
T are the collections of the electric field and induced point dipoles at each polarizable site, while the symmetric, positive definite matrix T , the size of which is 3N Pol × 3N Pol and which is usually called polarization or interaction matrix, is defined as in Equation ( 4), where (T ) ij is the 3 × 3 matrix defined as in Eq. ( 2).

T =          α -1 1 T 12 . . . T 1N P ol T 21 α -1 2 . . . T 2N P ol . . . . . . . . . . . . T N P ol 1 T N P ol 2 . . . α -1 N P ol          (4)
Equation ( 1) can now be rewritten in matrix form as

E P ol = 1 2 µ † T µ -µ † E (5)
It is possible to find the polarization energy as the minimum of E P ol with respect to the induced point dipoles, 34 which corresponds to the solution to the linear system

∂E P ol ∂µ = T µ -E = 0. ( 6
)
This is the general formulation of the polarization problem, and is clearly variational. In the case of AMOEBA, however, the polarization energy is no longer a variational functional of the induced dipoles. This is due to the fact that two sets of induced dipoles are generated by two electric fields, differing for the scaling of the local interactions. Particularly, one set of dipoles, indicated as µ d , is due to the so called "direct field" E d , which is produced by multipoles placed on all the other classical sites. The other set, indicated as µ p , is induced by a "polarization field" E p , where the contribution of the 1-2, 1-3 and 1-4 neighbors is scaled. A more detailed treatment on this point can be found in other works of some of us. 38 Here it suffices to say that the previous expression can be recast taking into account the different fields and induced dipoles. We recall that

E Pol A = - 1 2 µ † d E p , (7)
is the proper expression for the AMOEBA polarization energy. Imposing the stationarity conditions for both sets of dipoles

       T µ d = E d T µ p = E p , (8)
where each set of induced dipoles is the variational minimizer of the corresponding energy functional, it is possible to reformulate the AMOEBA polarization energy as the combination of three variational expression, finally obtaining

E Pol A (µ d , µ p ) = 1 2 µ † d T µ p - 1 2 µ † p E d + µ † d E p (9)

An SCF-QM/AMOEBA formulation

The coupling of a classical description of the environment with SCF-based methods has already been discussed, especially for polarizable continuum solvation models in a variational scheme. 39,40 Generalizing those formulations to the QM/AMOEBA approach, the global variational energy functional can be written as the sum of three terms: (i) a purely QM energy functional, i.e., the SCF energy functional E QM , (ii) a purely MM term given by the sum of the bonding, dispersion/repulsion, electrostatic and polarization terms according to their definition within the AMOEBA FF; and (iii) a coupling term; namely we have

E(P , µ) = E QM (P ) + E MM (µ) + E Coup (P , µ) = E QM (P ) + E Env (P , µ), (10)
where in the RHS of the equation we added together the coupling and MM energy functionals in E Env . Notice how the variational strategy, by introducing the global energy functional in Eq. ( 10), automatically takes into account the mutual polarization effects between the QM density and the induced dipoles. 39,41 The variational environment term, E Env , can be written as the sum of a constant MM contribution E FF +E (0) which does not depend on either the electronic density or the induced dipoles, a polarization energy, and a QM/MM coupling part. For the following manipulations, it is convenient to write separately the interaction of the QM density with the static multipoles as E QM/MM (P), and to add the interaction between the induced dipoles and the QM density to the polarization energy, from now on named E Pol (P, µ); the resulting expressions of the various terms are:

E Env (P , µ) = E FF + E (0) + E Pol (P, µ) + E QM/MM (P ) (11)
E Pol (P , µ) = 1 2 µ † d T µ p - 1 2 µ † p E d + µ † d E p - 1 2 µ p + µ d † E QM (P) (12)
E QM/MM (P ) = q † V QM (P ) -µ † s E QM (P ) -Θ † ∇E QM (P ). ( 13
)
where q i , µ s,i and Θ i are the fixed charges, dipoles and quadrupoles, respectively.

The electronic terms of E QM/MM in Eq. ( 13) can be expressed as functions of the one-particle electron density matrix elements:

                       V QM ( r i ; P ) = N b µν P µν V µν ( r i ) = - N b µν P µν R 3 χ µ ( r)χ ν ( r) | r -r i | d 3 r E QM ( r i ; P ) = N b µν P µν E µν ( r i ) = N b µν P µν R 3 χ µ ( r)χ ν ( r)( r -r i ) | r -r i | 3 d 3 r ∇ E QM ( r i ; P ) = N b µν P µν G µν ( r i ) = N b µν P µν R 3 χ µ ( r)χ ν ( r)J E( r i ) d 3 r. ( 14
)
where the electron density P has been expanded in a basis of atomic orbitals and J is the Jacobian of the electric field with elements (J E) αβ ( r) = ∂E α ( r) ∂r β , being α and β generic Cartesian components. The corresponding terms induced by the nuclear distribution are trivial and are not reported.

Imposing the stationarity conditions for the global functional in Eq. (( 10)), taking into account the constraints on the electronic density matrices, gives the coupled QM/AMOEBA equations.

The QM/AMOEBA Fock matrix F is obtained by differentiating Eq. 10 with respect to the density matrix:

F = ∂E(P , L, µ) ∂P = ∂E 0 (P ) ∂P + ∂E Env (P , µ) ∂P = F (0) + F Env . ( 15
)
where the term F (0) of ( 15) corresponds to the Fock matrix of the pure QM electronic problem, while F Env is the contribution to F from the classical environment. The elements of the latter can be written as

F Env µν = q † V µν -µ † s E µν -Θ † G µν - 1 2 (µ p + µ d ) † E µν ( 16
)
where the matrices V µν , E µν and G µν are those appearing in Eq. ( 14). Note that the E µν elements involved in the second term on the RHS of Eq. ( 16) are formally identical to those appearing in the first one, but are computed over a different set of atomic sites (the polarizable ones only).

The linear equations for the dipoles are obtained by setting the derivatives of Eq. (10) with respect to both sets of dipoles to zero:

T µ p = E p + E QM (P ) T µ d = E d + E QM (P ). ( 17
)
Notice that the Fock matrix depends on the induced dipoles and that the RHSs of the dipole equations depend on the density matrix: the two sets of equations are therefore coupled, accounting correctly for the mutual polarization of the QM and classical charge densities. The coupled equations need to be solved iteratively, which is not a problem in practice as the SCF equations are already solved with an iterative algorithm. We note that the matrix T depends only on geometrical parameters, and can be computed and inverted at the first SCF cycle, thus reducing the computational requirements at the following steps.

Alternatively, this problem can be also efficiently solved with an iterative procedure (see Section 2.4 for more details).

Extension to electronic excitations

The discussion so far has involved the energies and properties of the electronic ground state.

To extend the treatment to excited states, we follow a linear response theory for SCF methods. A complete derivation of the polarizable QM/MM LR response equations can be found elsewhere, 14,41 here it suffices to say that the electronic transition energies ω K and densities X K , Y K are found by solving the modified Casida's equations 42

   A B B A       X K Y K    = ω K    1 0 0 -1       X K Y K    . (18)
The matrices A and B are defined as

A ai,bj = δ ab δ ij ( a -i ) + (ai|bj) -c x (ab|ij) + c l f xc ai,bj + C Pol ai,bj (19)
C Pol ai,bj = - N Pol P R 3 φ a ( r) ( r -r p ) | r -r p | 3 φ i ( r)d r • µ T P (φ b , φ j ). ( 20
)
The matrix C Pol depends on the dipoles µ T induced by the transition density

P T K = X K +Y K ,
which are obtained by solving, for each couple of transition vectors X K , Y K , the response linear equations

Tµ T = E(P T K ), (21)
where

E P (P T K ) = ai P T K,ai R 3 φ a ( r) ( r -r p ) | r -r p | 3 φ i ( r)d r. (22)
Notice that in the right hand side of eq. 21 only the electric field due to the transition density appears: as the classical multipoles do not contribute to the transition dipoles, the p and d response dipoles introduced in section 2.1 are identical. Again, the modified Casida's equations depend (linearly) on the transition induced dipoles, which in turn depend (linearly) on the transition vectors X, Y: although the coupling is linear, it is still convenient to solve the couple equations iteratively, by solving the transition dipoles equations at each step of the iterative procedure used to solve Casida's equations.

The great advantage of the LR scheme is that it allows one to obtain a whole spectrum of transition energies at once. If the system is isolated, it provides the same results than state specific (SS) approaches, where instead the wavefunction of each excited state is explicitly computed. By contrast, the introduction in the Hamiltonian of a non-linear term due to the presence of a polarizable environment results in a formal difference between the two approaches. 25,43 For this reason, a correction to the LR scheme (called corrected Linear Response, c-LR) has been proposed. 28 The difference between the two approaches (LR and c-LR) can be depicted by a two-step process for excitation in solution: first the molecule in its ground state, in equilibrium with the solvent, is excited to the state K in the presence of a solvent response frozen to the one of the solute ground state. The description for this step of the process is the same in both the approaches. In the second step, the fast electronic degrees of freedom of the solvent equilibrate with the electron density of the solute excited state, and in this part of the process the LR framework lacks in taking explicitly into account the energy variation accompanying this relaxation, whereas it accounts for a correction which, being originated by the dynamic solute-solvent interactions, can be classified as a part of dispersion. The working equation, namely

ω c-LR K = ω 0 K - 1 2
N pol P µ P ( r P ; P ∆ K ) E( r P , P ∆ K ) (23) corresponds to the expression of the c-LR transition energy for the K-th excited state. In this formulation ω 0 K represents the response of the system in a solvent frozen in its initial configuration (corresponding to the system in its ground-state), and it is the solution of a non-Hermitian eigensystem as in Eq. ( 18), where the C Pol ai,bj = 0 but the orbitals and their corresponding energies, employed to build the A and B matrices have been obtained by solving the SCF equation for the solvated system. The P ∆ K is the so called relaxed-density matrix, computed through the so-called Z-vector approach 44 as

P ∆ K = P T K + Z K (24)
where P T K is the unrelaxed density matrix with elements given in terms of the transition vectors |X K , Y K , and the Z-vector contribution Z K accounts for orbital relaxation effects.

Implementation

The hybrid QM/AMOEBA method has been implemented in a development version of the Gaussian software. [START_REF] Frisch | Gaussian 09 Revision A.1[END_REF] In the present implementation, the program computes all the terms that depend on the QM density, including the sets of induced dipoles. The E F F + E (0) energy terms depend neither on the QM density nor on the induced dipoles and, in a static picture, are constant.

In order to compute the induced dipoles, both an inversion approach and an iterative one have been implemented. In the former case the polarization matrix (Eq. 4) is computed and inverted only once during the first step of the SCF, then stored and multiplied for the appropriate field to obtain the solutions to the polarization problem. At each SCF cycle only two matrix-vector products (one per set of dipoles) need to be performed, concentrating the computational effort almost entirely in the first SCF step.

By contrast, in the iterative approach, the polarization equations are solved iteratively at each SCF step or at each step of the iterative solution of Casida's equations, using as a guess the induced dipoles of a previous step, when available. Two algorithms have been implemented: Jacobi Iterations coupled with direct inversion in the iterative subspace (JI/DIIS) and the Preconditioned Conjugate Gradient (PCG). 34,[START_REF] Caprasecca | [END_REF] Both are available for the ground and LR excited state problems.

The advantage of the inversion approach is that the most expensive part of the procedure, i.e., the matrix inversion, needs only be performed once. However, the computational cost of such an operation grows with the cube of the number of polarizable sites: for large systems, iterative alternatives become mandatory. Furthermore, iterative methods can be easily combined with fast summation techniques, such as the Fast Multipole Method, [START_REF] Caprasecca | [END_REF]47 in order to achieve linear scaling in computational cost for the solution of the polarization equations.

We here present the relative time requirements for the computation of the induced dipoles in the whole SCF procedure, using the inversion and iterative algorithms. The calculations have been performed on the test system whose results will be discussed in Section 3. The dimension of the polarization problem has been artificially varied, by considering various polarization cutoffs, i.e., by including a different number of induced dipoles. The dimension of the QM system, as well as that of the non-polarizable classical environment, have been kept fixed. Figure 1 reports the dependence of the relative times on the number of induced dipoles. From the plots reported in the figure, it is clear that the iterative procedures are more convenient than the inversion one. This becomes more evident as the size of the problem increases, since both JI and PCG algorithms scale quadratically with the number of induced dipoles, while the inversion algorithm, LU decomposition, is characterized by a cubic scaling.

In particular, when the polarization radius is increased from 7 to 20 Å (number of dipoles increasing by ∼ 10× from 765 to 7338), the inversion algorithm increases by ∼ 300×. For both JI and PCG algorithms, the relative increase is ∼ 60×.

Furthermore when the number of polarizable dipoles is small (polarization radius 7 Å, 765 dipoles), even if the JI/PCG methods are faster than the inversion one, the choice of the solution scheme does not really matter since the relative time to solve for the induced dipoles with respect to the total time is, for each of the three cases, less than the 0.1%. This is no longer true in the most expensive case (polarization radius 20 Å, 7338 dipoles), where almost 16, 5% of the time for the calculation is spent in computing the induced dipoles with the inversion method (the iterative solutions take less than 1% of the total time).

A test case: the solvatochromism of betaine-30

In this section we present an application of the QM/AMOEBA implementation to the simulation of the excitation energies of a well-known solvatochromic probe, the pyridinium N-phenolate betaine dye, from now on indicated as "betaine-30" (see Figure 2). The solvent we selected for this test case is water, since it constitutes one of the most interesting cases due to its high polarity combined with an hydrogen-bonding character. In particular, we will try to dissect these two components of the solvent effect by comparing QM/AMOEBA with a continuum description using the Polarizable Continuum Model (PCM) within its Integral Equation Formalism (IEF). 48 We recall that this model describes the environment as a structureless continuum, characterized by its macroscopic dielectric function. A cavity containing the solute is built around it, and the solvent polarization as a response to the solute charge density is represented by an induced surface charge distribution on the cavity.

In order to obtain a realistic sampling of the solute-solvent interaction, a well-established procedure has been applied, which makes use of snapshots extracted from a classical molecular dynamics (MD) simulation with the AMBER General Force Field for organic molecules.

An optimized geometry of the dye was solvated with a cubic box containing 11380 water molecules, described at TIP3P level, 49 with dimensions 74 × 72 × 78 Å. The system was heated from 0 to 300 K for 100 ps with a 2 fs time step, employing the SHAKE algorithm 50 as in all the following steps. The Berendsen thermostat 51 was used to control the temperatue. A 200 ps equilibration in the NVT ensemble was then performed. Before the actual production, the system was further equilibrated in the NPT ensamble at 1 atm for 5 ns.

Here, the Monte Carlo barostat implemented in Amber14 52 was employed. In the 30-ns long MD simulation, the betaine-30 was kept frozen in its ground state equilibrium geometry (computed at QM/PCM level). A set of 100 uncorrelated snapshots was extracted, on which the QM/AMOEBA calculations were performed. Since the effect of polarization is short-ranged, and the computational cost increases markedly with the number of induced dipoles, only the classical atoms within a certain radius (the polarization radius, R pol ) were allowed to polarize. The optimal value of the radius was chosen after performing a convergence test, where the first three excited states of the solvated betaine dye were computed at increasing values of R pol . The results, reported in Figure 3, show that convergence can be observed starting from 10-12 Å. A safe value R pol = 15 Å was chosen and employed in all the calculations presented. Blue, red and yellow curves refer to the first, second and third excited states, respectively, calculated on one snapshot at TDDFT level, within the linear response approach. A value of R pol = 0 indicates that all the solvent molecules are non polarizable, and are only described in terms of fixed multipoles.

All calculations (geometry optimizations and excitation energies) have been performed at the (TD)DFT level of theory using the CAM-B3LYP exchange-correlation functional 53 together with the 6-31+G(d,p) gaussian basis set. We particularly focused on the properties of the lowest (bright) excitation, on which the E T (30) scale is based. The excitation shows a strong charge-transfer character, and for this reason the corrected linear response (c-LR) approach 28 is expected to give a more accurate picture than the standard LR.

The results obtained with the three different QM/AMOEBA responses described in section 2.3 (namely ω 0 , LR and c-LR) are summarized in Figure 4. From the comparison, it is evident that including the effects of the relaxation of the solvent response through a c-LR formulation a significant redshift (of about 0.2 eV) is found with respect to the frozen solvent approximation (ω 0 ). As expected due to the CT character of the excitation, the LR formulation does not change the ω 0 value as it cannot account for the effects of the rearrangement of electron density undergone by the dye upon excitation.

The comparison with experimental data is, however, the most interesting test. In order to be able to provide an analysis which is not biased by systematic errors in the excitation energies due to the QM level of theory, we compare the computed and experimental gas-towater solvatochromic shifts instead than the absolute excitation energies. To better elucidate the various effects that determine the observed solvatochromic shift, i.e., short-range and specific interactions and the bulk effects, it is useful to compare with two different solvation models. The first, QM/PCM, employs a purely continuum description of the solvent whereas the second, QM(ME)/PCM, includes a "minimal environment" in the QM system, i.e., the two water molecules hydrogen-bonded to the oxygen of the dye, while the rest of the solvent is still treated at PCM level. The latter model should better describe the combination of short-range specific and bulk interaction with respect to the QM/PCM analog. In the QM(ME)/PCM calculation the configuration of the QM water molecules has been optimized at CAM-B3LYP/6-31+G(d,p). All calculations refer to c-LR responses. as the use of a calculated vertical excitation and an experimental band maximum, and the extrapolated value used for the experimental gas-phase value. However, the additional +0.4 eV shift obtained going from the QM(ME)/PCM to the QM/AMOEBA model clearly suggests that the strong hydrogen-bonding effect is not the only source of differences with respect to a purely bulk description. Looking more deeply into the results, it is interesting to point out that the difference between QM/AMOEBA and QM(ME)/PCM solvatochromic shifts is already present at ω 0 level, which means that the addition of the c-LR response is almost equally described by the two models. We can thus speculate that the zwitterionic nature of betaine-30 ground state has a large, specific, long range effect in orienting the water molecules and this leads to a final solvation effect that cannot be accurately described with a pure continuum or a partial atomistic/continuum approach.

Conclusions

We presented here the theoretical development and the computational implementation of a polarizable embedding QM/MM within the AMOEBA framework. The implementation accounts for a fully relaxed evaluation of the QM/MM energy both for ground and excited states, in the framework of a DFT/TD-DFT theory. In particular for the excited state, a state-specific (SS) formulation of the response of the classical part of the system has been used following the c-LR approach originally developed within continuum models. 28 The test case of betaine-30 is presented and discussed in comparison with purely continuum and mixed atomistic/continuum models: the obtained results in the simulation of the large gas-to-water solvatochromism of this dye, show that the QM/AMOEBA approach can describe with equal accuracy the effect of water on the zwitterionic ground state and the modification induced by the excitation.

This work represents the first step of a series of new developments towards a highperformance, parallel implementation of polarizable QM/MM molecular dynamics, which stems from the recent developments in the implementation of AMOEBA for classical MD simulations in the newly developed Tinker-HP 34,55-57 code and with the implementation of a versatile and transparent Tinker HP/Gaussian interface for energy, gradients and properties.

The implementation of analytical gradients for the QM/AMOEBA model and the development of the Gaussian/AMOEBA interface will allow us to exploit not only parallelism, but also linear scaling techniques [START_REF] Caprasecca | [END_REF] in order to further reduce the overall computational cost due to the polarizable embedding. The addition of a further layer to the model, namely, of a polarizable continuum solvation model, will finally allow to deal with long range, bulk effects while, at the same time, reducing the portion of environment to be treated explic-itly. An innovative, parallel, linear-scaling implementation of the Conductor-like Screening Model 38,58,59 has already been coupled to AMOEBA for classical MD simulations 38 and to the MMpol dipole-based polarizable force field for QM/MM calculations: 46 the combination of all the aforementioned developments will result in a powerful and efficient tool to investigate dynamic properties and reactivity in complex environments.

B

  ai,bj = (ai|bj) -c x (aj|ib) + C Pol ai,bj , where are the orbital energies, (pq|rs) are two-electron integrals in Mulliken notation and we assume that the orbitals are real and the coefficients c x and c l define whether we are considering the Hartree-Fock theory (c x = 1, c l = 0), pure DFT (c x = 0, c l = 1) or hybrid DFT. The elements of matrix C Pol are due to the environment polarization:

Figure 1 :

 1 Figure1: Relative times needed to solve for the induced dipoles, using different methods. Green line: inversion; Red line: Jacobi Iterations; Blue line: Preconditioned Conjugate Gradient. The times are cumulative and refer to the sum of all SCF cycles until convergence is reached. The number of SCF cycles considered is the average one across the various structures considered. The values have been normalized for comparison with respect to the time required to solve the smallest polarisation problem (7 Å polarization radius, 765 induced dipoles) with the inversion procedure.

Figure 2 :

 2 Figure 2: The betaine dye studied.

1 Figure 3 :

 13 Figure3: Lowest three excitation energies of betaine-30 as functions of the polarization radius, R pol , with the total radius for the inclusion of classical environment fixed at 25 Å. Blue, red and yellow curves refer to the first, second and third excited states, respectively, calculated on one snapshot at TDDFT level, within the linear response approach. A value of R pol = 0 indicates that all the solvent molecules are non polarizable, and are only described in terms of fixed multipoles.

Figure 4 :

 4 Figure 4: Distributions of the lowest (bright) excitation of betaine-30 calculated with the three different QM/AMOEBA responses (LR: blue line, ω 0 : yellow line, c-LR: red line). The curves are gaussian fits of the histograms.

Table 1 :

 1 Gas-to-water solvatochromic shift (in eV) for the lowest (bright) excitation of betaine-30, computed with different solvation models. The QM/AMOEBA value is obtained from the average over 200 snapshots extracted from the MD. The confidence interval is indicated. All calculations were performed at TDCAM-B3LYP/6-31+G(d,p) level.As expected, the QM/PCM model markedly underestimates the solvatochromic shift with respect to the two other models (and to experiments). This clearly shows that including an atomistic description of the most strongly interacting water molecules is fundamental to account for the differential solvation effects in the ground and the excited states of the

	QM/PCM	0.79
	QM(ME)/PCM	1.17
	QM/AMOEBA 1.58 ± 0.02
	Exp 54	1.56

betaine-30. However, it is worth noting that the inclusion of the two hydrogen bonded water molecules in the QM region in combination with a continuum description (through the QM(ME)/PCM model) still fails to recover a large portion of the observed solvatochromic effects (namely 0.4 eV are still missing to reproduce the experiments). The inclusion of a larger number of explicit solvent molecules with AMOEBA instead leads to a solvatochromic shift which is almost exactly equivalent to the experiments. Such an excellent agreement might be fortuitous due to the many approximations introduced in the comparison, such

Acknowledgement

This work was supported in part by French state funds managed by CalSimLab and the ANR within the Investissements dAvenir program under reference ANR-11-IDEX-0004-02. F. L. gratefully acknowledges the Alexander von Humboldt foundation for funding. E. P. and J.-P. P. are grateful for support by the Direction Generale de l'Armement (DGA) -Maitrise NRBC of the French Ministry of Defense. S.C. and B.M. acknowledge the European Research Council (ERC) for financial support in the framework of the Starting Grant (EnLight-277755).


Vito Acito

Michele Ciavarella

Alexis M Prevost

Antoine Chateauminois

Adhesive contact of model randomly rough rubber surfaces

Keywords: numbers: 46.50+d Tribology and Mechanical contacts; 62.20 Qp Friction, Tribology and Hardness Keywords: Adhesive contact, Randomly rough surfaces, Silicone rubber, Fuller and Tabor theory, JKR theory

We study experimentally and theoretically the equilibrium adhesive contact between a smooth glass lens and a rough rubber surface textured with spherical microasperities with controlled height and spatial distributions. Measurements of the real contact area A versus load P are performed under compression by imaging the light transmitted at the microcontacts. A(P ) is found to be non-linear and to strongly depend on the standard deviation of the asperity height distribution. Experimental results are discussed in the light of a discrete version of Fuller and Tabor's (FT) original model (Proceedings of the Royal Society A 345 (1975) 327), which allows to take into account the elastic coupling arising from both microasperities interactions and curvature of the glass lens. Our experimental data on microcontact size distributions are well captured by our discrete extended model. We show that the elastic coupling arising from the lens curvature has a significant contribution to the A(P ) relationship. Our discrete model also clearly shows that the adhesioninduced effect on A remains significant even for vanishingly small pull-off forces. Last, at the local asperity length scale, our measurements show that the pressure dependence of the microcontacts density can be simply described by the original FT model.

I. INTRODUCTION

Surface roughness has long been recognized as a major ingredient that cannot be ignored when considering adhesion and friction between macroscopic bodies. As emphasized by the pioneering work of Bowden and Tabor [1], a key component in the description of the tribological properties of rough bodies is the determination of the actual contact area formed by microasperity contacts of typical length scales distributed over orders of magnitude. This problem involves a still poorly understood interplay between the mechanical properties of the solid bodies, the geometry of the surfaces and adhesive forces. In the case of rubber-like materials, the effects of adhesive forcessuch as Van der Waals forces-on contact deformation have long been recognized with the pioneering work of Johnson, Kendall and Roberts (JKR) [2]. Based on the JKR model for the adhesion of a smooth sphere on a rigid flat, a first theoretical and experimental investigation of the effect of roughness on elastic contact with adhesion was carried out by Fuller and Tabor (FT) [3]. FT's model is based on Greenwood and Williamson's (GW) [4] description of roughness which considers identical spherical asperities of radius R whose summit heights is obtained from a Gaussian distribution with a standard deviation σ. Using JKR's solution for each asperity contact, FT obtained a solution for the rough adhesive contact problem, yielding an expression for the pull-off force which depends on the dimensionless parameter

θ F T = σ 3/2 E * R 1/2 w (1)
where E * = E/(1 -ν 2 ) is the reduced Young's modulus (ν being Poisson's ratio), R is the radius of the microasperities and w is the adhesion energy. According to FT's model, the pull-off force decays exponentially with θ F T . In the same work, FT also conducted experiments where molded rubber spheres were put in contact with a Perspex flat with a random roughness obtained by bead-blasting or abrasion. However, such randomly rough surfaces are too complex to be properly amenable to the derived adhesive contact model. In addition, the lens curvature itself raises some issues. Even though FT's experiments show that lens curvature has a negligible role on the pull-off force, one can still question how it would affect local microcontact sizes and pressure spatial distributions.

In this paper, we revisit FT's model for adhesive contacts of statistically rough surfaces in the light of experiments where we take advantage of recent micromilling techniques to engineer randomly rough silicone surfaces with prescribed distributions of both asperity lateral positions and heights. They consist of spherical asperities whose radii of curvature (≈ 100 µm) allow for an optical measurement of the spatial distribution of the microcontact areas when the patterned silicone substrate is indented by a smooth glass sphere. As compared to two previous studies [5,6] using similar surfaces, the effects of adhesion were here enhanced by reducing both the elastic modulus and the small scale roughness of the asperities. Such surfaces are reminiscent of the regular arrays of spherical rubber caps investigated by Verneuil and coworkers [7]. In this work, we get access to local quantities such as microcontact radii and contact pressures as a function of height distributions.

To analyze our measurements within the framework of FT's model, we do focus on adhesive contact close to equilibrium conditions using a contact loading procedure for which viscoelastic dissipative effects are minimized. Adhesive effects at equilibrium are discussed by comparing experimental results for the load dependence of the real contact area to predictions of an adhesive multiasperity contact model which takes into account the effects of the elastic coupling arising from both lens curvature and microcontacts interactions. We also show that the use of textured surfaces allows for an accurate discussion of FT's adhesive contact model based on our investigation of the statistical distributions of microcontact radii and pressures.

II. EXPERIMENTAL DETAILS

Normal contact experiments are carried out between an optically smooth glass lens (radius of curvature 103.7 mm, BK7, Melles-Griot) and a nominally flat Poly-DiMethylSiloxane (PDMS) slab decorated with micrometer sized spherical caps of equal radius of curvature, distributed randomly both spatially and in height. As detailed in [5,6], these patterned surfaces are obtained by crosslinking PDMS in Poly(MethylMethAcrylate) (PMMA) molds milled with a desktop CNC Mini-Mill machine (Minitech Machinary Corp., USA) using ball end mills of radius 100 µm. Prescribed random distribution of asperities were distributed over 1 cm 2 with a non-overlapping constraint and an asperity density of 2.10 7 m -2 . Three different patterns have been made on the same PMMA mold, one reference area covered with a square network of spherical caps having all the same height of 40 µm and two areas decorated with Gaussian height distributions of spherical caps with standard deviations σ = 5 µm and σ = 10 µm respectively. In addition, a fourth part of the PMMA mold has been kept smooth to allow measurements of the elastic modulus and the adhesion energy of the PDMS sample using a JKR contact configuration.

As shown in Fig. 1a, micromilled spherical caps on the PMMA mold inherently present a micrometric scale roughness which can induce a reduction of adhesion of the replicated patterned PDMS substrate. To minimize these effects, PMMA molds have thus been exposed to a saturated CHCl 3 vapor for 15 minutes. As a result of surface plasticization of the glassy acrylate polymer, surface tension effects were found to result in a smoothening of the surface of the spherical cavities of the mold (Fig. 1a). It was verified with optical profilometry that this treat- ment did not induce any change in the standard deviation of asperity height distribution. However, it resulted in a roughly 10% increase in the radius of curvature of the asperities up to 110 µm which was accounted for in both data analysis and contact simulations. Figure 1b shows for both random vapor treated PMMA molds, the Probability Density Functions (PDF) of the asperity height distributions measured over more than 400 asperities. PDFs are still well fitted by Gaussians with standard deviations very close to the prescribed ones, namely σ = 5.13 ± 0.19 µm and σ = 9.22 ± 0.15 µm.

PDMS substrates are obtained by crosslinking at 70 • C for 48 hours a mixture of commercially available Sylgard 184 and Sylgard 527 liquid silicones (Dow Corning, Midland, USA). As detailed in [8], mixing both silicones in various proportions allows tuning the elastic modulus in the [kPa-MPa] range. Here, Sylgard 527 and Sylgard 184 were mixed in a 0.75:0.25 weight ratio to achieve a Young's modulus of about 0.5 MPa. The elastic modulus E and the adhesion energy w of the flat PDMS substrate were obtained with standard JKR measurements with an optically smooth glass lens. These measurements were performed using the setup and the step-by-step contact loading procedure described below. They yielded E = 0.51 ± 0.11 MPa and w = 40.8 ± 2.1 mJ m -2 .

Normal contact experiments are carried out using a custom contact setup which is fully described elsewhere [9]. The glass lens indenter is fixed to a motorized vertical translation stage by means of a double cantilever beam with a stiffness of 290 N m -1 . During loading, the normal load is determined with a 0.1 mN resolution from the measurement of the deflection of the cantilever using a high resolution optical sensor (Philtec D64-L). In addition, the vertical position of the lens is monitored with submicrometer accuracy using a laser sensor (Keyence LK-H057). The PDMS substrate is fixed to two crossed motorized linear translation stage allowing to change the relative position of the patterned surface with respect to the glass lens. Contact pictures are recorded through the transparent PDMS substrate using a zoom objective (Leica APO Z16) and a high resolution CMOS camera (SVS Vistek Exo, 2048 x 2048 pixels 2 , 8 bits). Once illuminated in transmission with an LED light spot, microcontact appear as bright disks (Fig. 2) whose center and area (and radius a i for a given asperity i) can be measured with standard thresholding techniques [5,6]. In the absence of any contact with the glass lens, images of the LED light source result in a small residual bright point at the apex of the microlenses. The radius of these bright spots thus sets to about 8 µm the detection threshold of microcontacts.

As already mentioned in the introduction, our study focuses on adhesive contacts close to equilibrium. For glass-rubber single asperity contacts (including glass-PDMS), it is known that adhesive contact equilibrium is achieved more readily during contact loading than unloading [10,11]. We thus privileged a step-by-step contact loading procedure where the normal load P is increased in discrete steps of 4 mN from 4 mN to 20 mN. At each loading step, it was however observed that the equilibrium adhesive state is not reached instantaneously. Indeed, the measured total contact area slowly increased over a duration that however never exceeded ∼ 1000 s within the investigated range of normal loads. Consequently, in our experiments both the normal force and contact pictures are recorded after a dwell time of 1000 s following each loading step. At the end of the loading procedure, the contact is unloaded at an imposed displacement rate between 2 and 20 µm s -1 from the maximum applied load of 20 mN. For each patterned surface, this procedure is repeated over 36 independent contact realizations by moving laterally the PDMS substrate in 0.5 mm long steps with respect to the fixed glass lens. Within the selected load range, no contact was found to occur in the smooth regions between spherical caps.

III. CONTACT MODEL

Modeling the adhesive contact response was performed by extending to the adhesive case a discrete numerical model developed previously for non-adhesive randomly rough contacts [6]. As in FT's original model, this approach assumes that JKR theory holds at the microcontact length scale. Compared to FT's model, our model also takes into account to first-order elastic interactions arising from both the indenting lens curvature and microcontacts interactions. Following Ciavarella et al. [12], one imposes a displacement which is sensitive to the effect of the spatial distribution of contact pressure in the neighboring asperities. More precisely, for each microasperity contact, a shift of the position of the deformable surface is introduced, which results from the vertical displacement caused by the neighboring ones. Accordingly, the inden-tation depth δ i of the i th microasperity contact writes

δ i = δ 0 i + N j =i α ij f (δ j ) , (2)
where δ 0 i > 0 is the indentation depth in the absence of any elastic coupling between microcontacts, and α ij f (δ j ) are the elements of the interaction matrix. δ 0 i is a purely geometrical term simply given by the difference between the positions of the two undeformed surfaces for the prescribed indentation depth ∆ [6]. The sum in the rhs of eqn (2) represents the interaction term derived from JKR contact theory. We take for this term an asymptotic expansion of the JKR solution for the vertical displacement of the surface, instead of its exact expression. As detailed in the Supplementary Information, elements α ij f (δ j ) of the interaction matrix thus read

[α ij f (δ j )] = - 4 3π 1 r ij f (δ j ) , i = j , (3)
where r ij is the distance from the center of asperity i to that of asperity j and R is the radius of curvature of the microasperities. The function f (δ j ) corresponds to the JKR relationship between load and contact radius a i and writes

f (δ j ) = 3Ra i (δ j )δ j -a 3 i (δ j ) 2R (4)
where a i (δ j ) can be deduced from the numerical inversion of the explicit JKR expression for δ i (a i ). For each value of the standard deviation σ of asperity heights, calculations are repeated over 500 different contact realizations generated from Gaussian sets of asperity heights. Asperities are spatially distributed according to a uniform distribution with a non-overlap constraint. Consistently with the experiments, calculations are carried out by increasing step by step the gap between surfaces until the prescribed maximum load is achieved.

Overall, this numerical model differs from FT's original model by 3 ingredients, (i) non-vanishing elastic interactions, (ii) a sphere-on-plane contact rather than plane-on-plane and (iii) a discrete number of asperities. Note that ingredients (i) and (ii) can be tuned to get a discrete version of FT's model, by simply setting [α ij f (δ j )] = 0 in the simulations and by replacing the curved lens by a flat.

IV. RESULTS AND DISCUSSION

A. Pull-off force and total contact area

We have first considered how integrated quantities, i.e. the pull-off force P c and the total area of contact A depend on both the loading parameters and the topography characteristics of the PDMS substrates. Figure 3 shows as solid lines the relationship between normal load P and vertical displacement d during unloading at imposed displacement rate after a dwell time of 1000 s at P =20 mN. Each curve is an average over 36 contact realizations with the standard deviation indicated by the error bars. For the constant height asperity pattern (σ = 0), the measured average pull-off force is P c = 6.6 ± 1.3 mN. For σ = 5µm and σ = 10µm, the pull-off force is reduced to P c = 3.4 ± 1.6 mN and P c = 0.8 ± 1.0 mN, respectively. As a reference, the contact unloading response of the smooth glass/PDMS contact is also shown in the inset of Fig. 3, yielding a pull-off force of 118 ± 5 mN. It turns out that roughness induces a decrease in the pull-off force by about two orders of magnitude. Such a strong decrease can be discussed in the light of the adhesion parameter θ F T introduced by FT in their model (Eq. 1). Physically, θ F T represents the ratio of the elastic force needed to push one microasperity to a depth σ to the adhesive pull-off force experienced by the same microasperity. In other words, it expresses the competition between compressive forces exerted by the higher asperities and the adhesive force between the lower asperities with respect to the mean nominal interface: when the roughness and hence the adhesion parameter increases, the adhesion falls to small values. Taking into account the elastic and adhesive properties of the PDMS substrate (E = 0.51 MPa and w = 40.8 mJ m -2 ) and the radius of the microasperities, one gets θ F T =13.7 and 37.7 for σ = 5 and 10 µm, respectively. For such high θ F T values, FT predict that the ratio of the pull-off force of the rough contact P c to the one measured with the same spherical surface and a smooth PDMS flat P s c , vanishes (cf Fig. 6 in [3]). Consistently with this prediction, we find experimentally that the ratio P c /P s c is very low (3.10 -2 and 6.10 -3 for σ = 5µm and σ = 10µm, respectively).

Noticeabily, the measured values of the pull-off force exhibit a rate dependence which is indicative of dissipative effects. When the unloading velocity is decreased to v = 2 µm s -1 , the pull-off force is reduced to 2.1 ± 0.3 mN for the random substrate with σ = 5 µm (results not shown). Conversely, increasing v to 200 µm s -1 yields an enhanced pull-off force of 9.3 ± 0.6 mN. Note that some hysteresis, consistent with the occurrence of dissipative effects, is also evidenced when step-by-step loading and constant rate unloading curves are compared (see open symbols versus solid lines on Fig. 3).

This adhesion hysteresis due to viscoelastic loading can be best seen when plotting the total contact area A, obtained by summing up all individual microcontact areas, versus P , for both step-by-step loading and constant rate unloading. This is shown in Fig. 4 with the example of the patterned surface with σ = 5 µm. A strong hysteresis is observed, that is enhanced with the unloading rate, consistently with the assumption of viscoelastic effects.

In the case of single-asperity adhesive contacts with FIG. 4: Total area of contact A vs normal load P during a step-by-step contact loading (black disks) and unloadings at a constant imposed displacement rate of 2 µm s -1 (blue circles) and 20 µm s -1 (red circles). The patterned surface shown here has σ = 5 µm. An adhesion hysteresis is observed between equilibrium loading and out-of equilibrium unloading.

rubbers, dissipative viscoelastic effects have long been reported to enhance adhesion (see e.g. [10,13,14]). Here, they are induced at the length scale of microasperity contacts. For the Sylgard silicone used in our experiments, viscoelastic dissipation may especially arise from the presence of free chains [15] and dangling bonds [16] within the polymer network. Such viscoelastic effects are obviously not taken into account in the FT's theory which is based on an equilibrium description of adhesive contacts. In what follows, we now focus on adhesive equilibrium during contact loading only.

Figure 5 shows A(P ) for all three asperity height distributions. As detailed in Section II, we recall that this data was obtained close to adhesive equilibrium by enforcing a dwell time subsequent to each load increment. For each distribution, data points correspond to 36 contact realizations. For a given load P , increasing σ results in a clear decrease of A. For a given σ, A(P ) also shows some non-linearity which departs from the observed behavior with Hertzian contacts with similar surfaces [6]. To further assess the effects of adhesion on the measured A(P ), contact simulations were carried out using our extended discrete model with either w = 0 or w = 40 mJ m -2 . As shown with the solid lines in Fig. 5, taking into account the elastic and adhesive properties of the PDMS substrate provides an accurate description of the experimental A(P ) within experimental errors. Interestingly, comparison with non-adhesive simulations (see the dotted lines in Fig. 5) shows that adhesion has a significant contribution to the load dependence of the contact area. The latter is indeed enhanced by a factor of roughly two depending on σ. Interestingly, such enhancement of the contact area is observed within a range of θ F T values for which FT's model predicts a vanishing small effect of adhesion on the pull-off force. This experimental observation is consistent with the findings of theoretical adhesive rough contact models (see e.g. [17]) which also predict an enhancement of the contact area with adhesion in the limiting case where no pull-off is can be detected. In addition, it is worthwhile noticing that the effects of adhesion on A(P ) decrease with σ, consistently with FT's model predictions.

Another issue is the contribution of the elastic coupling to the A(P ) relationship. This point is directly addressed in the contact simulations by turning on or off the elastic interaction term α ij f (δ j ). For the sphere-onflat contacts under consideration, both long range elastic interactions due to lens curvature and short-range interactions between neighboring microasperity contacts can be involved. In an attempt to quantify both effects, calculations were carried out either for a sphere-on-flat contact situation or for contacts between nominally flat surfaces. As shown in Fig. 6, the effect of the elastic interactions for sphere-on-flat contacts is clearly to decrease the real contact area at a given load. If only short range interactions were involved, the reverse trend would be expected as the effects of such interactions is to induce a shift of the position of the deformable surface which in turn results in a reduction of actual microasperity indentation depths. Indeed, the fact that simulations for nominally flat surfaces with the interaction term either turned on or off yields the same behavior (Fig. 7) shows that, for the microasperity distributions under consideration, short range elastic coupling between neighboring microcontacts are unimportant. As a consequence, it can be argued that elastic interactions mostly arise from long range effects due to the lens curvature.

B. Microcontact radius and contact pressure distributions

We have then considered local quantities, i.e. the radial distribution of microcontact radii a i , the asperity density η and the contact pressure p. Figure 8 shows the angularly averaged radial profiles of a i , < a i > (r) for the regular array of microasperities (σ = 0) measured at different applied normal loads P . Here, r is the radial coordinate from the center of the apparent contact, taken as the barycenter of all microasperity contacts. In Fig. 8, error bars correspond to the measured standard deviation over 36 different contact realizations. Also shown on Fig. 8 with the solid lines are the predictions of our discrete FT model, which are found to be in good agreement with < a i > (r) measured experimentally. Measured values of the microcontact radii can be compared to the theoretical values a 0 = (6πwR 2 /K) 1/3 and a min = (πwR 2 /6K) 1/3 (K = 4E * /3 being an elastic constant) that correspond, within the framework of JKR's model, to zero load and pull-off at fixed grips respectively. With the elastic and adhesive properties of the considered silicone substrate, a min = 6.58 µm, a value that is close to the detection threshold of our measurements (about 8 µm) and a 0 = 21.7 µm. Both values are marked in Fig. 8 with the dashed dotted and dotted lines, respectively. As can clearly be seen in Fig. 8, asperities lying at the periphery of the contact have microcontact radii that systematically fall between a min and a 0 . Accordingly, a significant portion of the microasperity contacts are under tensile load at the periphery of the apparent contact.

Radial contact pressure profiles p(r) were obtained by summing up local forces P i exerted on all microcontacts located within an annulus of width dr = 0.2 mm and radius r. For each microcontact, the value of P i was determined from the measured contact radius a i using the JKR relationship. Figure 9 shows p(r) for σ = 0 and increasing normal loads P . Such p(r) profiles are amenable to a direct comparison with JKR's predictions for a smooth sphere-on-flat contact. Tentative fits using JKR's model were thus performed taking the adhesion energy w as the only fitting parameter. As shown in the inset of Fig. 9 with the example of P = 8 mN, deviations between theory and experiments are clearly evidenced. Recently, Degrandi-Contraires and coworkers [18] have studied the adhesive contact between a flat PDMS decorated with cylindrical pillars of equal heights and a spherical indenter. In their analysis, they measured the apparent contact area as a function of the normal load. Using JKR's theory, they deduced from their measurements an effective adhesion energy attributed to the inherent roughness of the substrate. However, our results in Fig. 9 clearly show the limits of applying to our rough contact interfaces such an approach solely based on JKR's theory with an effective adhesion energy. In other words, this means that for the patterned surfaces used in our experiments, the separation of length scales between the microasperity contacts and the macroscopic contact region is insufficient to allow for the use of a JKR model based approach with an effective adhesion energy.

Figure 10 shows p(r) for surfaces with increasing σ, with the example of P = 12 mN. Clearly, increasing σ results in a decrease of the magnitude of tensile contact pressures at the contact periphery. For both σ = 5 µm and σ = 10 µm, such tensile contact pressures are barely detected. Average mean contact pressure measurements < p i = P i /πa 2 i > carried out at the scale of individual micro-asperity contacts (inset of Fig. 10) shows that some asperities at the periphery of the apparent contact are still under a slight (< p i >≈ -8 kPa) net tensile stress for σ = 5 µm, while all micro-contacts are under FIG. 9: Radial profiles of contact pressure p(r) for the patterned surface with σ = 0 at P = 4 mN (red curve), P = 8 mN (blue curve), P = 12 mN (green curve), P = 16 mN (light blue), and P = 20 mN (black curve). The dotted line corresponds to p = 0. Inset: attempt to fit the experimental contact pressure data at P = 8 mN to JKR's model with the adhesion energy w as a fitting parameter. The fit yields w = 1 mN.

compressive load for σ = 10 µm.

The dependence of adhesion on the standard devia- tion of asperity height distributions can be further analyzed in the light of measurements of the extension of the apparent contact area. The later can be determined quantitatively from the radial profiles of the microcontact density η(r). Such profiles are reported in Fig. 11 for P = 12 mN and for σ = 0, 5 and 10 µm. The vertical blue and red lines correspond to the calculated contact radius for non-adhesive and adhesive contacts of the spherical probe with the smooth PDMS substrate, respectively. From these profiles, one can determine an apparent contact radius a app by arbitrarily setting a threshold on η. Taking this threshold equal to 5 10 3 m -2 yields a app = 1.9, 2.5 and 2.9 mm for σ = 0, 5 and 10 µm, respectively. For σ = 0, the effective contact radius is clearly lower than the one obtained for a smooth contact thus reflecting the decrease in adhesion induced by roughness. When σ is increased, it turns out that the effective contact radius is increased and can even become close to that of the smooth JKR contact. Such feature can be explained by considering the occurrence of two opposite effects. On one side, the roughness induced reduction in adhesion leads to a decrease in the apparent contact size. On the other side, distributing microasperities along the vertical direction extends the size of the effective contact zone with respect to that achieved with all asperities at the same height (i.e. for σ = 0). New microcontacts are thus added in an annular region around the contact where the glass lens/PDMS separation is comparable to the surface roughness σ. Similar effects have already been reported for non-adhesive contacts [6,19] and the extent of the annular region has been shown to scale as R 5/9 σ 2/3 . Here, a similar effect could be invoked in the case of adhesive contacts to explain the extension of the contact at a given applied load when σ is increased.

We now examine how the microcontact density η de- pends on the local contact pressure p. Figure 12 shows, for all considered contact loads, a log-log plot of η for both σ = 5 and 10 µm as a function of p(r). Also shown in this graph are solid lines that correspond to the theoretical predictions of the original FT's model (see Appendix A). Clearly, FT's model is able to capture satisfactorily the measured η(p) over the investigated pressure range with no need to incorporate any elastic interactions between microasperity contacts. Yet, one can note that FT's model slightly underestimates the measured microcontact densities. A possible explanation for this discrepancy could be that the JKR theory embedded in FT's model neglects the "jumping-on" phenomenon. The latter occurs when the distance between interacting surfaces becomes of the order of magnitude of the range of surface forces. As a result, FT's model should underestimate the microasperity density as well as the total area of contact. As detailed in a recent paper by Greenwood [20], FT's model can be extended to take into account such effects. However, comparison between experiments and contact simulations suggests that such effects are limited for the surfaces considered here.

V. CONCLUSION

In this work, we have investigated the adhesive contact between a smooth glass lens and a silicone surface decorated with prescribed distributions of spherical microasperities differing in their standard deviations of asperity heights. Experimental data were analyzed in the light of a discrete adhesive contact model allowing to account for elastic coupling effects arising from lens curvature and interactions between neighboring asperities. As in the original FT's theory, the model developed in this work postulates a JKR constitutive law at the length scale of individual microasperity contacts. Calculations are carried out under the assumption that the elastic and adhesive properties measured in a single asperity contact with a lens of radius R = 103.7 mm can be transposed to the micrometer scale in order to describe the adhesive contacts between the soft micro-asperities and the glass surface. Our model allows to describe satisfactorily the experimental relationship between contact area and load for the various height distributions under consideration. Noticeably, it also indicates that adhesion increases significantly the contact area even if roughness results in a vanishingly small pull-off force as compared to a smooth contact interface. Such an effect would obviously be relevant to the description of the frictional properties of the rough contact interface which remains to be done. From a measurement of the microcontacts areas, we were able to determine the contact pressure distribution within the rough contact, which preserves one of the main feature of smooth JKR contacts, i.e. the occurrence of tensile stresses at the contact periphery. We also showed that the original FT's model allows for an accurate description of the dependence of the microcontacts density on contact pressure, without any need to take into account elastic coupling between neighboring asperities. Such a conclusion holds for the considered surface topographies with asperities a single length scale (R = 100 µm). However, theoretical calculations with more complex fractal surfaces indicate that the elastic coupling becomes important when surface topography is enriched with asperities at various length scales [21]. This issue could be adressed experimentally using more sophisticated patterned surfaces with hierarchical distributions of microasperities. For instance, as a first step, this could be done by removing the non-overlaping constraint used in the design of our patterned surfaces.

Appendix: FT model

In FT's model, the summits of the spherical asperities have a Gaussian height distribution given by φ(z) = 1 σ √ 2π exp -z 2 /2σ 2 . We consider here non-dimensional heights ζ ≡ z/σ so that the height distribution becomes φ(z) = 1 √ 2π exp -ζ 2 /2 . We follow the normalization suggested by Greenwood [20] and define the normalized load P , indentation depth δ and contact radius a as

P = P Rw ; δ = β δ R ; a = β a R (5)
with

β ≡ E * R w 1/3 (6)
Using this normalization, JKR equations now write δ = a 2 -√ 2πa; P = 4/3a 3 -8πa 3 (7) during loading. Assuming that the asperities do not jump into contact, the average load can be expressed as a function of the normalized mean plane separation h ≡ d/σ as

P (h) = N Z √ 2π ∞ 0 f δ exp -h + Zδ 2 /2 dδ (8)
where Z ≡ Rw 2 /E * 2 σ 3 1/3 = 4/3π 2/3 θ F T is FT's adhesion parameter without the numerical factor (3/4)/π 2/3 . As detailed in Appendix 2 of [20], the above expression can be turned to an explicit integral which can readily be evaluated numerically. Similarly, one can write the microcontact density η as

η(h) = η 0 ∞ h φ(ζ)dζ = -η 0 /2 erf h/ (2) -1 (9)
where η 0 is the microasperity density.

where r ij is the distance separating asperities i and j and where F (δ j ) is the normal indentation load of the j th asperity as obtained from JKR's theory F (δ j ) = 3RKa j (δ j )δ j -Ka 3 j (δ j ) 2R

(3)

where K = 16/9E is an elastic constant and a j (δ j ) is the contact radius of the j th asperity. The latter can be determined by numerically inverting the explicit δ(a) JKR's relationship given by

δ = a 2 R - 8πaw 3K ( 4
)
where w is the adhesion energy. From Eq 3, the interaction term α ij (δ j ) can be rewritten as

(α ij (δ j )) = - 4 3π 1 r ij f (δ j ) , i = j (5)
with f (δ j ) = 3Ra j (δ j )δ j -a 3 j (δ j ) 2R (6) This asymptotic expression of the interaction term is expected to be valid as long as the distance r ij between two neighboring microcontacts is large with respect to their respective contact radii. From a comparison with the exact JKR expression [2] for vertical displacement of the free surface (Fig. 2), we found that the asymptotic expression provides a very accurate value of the surface displacement since r ij > 3a i,j , which is always verified experimentally. Normalized vertical displacement of the surface u z R/a 2 calculated using the exact JKR relationship (red curve) and the asymptotic expression given in Eq. 1 (blue curve) for P/(3πwR) = 0.5 (R is the radius of the sphere, r is the radial coordinate and a is the contact radius).

Numerical implementation

Writing V i = z i -d, with d being the approach distance of the two surfaces at the location of the i th asperity, the contact problem is the solution of the following set of non-linear equations F (δ) = δ -V -A.ϕ (δ) = 0 (7) with A = (α ij (δ j )) and ϕ(u) = f (u) u 0 ; ϕ(u) = 0 u < 0 (8) This set of equations is solved using a Newton-Raphson algorithm as follows

δ n+1 = δ n -J -1 F (δ n ) .F (δ n ) ( 9
)
where J F is the Jacobian defined by

J F = ∂F i ∂δ k = δ ik -A.ϕ ′ (δ k ) (10)
According to Eq. ( 10)

J F (δ) = I -Q (11) with Q = α ij .ϕ ′ (δ k ) (12)
At each computation step, J -1 F is calculated and numerically inverted in order to compute δ n+1 from δ n . The calculation is repeated until |δ n+1 -δ n | is lower than a prescribed convergency criterion.

FIG. 1 :

 1 FIG. 1: Topography of the micromilled PMMA mold. (a) Typical optical profilometry profiles of the spherical cavities along a diameter before (red, lower curve) and after (blue, upper curve) exposure to a saturated CHCl3 vapor. (b) Probability density functions of the height distributions measured on both CHCl3 vapor treated PMMA random molds using optical profilometry. Solid lines are Gaussian fits with σ = 5.13 ± 0.19 µm (red) and σ = 9.22 ± 0.15 µm (blue).

FIG. 2 :

 2 FIG.2: Image in transmission of a typical contact under compressive load (P = 12 mN). Microcontacts appear as bright disks on top of the asperities. The white bar is 1 mm long.

FIG. 3 :

 3 FIG.3: Normal load P as a function of the displacement d (solid lines) during unloading at an imposed displacement rate of 20 µm s -1 , with σ = 0 (red lowest curve), σ = 5 µm (blue middle curve) and σ = 10 µm (green upper curve). Open symbols (with the same color code) correspond to step-bystep contact loading with the same patterned surfaces. Inset: P versus d for the smooth PDMS substrate.

FIG. 5 :FIG. 6 :

 56 FIG.5: Total contact area A as a function of the normal load P with σ = 0 (red disks), σ = 5 µm (blue disks) and σ = 10 µm (green disks). Solid and dotted lines (with the same color code) are the predictions of the adhesive and nonadhesive contact simulations, respectively.

FIG. 7 :

 7 FIG.7: Log-log plot of the calculated total area of contact A normalized by the number of asperities N as a function of the applied load P divided by the nominal apparent contact area Anom for nominally flat surfaces with σ = 5 µm. Calculations are carried out with (solid lines) and without (dotted lines) elastic interactions.

FIG. 8 :

 8 FIG.8: Radial profiles of average microcontact radius < ai >. Red: P = 4 mN, blue: P = 8 mN, green: P = 12 mN, light blue: P = 16 mN, black: P = 20 mN. Solid lines correspond to adhesive contact simulations. The dot-dashed and dashed lines correspond to the theoretical values of pull-off (at fixed grips) and zero load, respectively. Error bars correspond to standard deviations over 36 contact realizations.

FIG. 10 :

 10 FIG.10: Radial profiles of contact pressure p(r) at P =12 mN for σ = 0 (•), σ = 5 µm ( ) and σ = 10 µm (⋄). Inset: corresponding average microcontact pressure < pi > profiles for the same conditions.

FIG. 11 :

 11 FIG. 11: Radial profiles of microcontact density η(r) at P = 12 mN for σ = 0 (•), σ = 5µm ( ) and σ = 10µm (⋄). The vertical blue (resp. red) line marks the theoretical value of the contact radius according to Hertz (resp. JKR)'s theories with E = 0.51 MPa and w = 40.80 mJ m -2 . The horizontal dashed dotted line corresponds to the average microasperity density.

FIG. 12 :

 12 FIG. 12: Log-log plot of the dependence of the average microcontact density η with the local contact pressure p. Experimental data points for σ = 0 (resp. σ = 10 µm) are shown with the red circles (resp. blue squares). Solid lines are predictions of FT's model.

Figure 2 :

 2 Figure2: Normalized vertical displacement of the surface u z R/a 2 calculated using the exact JKR relationship (red curve) and the asymptotic expression given in Eq. 1 (blue curve) for P/(3πwR) = 0.5 (R is the radius of the sphere, r is the radial coordinate and a is the contact radius).

Acknowledgments

The authors wish to thank Guido Violano for his kind help in the obtention of some of the experimental data.

Supplementary Information Derivation and implementation of the multiasperity adhesive contact model 1 Derivation of the model

We consider the indentation of a nominally flat elastic plane decorated with a distribution of spherical caps by a rigid, smooth and spherical probe. The rough surface is described by a collection of N spherical asperities having all the same radius of curvature R and prescribed positions (x i , y i , z i ). Locally, the indentation of each asperity is assumed to obey JKR's theory. Elastic interactions between micro-contacts are accounted for by introducing a shift of the position of the deformable surface seen by each asperity due to the deflection caused by the neighboring ones. Accordingly, the actual indentation δ i of the i th asperity is thus given by

where δ 0 i > 0 is the indentation depth in the absence of any elastic coupling between microcontacts. As shown in Fig. 1,δ 0 i is a purely geometrical term simply given by the difference between the position of the two undeformed surfaces for the prescribed indentation depth ∆. The sum in the rhs of Eq. 1 represents the interaction term derived from JKR's theory. Instead Figure 1: Sketch of the geometric configuration for the indentation of the patterned surface by a smooth rigid lens. ∆ is the prescribed indentation depth taking as a reference for the vertical position of the indenting sphere the altitude at which the smooth surface is touching the uppermost asperity. of providing its exact expression, we derive α ij from an asymptotic expansion that is based on the Green's tensor [1] for a point loading on an elastic incompressible half-plane. It writes