Michał Derezi
email: mderezin@berkeley.edu

Daniele Calandriello
email: daniele.calandriello@iit.it

Michal Valko
email: valkom@deepmind.com

Exact sampling of determinantal point processes with sublinear time preprocessing

We study the complexity of sampling from a distribution over all index subsets of the set {1, ..., n} with the probability of a subset S proportional to the determinant of the submatrix L S of some n × n positive semidefinite matrix L, where L S corresponds to the entries of L indexed by S. Known as a determinantal point process (DPP), this distribution is used in machine learning to induce diversity in subset selection. When sampling from DDPs, we often wish to sample multiple subsets S with small expected size k E [|S|] n from a very large matrix L, so it is important to minimize the preprocessing cost of the procedure (performed once) as well as the sampling cost (performed repeatedly). For this purpose we provide DPP-VFX, a new algorithm which, given access only to L, samples exactly from a determinantal point process while satisfying the following two properties: (1) its preprocessing cost is n • poly(k), i.e., sublinear in the size of L, and (2) its sampling cost is poly(k), i.e., independent of the size of L. Prior to our results, state-of-the-art exact samplers required O(n 3 ) preprocessing time and sampling time linear in n or dependent on the spectral properties of L. We furthermore give a reduction which allows using our algorithm for exact sampling from cardinality constrained determinantal point processes with n • poly(k) time preprocessing. Our implementation of DPP-VFX is provided at https://github.com/guilgautier/DPPy/.

Introduction

Given a positive semi-definite (psd) n × n matrix L, a determinantal point process DPP(L), also known as an L-ensemble, is a distribution over all 2 n index subsets S ⊆ {1, . . . , n} such that Pr(S) det(L S ) det(I + L) , where L S denotes the |S| × |S| submatrix of L with rows and columns indexed by S. Determinantal point processes naturally appear across many scientific domains [START_REF] Macchi | The coincidence approach to stochastic point processes[END_REF][START_REF] Bardenet | On a few statistical applications of determinantal point processes[END_REF][START_REF] Guenoche | Random spanning tree[END_REF], and they have emerged as an important tool in machine learning [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF] for inducing diversity in subset selection and as a variance reduction approach. DPP sampling has been successfully applied in core ML problems such as recommender systems [START_REF] Gillenwater | Maximizing induced cardinality under a determinantal point process[END_REF][START_REF] Chen | Fast greedy map inference for determinantal point process to improve recommendation diversity[END_REF][START_REF] Brunel | Learning signed determinantal point processes through the principal minor assignment problem[END_REF], stochastic optimization [START_REF] Zhang | Determinantal point processes for mini-batch diversification[END_REF][START_REF] Zhang | Active minibatch sampling using repulsive point processes[END_REF], data summarization [CKS + 18], Gaussian processes [START_REF] Mariet | Kronecker determinantal point processes[END_REF][START_REF] Burt | Rates of convergence for sparse variational Gaussian process regression[END_REF], exact DPP k-DPP first sample subsequent samples [HKP + 06, KT11] n 3 nk2 [START_REF] Anari | Monte carlo markov chain algorithms for sampling strongly rayleigh distributions and determinantal point processes[END_REF] x x n • poly(k) n • poly(k) [START_REF] Li | Fast mixing markov chains for strongly rayleigh measures, dpps, and constrained sampling[END_REF] x x n 2 • poly(k) n 2 • poly(k) [START_REF] Launay | Exact Sampling of Determinantal Point Processes without Eigendecomposition[END_REF] x n 3 poly(k • (1 + L )) [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF] x n 3 poly(rank(L)) DPP-VFX (this paper) n • poly(k) poly(k)

Table 1: Comparison of DPP and k-DPP algorithms using the L-ensemble representation. For a DPP, k denotes the expected subset size. Note that k ≤ rank(L) ≤ n. We omit log terms for clarity.

experimental design [START_REF] Dereziński | Unbiased estimates for linear regression via volume sampling[END_REF][START_REF] Dereziński | Leveraged volume sampling for linear regression[END_REF], and many more. In these applications, we often wish to efficiently produce many DPP samples of small expected size 2 k E[|S|] given a large matrix L.

Sometimes, the distribution is restricted to subsets of fixed size |S| = k n, denoted k-DPP(L). [HKP + 06] gave an algorithm for drawing samples from DPP(L) distributed exactly, later adapted to k-DPP(L) by [START_REF] Kulesza | k-DPPs: Fixed-Size Determinantal Point Processes[END_REF], which can be implemented to run in polynomial time. In many applications, however, sampling is still a computational bottleneck because the algorithm requires performing the eigendecomposition of matrix L at the cost of O(n 3 ). In addition to that initial cost, producing many independent samples S 1 , S 2 , . . . at high frequency poses a challenge because the cost of each sample is at least linear in n. Many alternative algorithms exist for both DPPs and k-DPPs to reduce the computational cost of preprocessing and/or sampling, including many approximate and heuristic approaches. Contrary to approximate solutions, we present an algorithm which samples exactly from a DPP or a k-DPP with the initial preprocessing cost sublinear in the size of L and the sampling cost independent of the size of L.

Theorem 1 For a psd n × n matrix L, let S 1 , S 2 , . . . be i.i.d. random sets from DPP(L) or from any k-DPP(L). Then, there is an algorithm which, given access to L, returns a) the first subset, S 1 , in: n • poly(k) polylog(n) time, b) each subsequent S i in:

poly(k) time.

We refer to this algorithm as the Very Fast and eXact DPP sampler, or DPP-VFX. Table 1 compares DPP-VFX with other DPP and k-DPP sampling algorithms. In this comparison, we feature the methods that provide strong accuracy guarantees. As seen from the table, our algorithm is the first exact sampler to achieve sublinear overall runtime. Only the approximate MCMC sampler of [START_REF] Anari | Monte carlo markov chain algorithms for sampling strongly rayleigh distributions and determinantal point processes[END_REF] matches our n • poly(k) complexity (and only for a k-DPP), but for this method every next sample is equally expensive, making it less practical when repeated sampling is needed. In fact, to our knowledge, no other exact or approximate method (with rigorous approximation guarantees) achieves poly(k) sampling time of the present paper.

Our method is based on a technique developed recently by [START_REF] Dereziński | Leveraged volume sampling for linear regression[END_REF][START_REF] Dereziński | Correcting the bias in least squares regression with volume-rescaled sampling[END_REF] and later extended by [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF]. In this approach, we carefully downsample the index set [n] = {1, ..., n} to a sample σ = (σ 1 , ..., σ t ) ∈ [n] t that is small but still sufficiently larger than the expected target size k, and then run a DPP on σ. As the downsampling distribution we use a regularized determinantal point process (R-DPP), proposed by [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF], which (informally) samples σ with probability Pr(σ)

∼ det(I + L σ ),
where L is a rescaled version of L. We can summarize this approach as follows, where

|S| ≤ t n, {1, ..., n} σ∼R-DPP -----→ (σ 1 , ..., σ t ) S∼DPP ----→ S = {σ i : i ∈ S}.
The DPP algorithm proposed by [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF] follows the same diagram, however it requires that the size of the intermediate sample σ be Ω(rank(L) • k). This means that their method provides improvement over [HKP + 06] only when L can be decomposed as XX for some n × r matrix X, with r n. However, in practice, matrix L is often only approximately low-rank, i.e., it exhibits some form of eigenvalue decay but it does not have a low-rank factorization. In this case, the results of [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF] are vacuous both in terms of the preprocessing cost and the sampling cost, in that obtaining every sample would take Ω(n 3 ). We propose a different R-DPP implementation (see DPP-VFX as Algorithm 1)

where the expected size of σ is O(k 2 ). To make the algorithm efficient, we use new connections between determinantal point processes, ridge leverage scores, and Nyström approximation.

Definition 1 Given a psd matrix L, its ith λ-ridge leverage score (RLS) τ i (λ) is the ith diagonal entry of L(λI + L) -1 . The λ-effective dimension d eff (λ) is the sum of the leverage scores, i τ i (λ).

An important connection between RLSs and DPPs is that when S ∼ DPP(L), the marginal probability of index i being sampled into S is equal to the ith 1-ridge leverage score of L, and the expected size k of S is equal to the 1-effective dimension,

Pr(i ∈ S) = L(I + L) -1 ii = τ i (1), k E |S| = tr L(I + L) -1 = d eff (1).
Intuitively, if the marginal probability of i is high, then this index should likely make it into the intermediate sample σ. This suggests that i.i.d. sampling of the indices σ 1 , ..., σ t proportionally to 1-ridge leverage scores, i.e., Pr(σ 1 = i) ∝ τ i (1), should serve as a reasonable and cheap heuristic for constructing σ. In fact, we can show that this distribution can be easily corrected by rejection sampling to become the R-DPP that we need. Computing ridge leverage scores exactly costs O(n 3 ), so instead we compute them approximately by first constructing a Nyström approximation of L.

Definition 2 Let L be a psd matrix and C a subset of its row/column indices with size m |C|. Then we define the Nyström approximation of L based on C as the n

× n matrix L (L C,I ) L + C L C,I .
Here, L C,I denotes an m × n matrix consisting of (entire) rows of L indexed by C and (•) + denotes the Moore-Penrose pseudoinverse. Since we use rejection sampling to achieve the right intermediate distribution, the correctness of our algorithm does not depend on which Nyström approximation is chosen. However, the subset C greatly influences the computational cost of the sampling through the rank of L and the probability of rejecting a sample. Since rank( L) = |C|, operations such as multiplication and inversion involving the Nyström approximation will scale with m, and therefore a small subset increases efficiency. However, if L is too different from L, the probability of rejecting the sample will be very high and the algorithm inefficient. In this case, a slightly larger subset could improve accuracy and acceptance rate without increasing too much the cost of handling L. Therefore, subset C has to be selected so that it is both small and accurately represents the matrix L. Here, we once again rely on ridge leverage score sampling which has been effectively used for obtaining good Nyström approximations in a number of prior works such as [AM15, CLV17, RCCR18].

While our main algorithm can sample only from the random-size DPP, and not from the fixed-size k-DPP, we present a rigorous reduction argument which lets us use our DPP algorithm to sample exactly from a k-DPP (for any k) with a small computational overhead.

Related work Prior to our work, fast exact sampling from generic DPPs has been considered out of reach. The first procedure to sample general DPPs was given by [HKP + 06] and even most recent exact refinements [START_REF] Launay | Exact Sampling of Determinantal Point Processes without Eigendecomposition[END_REF][START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF][START_REF] Poulson | High-performance sampling of generic determinantal point processes[END_REF], when the DPP is represented in the form of an L-ensemble, require preprocessing that amounts to an expensive n × n matrix diagonalization at the cost O(n 3 ), which is shown as the first-sample complexity column in Table 1.

Nonetheless, there are well-known samplers for very specific DPPs that are both fast and exact, for instance for sampling uniform spanning trees [Ald90, Bro89, PW98], which leaves the possibility of a more generic fast sampler open. Since the sampling from DPPs has several practical large scale machine learning applications [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF], there are now a number of methods known to be able to sample from a DPP approximately, outlined in the following paragraphs.

As DPPs can be specified by kernels (L-kernels or K-kernels), a natural approximation strategy is to resort to low-rank approximations [KT11, GKT12, AKFT13, LJS16a]. For example, [START_REF] Hafiz Affandi | Nystrom approximation for large-scale determinantal processes[END_REF] provides approximate guarantee for the probability of any subset being sampled as a function of eigengaps of the L-kernel. Next, [START_REF] Li | Efficient sampling for k-determinantal point processes[END_REF] construct coresets approximating a given k-DPP and then use them for sampling. In their Section 4.1, [START_REF] Li | Efficient sampling for k-determinantal point processes[END_REF] show in which cases we can hope for a good approximation. These guarantees become tight if these approximations (Nyström subspace, coresets) are aligned with data. In our work, we aim for an adaptive approach that is able to provide a good approximation for any DPP.

The second class of approaches are based on Markov chain Monte-Carlo [START_REF] Metropolis | The Monte Carlo method[END_REF] techniques [Kan13, RK15, AGR16, LJS16b, GBV17]. There are known polynomial bounds on the mixing rates [START_REF] Diaconis | Geometric Bounds for Eigenvalues of Markov Chains[END_REF] of MCMC chains with arbitrary DPPs as their limiting measure. In particular, [START_REF] Anari | Monte carlo markov chain algorithms for sampling strongly rayleigh distributions and determinantal point processes[END_REF] showed them for cardinality-constrained DPPs and [START_REF] Li | Fast mixing markov chains for strongly rayleigh measures, dpps, and constrained sampling[END_REF] for the general case. The two chains have mixing times which are, respectively, linear and quadratic in n (see Table 1). Unfortunately, for any subsequent sample we need to wait until the chain mixes again.

Neither the known low-rank approximations or the known MCMC methods are able to provide samples that are exactly distributed according to a DPP (also called perfect sampling). This is not surprising as having scalable and exact sampling is very challenging in general. For example, methods based on rejection sampling are always exact, but they typically scale poorly to high-dimensional data and are adversely affected by the spikes in the distribution [START_REF] Erraqabi | Pliable rejection sampling[END_REF], resulting in high rejection rate and inefficiency. Surprisingly, our method is based on both low-rank approximation (a source of inaccuracy) and rejection sampling (a common source of inefficiency). In the following section, we show how to obtain a perfect DPP sampler from a Nyström approximation of the L-kernel. Then, to guarantee efficiency, in Section 3 we bound the number of rejections, which is possible thanks to the use of intermediate downsampling. Finally, we also allow C, D to be multisets or sequences, in which case each row/column is duplicated in the matrix according to its multiplicity (and in the case of sequences, we order the rows/columns as they appear in the sequence). Note that with this notation if

L = BB then L C,D = B C,I B D,I .
As discussed in the introduction, our method relies on an intermediate downsampling distribution to reduce the size of the problem. The exactness of our sampler relies on the careful choice of that intermediate distribution. To that end, we use regularized determinantal processes, introduced by [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF]. In the below definition, we adapt them to the kernel setting.

Definition 3 Given an n × n psd matrix L, distribution p (p 1 , . . . , p n ) and r > 0, let L denote an

n × n matrix such that L ij 1 r √ pipj L ij for all i, j ∈ [n]. We define R-DPP r p (L) as a distribution over events A ⊆ ∞ k=0 [n] k , where Pr(A) E σ 1 [σ∈A] det(I + L σ ) det(I + L) , for σ = (σ 1 , . . . , σ t ) i.i.d.
∼ p, t ∼ Poisson(r).

Since the term det(I + L σ ) has the same form as the normalization constant of DPP( L σ ), an easy calculation shows that the R-DPP can be used as an intermediate distribution in our algorithm without introducing any distortion in the sampling.

Proposition 1 (Der19, Theorem 8) For any L, p, r, and L defined as in Definition 3,

if σ ∼ R-DPP r p (L) and S ∼ DPP( L σ ) then {σ i : i ∈ S} ∼ DPP(L).
To sample from the R-DPP, our algorithm uses rejection sampling, where the proposal distribution is sampling i.i.d. proportionally to the approximate 1-ridge leverage scores l i ≈ τ i (1) (see Definition 1 and the following discussion), computed using any Nyström approximation L of matrix L. Apart from L, the algorithm also requires an additional parameter q, which controls the size of the intermediate sample. Because of rejection sampling and Proposition 1, the correctness of the algorithm does not depend on the choice of L and q, as demonstrated in the following result. The key part of the proof involves showing that the acceptance probability in Line 6 is bounded by 1. Here, we obtain a considerably tighter bound than the one achieved by [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF], which allows us to use a much smaller intermediate sample σ (see Section 3) while maintaning the efficiency of rejection sampling.

Theorem 2 Given a psd matrix L, any one of its Nyström approximations L and any positive q, DPP-VFX (Algorithm 1) returns S ∼ DPP(L).

An important implication of Theorem 2 is that even though the choice of L affects the overall execution of the algorithm, it does not affect the distribution of the output. Therefore we can reuse Algorithm 1 DPP-VFX sampling S ∼ DPP(L)

Input: L ∈ R n×n , its Nyström approximation L, q > 0 1: Compute l i = (L -L) + L(I + L) -1 ii ≈ Pr(i ∈ S) 2: Initialize s = i l i , z = tr L(I + L) -1 , L = s q 1 √ lilj L ij ij 3: repeat 4: sample t ∼ Poisson(q e s/q ) 5: sample σ 1 , . . . , σ t i.i.d. ∼ ( l1 s , • • • , ln s ), 6:
sample Acc ∼ Bernoulli e z det(I+ Lσ) e ts/q det(I+ L)

7: until Acc = true 8: sample S ∼ DPP L σ 9: return S = {σ i : i ∈ S} the same L to produce multiple independent samples S 1 , S 2 , ... ∼ DPP(L), which we prove in Appendix A. This is particularly important because, as we will see in Section 3, DPP-VFX can generate successive samples S 2 , S 3 , . . . much more cheaply than the first sample S 1 .

Lemma 1 Let C ⊆ [n] be a random set variable with any distribution. Suppose that S 1 and S 2 are returned by two executions of DPP-VFX, both using inputs constructed from the same L and

L = L I,C L + C L C,I .
Then S 1 and S 2 are (unconditionally) independent.

Before proceeding with the proof, we highlight for clarity that the Poisson r.v. t in DPP-VFX (Algorithm 1, Line 4) has a different role than the the Poisson r.v. used in the definition of R-DPPs. As a result, the definition of L (Line 2) is not exactly the one that would be suggested by Definition 3.

In particular, DPP-VFX uses q e s/q instead of just q as the mean parameter for t. The extra e s/q factor is necessary to correct for the rejection sampling step in Line 6 during the sampling loop, see the proof below.

Proof of Theorem 2

We start by showing that the Bernoulli probability in Line 6 is bounded by 1. Note that this is important not only to sample correctly, but also when we later establish the efficiency of the algorithm. If we showed a weaker upper bound, say c > 1, we could always divide the expression by c and retain the correctness, however it would also be c times less likely that Acc = 1.

Since L is a Nyström approximation for some C ⊆ I, it can be written as

L = L I,C L + C L C,I = BB C,I L + C B C,I B = BPB,
= exp t i=1 s ql σi B(I + PB BP) -1 B σiσi • e -z = e ts/q e -z ,
where the last equality follows because

B(I + PB BP) -1 B = B I -PB (I + BPB ) -1 BP B = L -L(I + L) -1 L = L -L + L(I + L) -1 .
Thus, we showed that the expression in Line 6 is valid. Let σ denote the random variable distributed as σ is after exiting the repeat loop. It follows that

Pr( σ ∈ A) ∝ E σ 1 [σ∈A] e z det(I + L σ ) e ts/q det(I + L) ∝ ∞ t=0
(q e s/q ) t e q e s/q t! • e -ts/q E σ 1

[σ∈A] det(I + L σ ) | t ∝ E t E σ 1 [σ∈A] det(I + L σ ) | t = t for t ∼ Poisson(q),
which shows that σ ∼ R-DPP q l (L) for l = ( l1 s , • • • , ln s ). The claim follows from Proposition 1.

Conditions for fast sampling

The complexity cost of DPP-VFX can be roughly summarized as follows: we pay a large onetime cost to precompute L and all its associated quantities, and then we pay a smaller cost in the rejection sampling scheme which must be multiplied by the number of times we repeat the loop until acceptance. We first show that if L is sufficiently close to L then we will exit the loop with high probability. We then analyze how accurate the precomputing step needs to be to satisfy this condition. Finally, we bound the overall computational cost, which includes also the final step where we sample S out of the intermediate sample σ using any off-the-shelf exact DPP sampler.

1) Bounding the number of rejections

The following result presents the two conditions needed for achieving efficient rejection sampling in DPP-VFX. First, the Nyström approximation needs to be accurate enough, and second, the intermediate sample size (controlled by parameter q) needs to be Ω(k 2 ). This is a significant improvement over the guarantee of [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF] where the intermediate sample size is Ω(rank(L) • k), which is only meaningful for low-rank kernels. The main novelty in this proof comes from showing the following lower bound on the ratio of the determinants of I + L and I + L, when L is a Nyström approximation: det(I + L)/ det(I + L) ≥ e k-z , where z tr( L(I + L) -1 ). Remarkably, this bound exploits the fact that any Nyström approximation of L = BB can be written as L = BPB , where P is a projection matrix. Note that while our result holds in the worst case, in deployement, the conditions on L and on q can be considerably relaxed.

Theorem 3 If the Nyström approximation L and the intermediate sample size parameter q satisfy tr L(I + L) -1 L -L(I + L) -1 L ≤ 1 and q ≥ max{s 2 , s}, then Pr(Acc = true) ≥ e -2 . Therefore, with probability 1 -δ, Algorithm 1 exits the rejection sampling loop after at most O(log δ -1 ) iterations and, after precomputing all of the inputs, the time complexity of the rejection sampling loop is O k 6 log δ -1 + log 4 δ -1 .

Proof We first bound the number of rejections, and then discuss the complexity of each loop iteration. Let σ be distributed as in Line 5. The probability of exiting the repeat loop at each iteration is

P E σ e z det(I + L σ ) e ts/q det(I + L) = ∞ t=0
(q e s/q ) t e q e s/q t! • e z-ts/q det(I + L)

E σ det(I + L σ ) | t = e q-q e s/q +z det(I + L) ∞ t=0
q t e q t! E σ det(I + L σ ) | t = e q-q e s/q +z det(I + L) det(I + L)

, where the last equality follows because the infinite series computes the normalization constant of R-DPP q l (L) given in Definition 3. If s ≥ 1 then q = s 2 and the inequality e x ≤ 1 + x + x 2 for x ∈ [0, 1] implies that q -q e s/q + z = s 2 (1 -e 1/s + 1/s) + z -s ≥ -1 + z -s. On the other hand, if q = s ∈ [0, 1], then q -qe s/q + z = (1 -e)s + z ≥ -1 + z -s. We proceed to lower bound the determinantal ratio. Here, let L = BB and L = BPB, where P is a projection matrix. Then, Putting all together, e q-q e s/q +z det(I + L)

det(I + L) ≥ exp -1 + z -s + k -tr( L) + tr L(I + L) -1 L ,
and using the definitions k = tr L(I+L) -1 , s = tr(L-L+ L( L+I) -1 ), and z = tr( L( L+I) -1 ) on the middle term z -s + k -tr( L) we have

tr L( L + I) -1 + L(I + L) -1 -L + L -L( L + I) -1 -L = tr L(I + L) -1 -L = tr L(I + L) -1 -L(I + L) -1 (I + L) = -tr L(I + L) -1 L ,
and therefore, e q-q e s/q +z det(I + L)

det(I + L) ≥ exp -1 + tr L(I + L) -1 L -tr L(I + L) -1 L ,
and we obtain our condition. Thus, with probability 1 -δ, the main loop will be repeated O(log δ -1 ) times. We now quantify the cost of a single loop iteration. First note that since the number of samples t i drawn from l in the ith iteration of the loop is a Poisson distributed random variable, a standard Poisson tail bound implies that with probability 1 -δ, all iterations will satisfy t i = O(k 2 + log δ -1 ).

Then, drawing the Poisson r.v. t i (Line 4) requires O(k 2 + log δ -1 ) time. Drawing σ from the multinomial ( l1 s , • • • , ln s ) can be done by first sorting the l i , which is a one-time O(n log(n)) cost, and then using specialized samplers which require O(1) time [START_REF] Bringmann | Efficient sampling methods for discrete distributions[END_REF]. Therefore, the dominant cost is computing the determinant of the matrix I + L σ in O(t 3 i ), and the result follows. 2) Bounding the precompute cost All that is left is to control the cost of the precomputation phase. We separate the analysis into two steps: how much it costs to choose L to satisfy the assumption of Theorem 3, and how much it costs to compute everything else given L, see Appendix A.

Lemma 2 Let L be constructed by sampling m = O(k 3 log n δ ) columns proportionally to their RLS. Then, with probability 1 -δ, L satisfies the assumption of Theorem 3.

There exist many algorithms to sample columns proportionally to their RLSs. For example, we can take the BLESS [START_REF] Rudi | On fast leverage score sampling and optimal learning[END_REF] with the following guarantee.

Proposition 2 (RCCR18, Theorem 1) There exists an algorithm that with probability 1 -δ samples m columns proportionally to their RLSs in O(nk 2 log 2 n δ + k 3 log 4 n δ + m) time. We can now compute the remaining preprocessing costs, given a Nyström approximation L.

Lemma 3 Given L with rank m, we can compute l i , s, z, and L in O(nm 2 + m 3 ) time.

3) Bounding the overall cost We can now fully characterize the computational cost.

Theorem 1 (restated for DPPs only) For a psd n × n matrix L, let S 1 , S 2 be i.i.d. random sets from DPP(L). Denote with L a Nyström approximation of L obtained by sampling m = O(k 3 log(n/δ)) of its columns proportionally to their RLSs. If q ≥ max{s 2 , s}, then w.p. 1 -δ, DPP-VFX returns a) subset S 1 in: O(nk

6 log 2 n δ + k 9 log 3 n δ + k 3 log 4 n δ ) time, b) then, S 2 in: O k 6 log 1 δ + log 4 1 δ time.
Discussion The above result follows from the bounds on the precompute costs, on the number of iteratios, on the iterations cost, and on the fact that the final DPP sampling step S ∼ DPP L σ requires O(t 3 ) ≤ O((k 2 +log δ -1 ) 3 ) time. Note however that due to the nature of rejection sampling, as long as we exit the loop, i.e., we accept the sample, the output of DPP-VFX is guaranteed to follow the DPP distribution for any value of m and q. In Theorem 1 we set m = (k 3 log n δ ) and q ≥ max{s 2 , s} to satisfy Theorem 3 and guarantee a constant acceptance probability in the rejection sampling loop, but this might not be necessary or even desirable in practice. Experimentally, much smaller values of m, starting from m = Ω(k log n δ ) seem to be sufficient to accept the sample, while at the same time a smaller m greatly reduces the preprocessing costs. In general, we recommend to separate DPP-VFX in three phases. First, compute an accurate estimate of the RLS using off-theshelf algorithms in O(nk 2 log 2 n δ + k 3 log 4 n δ ) time. Then, sample a small number m of columns to construct an explorative L, and try to run DPP-VFX. If the rejection sampling loop does not terminate sufficiently fast, then we can reuse the RLS estimates to compute a more accurate L for a larger m. Using a simple doubling schedule for m, this strategy quickly reaches a regime where DPP-VFX is w.h.p. guaranteed to accept, resulting in faster sampling.

Reduction from DPPs to k-DPPs

We next show that with a simple extra rejection sampling step we can efficiently transform any exact DPP sampler into an exact k-DPP sampler.

A common heuristic to sample S from a k-DPP is to first sample S from a DPP, and then reject the sample if the size of S is not exactly k. As we show in this section, the success probability of this procedure can be improved by appropriately rescaling L by a constant factor α,

sample S α ∼ DPP(αL), accept if |S α | = k.
Note that rescaling the DPP by a constant α as above only changes the expected size of the set S α , and not its distribution. Therefore, if we accept only sets with size k, we will be sampling exactly from our k-DPP. Moreover, if k α = E[|S α |] is close to k, the success probability will improve. With a slight abuse of notation, in the context of k-DPPs we will indicate with k the desired size of S α , i.e., the final output size at acceptance, and with k α = E[|S α |] the expected size of the scaled DPP.

While the above rejection sampling heuristic is widespread, until now there has been no proof that this heuristic can provably succeed in few rejections We solve this open question with two new results. First, we show that for an appropriate rescaling α we only reject S α roughly O( √ k) times. Then, we show how to find such an α with a O(n • poly(k)) time preprocessing step.

Theorem 4 There exists constant C > 0 such that for any rank n psd matrix L and k ∈ [n], there is α > 0 with the following property: if we sample

S α ∼ DPP(α L), then Pr(|S α | = k) ≥ 1 C √ k •
The proof in Appendix B. relies on a known Chernoff bound for the sample size |S α | of a DPP. When applied naïvely, the inequality does not offer a lower bound on the probability of any single sample size. However we show that the probability mass is concentrated on O( √ k α ) sizes. This leads to a lower bound on the sample size with the largest probability, i.e., the mode of the distribution. Then, it remains to observe that for any k ∈ [n] we can always find α for which k is the mode of |S α |. We conclude that given α , the rejection sampling scheme described above transforms any poly(k) time DPP sampler into a poly(k) time k-DPP sampler. It remains to efficiently find α , which once again relies on using a Nyström approximation of L.

Lemma 4 If k ≥ 1, then there is an algorithm that finds α in O(n • poly(k)) time.

While the existence proof of α is general and based on simple unimodality, this characterization is not sufficient to control the mode of the DPP when α is perturbed, as it happens during an approximate optimization of α. However, for DPPs the mode can be expressed as a Poisson binomial random variable [H + 56] based on the spectrum of L, which can be controlled [D + 64] using the approximate spectrum of L. While our analysis shows that in the worst case a much more accurate L is necessary compared to Theorem 3, in practice the same m = Ω(k log n δ ) seems to suffice.

Experiments

In this section, we experimentally evaluate3 the performance of DPP-VFX and compare it to exact sampling [HKP + 06] and MCMC-based approaches [START_REF] Anari | Monte carlo markov chain algorithms for sampling strongly rayleigh distributions and determinantal point processes[END_REF]. In particular, we are only interested in evaluating computational performance, i.e., how DPP-VFX and baselines scale with the size n of the matrix L. Note that an inexact DPP sampler, e.g., MCMC samplers or previous approximate rejection samplers [START_REF] Dereziński | Fast determinantal point processes via distortion-free intermediate sampling[END_REF], also had another metric to validate, i.e., they had to empirically show that their samples were close enough to a DPP distribution. However, Section 2 proves that DPP-VFX's output is distributed exactly according to the DPP and therefore it is strictly equivalent to any other exact DPP sampler, e.g., the Gram-Schmidt [START_REF] Jennifer | Approximate inference for determinantal point processes[END_REF] or dual [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF] samplers implemented in DPPy.

We use the infinite MNIST digits dataset (i.e., MNIST8M [START_REF] Loosli | Training invariant support vector machines using selective sampling[END_REF]) as input data, where n varies up to 10 6 and d = 784. All algorithms, including DPP-VFX, are implemented in python as part of the DPPy library [START_REF] Gautier | DPPy: Sampling determinantal point processes with Python[END_REF]. All experiments are carried out on a 24-core CPU and fully take advantage of potential parallelization. For each experiment we report the mean and a 95% confidence interval over 10 runs of each algorithm's runtime. This first experiment compares the time required to generate a first sample. We consider subsets of MNIST8M that go from n = 10 3 to n = 7 • 104 , i.e., the original MNIST dataset, and use an RBF kernel with σ = √ 3d to construct L. For the Nyström approximation we set m = 10d eff (1) ≈ 10k. While this is much lower than the O(k 3 ) value suggested by the theory, as we will see it is already accurate enough to result in drastic runtime improvements over exact and MCMC. Following the strategy of Section 4, for each algorithm we control 4 the size of the output set by rescaling the input matrix L by a constant α such that E[|S α |] = 10. Results are reported in Figure 1.

Exact sampling is performed using eigendecomposition and DPPy's default Gram-Schmidt [START_REF] Jennifer | Approximate inference for determinantal point processes[END_REF] sampler. It is clearly cubic in n, and we could not push it beyond n = 1.5 • 10 4 . For MCMC, we enforce mixing by runnning the chain for nk steps, the minimum recommended by [START_REF] Anari | Monte carlo markov chain algorithms for sampling strongly rayleigh distributions and determinantal point processes[END_REF]. However, for n = 7 • 10 5 the MCMC runtime is 160s and exceeds the plot's limit, while DPP-VFX completes in 16s, an order of magnitude faster. Moreover, DPP-VFX rarely rejects more than 10 times, and the mode of the rejections up to n = 7 • 10 5 is 1, i.e., we mostly accept at the first iteration.

Since in Figure 1, the resampling time of both exact and DPP-VFX are negligible, we investigate the cost of a second sample, i.e., of resampling, on a much larger subset of MNIST8M up to n = 10 6 , normalized to have maximum row norm equal to 1. In this case it is not possible to perform an eigendecomposition of L, so we replace RBFs with a linear kernel, for which the input X represents a factorization L = XX and we use an exact dual sampler on X [START_REF] Kulesza | Determinantal Point Processes for Machine Learning[END_REF]. However, as Figure 2 shows, even with a special algorithm for this simpler setting the resampling process still scales with n. On the other hand, DPP-VFX's complexity (after preprocessing) is light years better as it scales only with k and remains constant regardless of n. ∈ R n×m such that L = BB . We also introduce L m B B to act as a R m×m counterpart to L. Denote with e i the i-th indicator vector. Then, exploiting the fact that BB (I + BB ) -1 = B(I + B B) -1 B for any matrix, we can compute l i as

l i = [L -L + BB (I + BB ) -1 ] ii = [L -L] ii + (I + L m ) -1/2 B e i 2 2 .
Computationally, this means that we first need to compute B, which takes O(m 3 ) time to compute L +/2 C,C , and O(nm 2 ) time for the matrix multiplication. Then, [ L] ii is the 2 norm of the i-th row of B which can be computed in nm time. Similarly, (I + L m ) -1/2 B e i 2 2 requires O(m 3 + nm 2 ) time. To compute s, we simply sum l i , while to compute z we first compute the eigenvalues of L m , a i = λ( L m ) i , in O(m 3 ) time, and then compute z = i a i /(a i + 1). We can also recycle the eigenvalues to precompute log det(I + L) = log det(I + L m ) = i log(a i + 1).

B Omitted proofs for the reduction to k-DPPs

In this section we present the proofs omitted from Section 4. Recall that our approach is based on the following rejection sampling strategy,

sample S α ∼ DPP(αL), accept if |S α | = k.
First, we show the existence of the factor α for which the rejection sampling is efficient.

Theorem 5 (restated Theorem 4) There exists constant C > 0 such that for any rank n PSD matrix L and k ∈ [n], there is α > 0 with the following property: if we sample S α ∼ DPP(α L), then

Pr(|S α | = k) ≥ 1 C √ k • (
| = M α ) = Pr(|S α | = k)
is large enough. In other words, rescaling DPP(αL) to make sure that its mean k α is close to k is sufficient to guarantee that |S α | = k with high enough probability.

Our starting point is a standard Chernoff bound for |S α |.

Proposition 4 (PP14) Given any PSD matrix L, if S α ∼ DPP(αL), then for any a > 0, we have

Pr |S α | -E[|S α |] ≥ a ≤ 5 exp - a 2 16(a + 2 E[|S α |]) •
Note that Proposition 4 is not sufficiently strong by itself, i.e., if we tried to bound the distance

||S α | -E[|S α |]
| to be smaller than 1 we would get vacuous bounds. However, Proposition 4 implies that there is a constant

C > 0 independent of L such that Pr ||S α | -k α | ≥ C √ k α + 1 ≤ 1 2 for all α > 0. In particular, this means that the mode of |S α |, i.e. M α = argmax i Pr(|S α | = i) satisfies Pr(|S α | = M α ) ≥ 1 2C √ k α C √ kα i=-C √ kα Pr(|S α | = k α + i) ≥ 1 4C √ k α • (2)
The distribution of |S α | is given by Pr(|S α | = i) ∝ e i (αL), where e i (•) is the ith elementary symmetric polynomial of the eigenvalues of a matrix. Denoting λ 1 , . . . , λ n as the eigenvalues of L, we can express the elementary symmetric polynomials as the coefficients of the following univariate polynomial with real non-positive roots,

n i=1 (x + αλ i ) = n k=0
x k e n-k (αL).

The non-negative coefficients of such a real-rooted polynomial form a unimodal sequence (Lemma 1.1 in [START_REF] Petter | Unimodality, log-concavity, real-rootedness and beyond[END_REF]), i.e., e 0 (αL) ≤ • • • ≤ e Mα (αL) ≥ • • • ≥ e n (αL), with the mode (shared between no more than two positions k, k+1) being close to the mean k α : |M α -k α | ≤ 1 (Theorem 2.2 in [START_REF] Petter | Unimodality, log-concavity, real-rootedness and beyond[END_REF]). Moreover, it is easy to see that M 0 = 0 and M α = n for large enough α, so since the sequence is continuous w.r.t. α, for every k ∈ [n] there is an α such that Pr(|S α | = k) = Pr(|S α | = M α ) (every k can become one of the modes). In light of (2), this means that

Pr(|S α | = k) ≥ 1 4C √ k α ≥ 1 4C √ k + 1
, where the last inequality holds because |k -k α | ≤ 1.

Finally, we show how to find α efficiently.

Lemma 8 If k ≥ 1 there is an algorithm that finds α in O(n • poly(k)) time.

Proof In order to leverage Theorem 4, we need to find an α such that k = M α , that is such that the mode of DPP(α L) is equal to k. Unfortunately simple unimodality is not sufficient to control M α when α is perturbed, as it happens during an approximate optimization of α. We now characterize more in detail the distribution of |S α |.

In

M α =      k if k ≤ k α < k + 1 k+2 , k or k + 1 if k + 1 k+2 ≤ k α ≤ k + 1 -1 n-k+1 , k + 1 if k + 1 -1 n-k+1 < k α ≤ k + 1.
Therefore it is sufficient to find any constant α that places k α in the interval [k, k + 1 k+2 ). Unfortunately, while the formula for k α = n i=1 λ i (L)/(1/α+λ i (L)) is a unimodal function of the eigenvalues of L which is easy to optimize, the eigenvalues themselves are still very expensive to compute. For efficiency, we can optimize it instead on the eigenvalues of a Nyström approximation L, but we have to be careful to control the error. In particular, remember that k α = E[|S α |] when S ∼ DPP(αL), so given a Nyström approximation L we can define s α tr(α(L -L) + L( L + I/α) -1 ) as a quantity analogous to s from DPP-VFX. Then, we can strengthen Lemma 2 as follows.

Lemma 9 (see also Lemma 2) Let L be constructed by sampling m = O((k α /ε 2 ) log(n/δ)) columns proportionally to their RLS. Then with probability 1 -δ

1 1 + ε k α ≤ s α ≤ 1 1 -ε k α .
Proof of Lemma 9 We simply apply the same reasoning of Lemma 2 on both sides.

Let (1 -ε)s α = k, with ε that we tune shortly. Then proving the first inequality to satisfy Proposition 5 is straightforward, k = (1 -ε)s α ≤ k α . To satisfy the other side we upper bound k α ≤ (1 + ε)s α = (1 -ε)s α + 2εs α . We must now choose ε such that 2εs α = 1/(k + 3) < 1/(k + 2). Substituting, we obtain ε = 1 2(k+3)s α • Plugging this in the definition of s α we obtain that α must be optimized to satisfy

s α = 2k 2 + 6k + 1 2k + 6
, which we plug in the definition of ε obtaining our neccessary accuracy ε = 1/(2k 2 + 6k + 1). Therefore, sampling m = O(k α k 4 ) columns gives us a s α sufficiently accurate to be optimized. However, we still need to bound k α , which we can do as follows using Lemma 9 and k ≥ 1, Finally, note that the bounds on the accuracy of L are extremely conservative. In deployment, it is much faster to try to optimize α on a much coarser L first, e.g., for m = O(k 1 ), and only if this approach fails, then to increase the accuracy of L.

k α ≤ 1 + 1 2k 2 + 6k + 1 s α ≤ 1 + 1 9 s α =

2

  Exact sampling using any Nyström approximation Notation We use [n] to denote the set {1, . . . , n}. For a matrix B ∈ R n×m and index sets C, D, we use B C,D to denote the submatrix of B consisting of the intersection of rows indexed by C with columns indexed by D. If C = D, we use a shorthand B C and if D = [m], we may write B C,I .

  for any B such that L = BB, where P B C,I L + C B C,I is a projection, so that P 2 = P. Let L B B, where the ith row of B is the rescaled ith row of B, i.e. b i s qli b i . Then, we have det(I + L σ ) det(I + L) = det(I + B σ,I B σ,I ) det(I + PB BP) = det (I + B σ,I B σ,I )(I + PB BP) -1 = det I + ( B σ,I B σ,I -PB BP)(I + PB BP) -1 ≤ exp tr ( B σ,I B σ,I -PB BP)(I + PB BP) -1

  det(I + L) det(I + L) = det(I + B B) det(I + PB BP) = det I -(B B -PB BP)(I + B B) -1 -1 ≥ exp tr (B B -PB BP)(I + B B) -1 = exp tr B(I + B B) -1 B -tr BP(I + B B) -1 PB , where we can reformulate k = tr B(I + B B) -1 B and the other elements in terms of L and L as exp k -tr BP(I + B B -B B)(I + B B) -1 PB = exp k -tr( L) + tr BPB B(I + B B) -1 PB = exp k -tr( L) + tr BPB (I + BB ) -1 BPB = exp k -tr( L) + tr L(I + L) -1 L .

Figure 1 :

 1 Figure 1: First-sample cost for DPP-VFX against other DPP samplers (n is data size).

Figure 2 :

 2 Figure 2: Resampling cost for DPP-VFX compared to the exact sampler of [HKP + 06].

Lemma 7 (

 7 restated Lemma 3) Given L and an arbitrary Nyström approximation L of rank m, computing l i , s, z, and L requires O(nm 2 + m 3 ) time. Proof Given the Nyström set C, let us define the matrix B L I,C L +/2 C

  particular, |S α | can be defined as the sum of Bernoullis |S α | = n i=1 b α,i each distributed according to b α,i ∼ Bernoulli (λ i (αL)/(1 + λ i (αL))) [HKP + 06]. The sum of independent but not identically distributed Bernoullis is a so-called Poisson binomial random variable [H + 56]. More importantly, the following result holds for Poisson binomial random variable. Proposition 5 (D + 64, Theorem 4) Given a Poisson binomial r.v. |S α | with mean k α , let k k α . The mode M α is

  Therefore m = O(k α k 4 ) ≤ O(k 5 ) suffices for the accuracy. Moreover, since s α is parametrized only in terms of the eigenvalues of L, which can be found in O(nm 2 + m 3 ) time, we can compute an α such thats α = 2k 2 +6k+1 2k+6 in O(nk 10 + k 15 ) time, which guarantees k ≤ k α < k + 1 k+2 •

To avoid complicating the exposition with edge cases, we assume k ≥ 1. Note that this can be always satisfied without distorting the distribution by rescaling L by a constant, and is without loss of generality, as our analysis can be trivially extended to the case 0 ≤ k < 1 with some additional notation.

The code used for these experiments is available at https://github.com/LCSL/dpp-vfx.

For simplicity we do not perform the full k-DPP rejection step, but only adjust the expected size of the set.

Acknowledgements

MD thanks the NSF for funding via the NSF TRIPODS program. This material is based upon work supported by the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216, and the Italian Institute of Technology. We gratefully acknowledge the support of NVIDIA Corporation for the donation of the Titan Xp GPUs and the Tesla k40 GPU used for this research.

A Omitted proofs for the main algorithm's guarantees

In this section we present the proofs omitted from Sections 2 and 3, which regarded the correctness and efficiency of DPP-VFX. We start by showing that multiple samples drawn using the same Nyström approximation are independent.

Lemma 5 (restated Lemma 1) Let C ⊆ [n] be a random set variable with any distribution. Suppose that S 1 and S 2 are returned by two executions of DPP-VFX, both using inputs constructed from the same L and L = L I,C L + C L C,I . Then S 1 and S 2 are (unconditionally) independent. Proof Let A and B be two subsets of [n] representing elementary events for S 1 and C, respectively. Theorem 2 implies that

Now, for any

representing elementary events for S 1 and S 2 we have that

Pr(C = B).

Since B∈[n] Pr(C = B) = 1, we get that S 1 and S 2 are independent.

We now bound the precompute cost, starting with the construction of the Nyström approximation L.

Lemma 6 (restated Lemma 2) Let L be constructed by sampling m = O(k 3 log n δ ) columns proportionally to their RLS. Then, with probability 1 -δ, L satisfies the assumption of Theorem 3.

Proof Let L = BB and L = BPB, where P is a projection matrix. Using algebraic manipulations we get s = tr(L -L + L(I + L) -1 ) = tr(B(PB BP + I) -1 B).

The PB BP matrix in the above expression been recently analyzed by [CCL + 19] in the context of RLS sampling who gave the following result that we use in the proof.

Proposition 3 (CCL + 19, Lemma 6) Let the projection matrix P be constructed by sampling O(k log( n δ )/ε 2 ) columns proportionally to their RLS. Then,

We decompose the condition into two parts. In the first part, we bound k -s ≤ 1/2, and then bound z -tr BP(I + B B) -1 PB ≤ 1/2, see the proof of Theorem 3 for more details. It is easy to see that applying the construction from Proposition 3 leads to the following bound on s, s = tr(B(PB

Tuning ε = 1/(2k + 2), we obtain s ≤ k + 1/2 and reordering gives us the desired accuracy result. Similarly, we can invert the bound of Proposition 3 to obtain

and therefore, tr

where the last inequality is due to the fact that L L since it is a Nyström approximation and that the operator tr L(L + I) -1 is monotone. With the same ε as before, we obtain the bound. Summing the two 1/2 bounds gives us the result.

Finally, we show how to compute the remaining quantities needed for DPP-VFX (Algorithm 1).


G Haiat

M C Phan Huy

E Barthel

The adhesive contact of viscoelastic spheres

Keywords: Adhesion and adhesive, Crack propagation and arrest, Contact mechanics, Viscoelastic material, Integral transforms

We have formulated the restricted self-consistent model for the adhesive contact of linear viscoelastic spheres. This model is a generalization of both the Ting (J. Appl. Mech. 33 (1966) 845) approach to the viscoelastic contact of adhesionless spheres and the restricted self-consistent model for adhesive axisymmetric bodies. We also show how the model can be used in practice by giving a few examples of numerical solutions.

Introduction

The contact problem as we consider it nowadays seems to have ÿrst taken shape in the last quarter of the 19th century when [START_REF] Hertz | Ueber die beruehrung fester elastische koerper[END_REF] and [START_REF] Boussinesq | Application des Potentiels  a l'Etude de l'Equilibre et du Mouvement des Solides Elastiques[END_REF] presented solutions to the adhesionless contact of linear elastic spheres.

In the 1930s, the ÿeld of surface forces was mature enough for the problem of the adhesion of rigid bodies and its relation to the surface interactions to be tackled by [START_REF] Bradley | The cohesive energy between solid surfaces and the surface energy of solids[END_REF] and [START_REF] Derjaguin | Untersuchungen ueber die reibung und adhaesion[END_REF].

In the 1960s, the problem of the adhesionless contact of viscoelastic spheres arose. The additional di culty over the elastic case is due to the dependence of the solution upon the history of the system. In particular, the outward run (decreasing contact radius) proved di cult to solve, although a number of approaches were proposed. The deÿnitive solution was proposed by [START_REF] Ting | The contact stresses between a rigid indenter and a viscoelastic half-space[END_REF], and interest in the problem waned soon afterwards.

The 1970s saw a renewed interest in the adhesive contact of elastic bodies. The famous controversy between the Russian (DMT) [START_REF] Derjaguin | On the role of molecular forces in contact deformation[END_REF] and English (JKR) [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF]) schools ended when it was ÿnally recognized that both the models were indeed limiting cases within a more general model-the restricted self-consistent model, a terminology coined by [START_REF] Hughes | Some applications of continuum mechanics in colloid and interface science[END_REF], to our knowledgein which the interactions between the surfaces are explicitly accounted for in a more physical manner [START_REF] Tabor | Surface forces and surface interactions[END_REF]. The transition between these limit cases as the range of the interactions goes from inÿnity to zero was well demonstrated by numerical calculations [START_REF] Muller | On the in uence of molecular forces on the deformation of an elastic sphere and its sticking to a rigid plane[END_REF][START_REF] Muller | General theoretical considerations on the in uence of surface forces on contact deformations[END_REF][START_REF] Greenwood | Adhesion of elastic spheres[END_REF].

In the 1990s, the ÿrst analytical model of the transition was proposed by [START_REF] Maugis | Adhesion of spheres: the JKR-DMT transition using a dugdale model[END_REF], who included the interaction between surfaces through a Dugdale model; in this model, the DMT to JKR transition is spanned when the normalized interaction parameter goes from zero to inÿnity. Soon afterwards, the usefulness of the Maugis model was evidenced in a few experimental cases [START_REF] Maugis | JKR-DMT transition in the presence of a liquid meniscus[END_REF][START_REF] Barthel | Adhesion energy measurements in the presence of adsorbed liquid using a rigid surface force apparatus[END_REF][START_REF] Lantz | Atomic-force microscope study of contact area and friction on NbSe2[END_REF]. Simultaneously, the limitations of the Maugis model were pointed out, and the model was generalized in di erent ways (Barthel, 1998a[START_REF] Barthel | The adhesive contact of spheres: when the interaction is complex[END_REF][START_REF] Greenwood | An alternative to the Maugis model of adhesion between viscoelastic spheres[END_REF].

In the ÿnal years of the last century, the trend for the measurement of local mechanical properties of solids triggered interest in the adhesive contact of viscoelastic solids. Indeed, by reducing the size of the indenter, one unavoidably enters a regime in which the surface forces are signiÿcant. Thus, the coupling between mechanical response and surface interactions has to be taken into account.

It is interesting to note that the limit cases of the elastic adhesive contact problem turn out to be of limited interest in the case where the bulk response is viscoelastic. In the DMT limit (for small ), indeed, the adhesive interaction stress is negligible compared to the mechanical compliance of the solid. As a result, in this model, the coupling between mechanical compliance and surface interactions is limited to the addition of a global adhesive force, so that the richness of the basic phenomenon cannot be accounted for. In contrast, in the JKR limit (for large ), the gradient of the surface displacement is singular at the periphery of the contact zone, so that the strain rate is locally inÿnite, a shortcoming which deprives this limit of usefulness in the viscoelastic case.

Indeed, it is actually required that the viscoelastic response be taken into account to model the process zone. This was treated in depth by Schapery (1975a-c), with a view to dealing with crack extension in viscoelastic solids, and subsequently for crack healing [START_REF] Schapery | On the mechanics of crack closing and bonding in linear viscoelastic media[END_REF]. His results lie at the root of further work on viscoelastic cracks by [START_REF] Greenwood | The mechanics of adhesion of viscoelastic solids[END_REF], [START_REF] Johnson | Contact mechanics and the adhesion of viscoelastic spheres[END_REF] and [START_REF] Lin | Viscoelastic contract (sic), work of adhesion and the JKR technique[END_REF].

For the adhesive viscoelastic contact, however, the crux of the problem is to treat both the crack zone and the contact zone within a viscoelastic model. As in the adhesionless case, the formulation of the inward run (increasing contact radius) is reasonably straightforward. [START_REF] Hui | Contact mechanics and adhesion of viscoelastic spheres[END_REF] have given explicit results for the growing contact area case. In contrast, the decreasing contact radius (outward run), although of great signiÿcance, because it controls the contact rupture and thus the adherence, is still unsolved. Attard (2001a, b) has recently implemented a completely numerical model, a viscoelastic version of the full self-consistent model in the special case of an exponential relaxation function, based on the canonical relation between the exponential function and its derivative. [START_REF] Lin | The role of viscoelastic adhesive contact in the sintering of polymeric particles[END_REF] have proposed an analytical formulation, but this has up to now turned out to be of limited applicability.

In this paper, we show that the formulation of the adhesive contact of elastic bodies we proposed earlier lends itself to the description of the linear viscoelastic case as well. In addition to forming an extension of the restricted self-consistent model for the adhesion of elastic bodies, the formulas we thus obtain reduce to the results of [START_REF] Ting | The contact stresses between a rigid indenter and a viscoelastic half-space[END_REF] when the adhesive interaction is suppressed, thus providing a consistent description of adhesive contact problems. In particular, the decreasing contact radius case is treated on the same footing as the simpler increasing case.

We will ÿrst recall the basic equations for the restricted self-consistent model, in the linear elastic case. We then extend this approach to the linear viscoelastic case. Finally, we give an example of numerical solutions, for a strongly viscoelastic contact.

The restricted self-consistent model for the adhesion of elastic bodies

The self-consistent formulation of the adhesive contact problem has been reviewed recently [START_REF] Huguet | Surface forces and the adhesive contact of axisymmetric bodies[END_REF]. Here, we will only recall the main points:

• the gap between the surfaces geometrically depends upon the surface displacement;

• the surface displacement depends upon the surface stress through the mechanical behaviour of the solid; • the surface stress depends, among other parameters, upon the gap through the nature of the physical interaction between the surfaces.

These three interdependent relations must be satisÿed simultaneously, which results in a self-consistent solution.

The restricted self-consistent model

We now further assume that the range of the repulsive part of the interaction is much shorter than the range of the attractive part, as usually observed. This assumption leads to the restricted self-consistent adhesive contact problem. The boundary conditions are:

u(r) = -f(r) for r ¡ a;
(1) (r) known for r ¿ a;

(2) where u(r) is the normal surface displacement, is the penetration, a the radius of the contact zone and (r) the surface normal stress distribution. The shape of the indenter f(r) is convex and smooth but otherwise not speciÿed. The details of the physics of the adhesive process dictate the self-consistency equation. The diversity of the adhesive phenomena leads to di erent possible self-consistency equations. We here use a simple but approximate one, which is easily introduced if the adhesive process derives from an interaction potential V between the two surfaces, though it is more general. Then, the adhesion energy is

w ≡ V (+∞) -V (0) = +∞ 0 d z dV d z = - +∞ 0 d z (z): (3)
Introducing the gap

h(r) = u(r) -+ f(r) (4)
and applying the change of variable z = h(r), the self-consistency equation (Eq. ( 3)) is

w = - +∞ a (r) dh(r) dr dr: (5)
The strategy for ÿnding the solution to the self-consistent adhesive problem is to converge to a stress distribution outside the contact zone (Eq. ( 2)), which veriÿes both the boundary conditions inside the contact zone (Eq. ( 1)) and the self-consistency equation (Eq. ( 5)).

For a more detailed discussion of self-consistency equations, refer to Barthel (1998b[START_REF] Barthel | The adhesive contact of spheres: when the interaction is complex[END_REF].

Clearly, the calculation of the gap h(r) from the boundary conditions (Eqs. ( 1) and ( 2)) is central to this self-consistent approach. We will now recall the main results relevant to this point.

Surface elasticity

The necessary relations for the elastic case are essentially derived from the work of [START_REF] Sneddon | Fourier Transform[END_REF][START_REF] Sneddon | The relation between load and penetration in the axisymmetric boussinesq problem for a punch of arbitrary proÿle[END_REF] and [START_REF] Lowengrub | The e ect of internal pressure on a penny-shaped crack at the interface of two bonded dissimilar elastic half-spaces[END_REF]. In contrast to our recent papers, we here introduce a new function, Â, which is necessary in the viscoelastic case, but redundant with the more usual g function in the elastic case. For axisymmetric bodies, the auxiliary function g [resp. Â] is a suitable transform (of the Abel type) of the normal surface stress distribution [resp. the normal surface displacement u]:

g(s) ≡ - +∞ s r (r) (r 2 -s 2 ) dr; (6) Â(s) ≡ d ds s 0 ru(r) √ s 2 -r 2 dr (7) = u(s = 0) + s s 0 dr u (r) (s 2 -r 2 ) : (8)
These transforms are characterized by several properties:

Prop I: they can be analytically inverted (Appendix A.1); Prop II: g(s) depends on (r) for r ¿ s only, and Â(s) depends on u(r) for r ¡ s only;

Prop III: under the conditions of linear elastic behaviour and the absence of shear stresses at the interface, mechanical equilibrium leads to g(s) = KÂ(s) for all s;

where K depends upon the Young modulus E and the Poisson ratio as

K = E 2(1 -2 ) : (10)
Thus, from the boundary condition Eq. ( 1) and Prop. II (Eq. ( 8)), one obtains Â(r) inside the contact zone (r ¡ a):

Â(r) = -0 (r); (11)
where

0 (r) = d dr r 0 sf(s) (r 2 -s 2 ) ds = r r 0 ds f (s) (r 2 -s 2 ) (12)
is a signature of the shape of the indenting bodies only. Similarly, we obtain g(r) outside the contact zone (r ¿ a) from the boundary condition Eqs.

(2) and Prop. II (Eq. ( 6)). The main step in the calculation of the solution to the elastic problem is to obtain the extensions of  and g to, respectively, outside and inside the contact zone. This step is straightforward in the elastic case because of the simplicity of the equilibrium relation Eq. (9) (Prop. III), except for the behaviour of g and  at the contact zone boundary a. We will here require continuity of the g(r) function (and thus of the Â(r) function), which implies continuity of the stress distribution (r) (and therefore of the derivative of the displacement ÿeld u(r)) at r = a. This condition actually provides an additional equation, which we now discuss.

Penetration

From the speciÿc form of  given by Eq. ( 11), and the equilibrium relation Eq. ( 9), there results that the assumption of continuity of the stress distribution at the border of the contact zone a determines the penetration as a function of the shape of the bodies f(r), the contact radius a and the attractive stress distribution outside the contact zone by

= 0 (a) + 1 K g(a): (13)
A problem seems to arise in the JKR case, where it is well known that the stress is discontinuous at r = a. However, it only appears as a limit case of models in which the continuity is obeyed [START_REF] Maugis | Adhesion of spheres: the JKR-DMT transition using a dugdale model[END_REF]Barthel, 1998a). Now that  is completely determined outside the contact zone, we can tackle the calculation of the gap.

The gap and the self-consistency equation

From the inverse of Eq. ( 7), Eq. (A.3), in Appendix A, the surface displacement can be split into

u(r) = 2 a 0 ds Â(s) √ r 2 -s 2 + 2 r a ds Â(s) √ r 2 -s 2 : ( 14
)
It can be shown that a convenient form for the gap (see Appendix A.2) is

h(r) = h Hertz (r; a) + 2 r a Â(s) -Â(a) √ r 2 -s 2 ds (15) with h Hertz (r; a) = 2 1 R f H (r; a); (16)
where

f H (r; a) ≡ r a ds (s 2 -a 2 ) √ r 2 -s 2 = a 2 r 2 -a 2 + r 2 2 -a 2 arccos a r : (17)
This form of the gap will be used below to compute the self-consistency equation in a special case, the general expression being too complex to be of real use.

In addition, from this formulation (Eq. ( 15)), it becomes apparent that, as noted above, the continuity of  entails the suppression of the square root singularity of the gap h at a (cf. Appendix A.2).

Force

Similarly, the total force,

F = 4 +∞ 0 dr g(r) (18)
can be split into two terms. The outer term is

F ext = 4 +∞ a dr g(r): (19)
The inner term is

F int = 4 a 0 dr g(r) = 4K(a -0 (a)) (20) = F H (a) + F JKR (a) (21)
with

F H (a) = 4K(a 0 (a) -0 (a)); (22) F JKR (a) = 4ag(a); ( 23
)
where F H is the non-adhesive Hertz force which, like the function

0 (r) ≡ r 0 ds 0 (s) (24)
depends only upon the shape of the indenter, and F JKR is the corrective JKR-like at punch term.

Resolution of the adhesive contact problem

In the experiments, two main possibilities are met, where the displacement (ÿxed grip) or the force (ÿxed load) are prescribed. Assuming a self-consistent model as in Eq. ( 5), the problem is thus to determine the contact radius a and the attractive stress distribution by (a) in a ÿxed grip experiment, solving the displacement (Eq. ( 13)) and self-consistency (Eq. ( 5)) equations simultaneously, and then calculating the force; (b) in a ÿxed load experiment, solving the force (Eq. ( 18)) and self-consistency (Eq.

( 5)) equations simultaneously and then calculating the displacement.

Examples of such procedures can be found in Barthel (1998b[START_REF] Barthel | The adhesive contact of spheres: when the interaction is complex[END_REF] and [START_REF] Greenwood | An alternative to the Maugis model of adhesion between viscoelastic spheres[END_REF].

Adhesion of viscoelastic bodies

Linear viscoelastic behaviour

Assuming the simplest viscoelastic behaviour-a constant Poisson ratio-the stress and deformation now obey

(t) = t 0 d (t -) d d ( ); ( 25
) (t) = t 0 d (t -) d d ( ); ( 26
)
where the relaxation function (t) and the creep function (t) are inverse for this product of convolution. Therefore we have

(t) = (t 0 ) + t t0 d (t -) 9 9 t0 d ( -) 9 9 ( ) (27) provided (t) is independent of time for t ¡ t 0 .
Since there is no coupling through the response functions and between the temporal and spatial variations of the physical quantities, the equilibrium equation (Eq. ( 9)) for the adhesive contact between a rigid indenter and a viscoelastic solid reads

g(t) = t 0 d (t -) d d Â( ) (28)
or its inverse

Â(t) = t 0 d (t -) d d g( ): (29)
The essential Eqs. ( 5), ( 13), (15), and (18) remain unaltered, except for the introduction of the time dependence. For instance, the self-consistency equation Eq. ( 5) now reads 31)) [resp. Eq. ( 4)] for t ¡ tmax [resp. tmax ¡ t]) and the self-consistency equation (Eq. ( 30)). The shaded rectangles exemplify the r-t zones relevant for the calculation of the gap at r 1 (inward) and r 2 (outward). Like the penetration, calculation of the stress inside the contact zone ( point) is done along a horizontal line. The functions ta -(r) and ta + (r) are deÿned in the text (Eqs. ( 34) and ( 43)).

w = +∞ a(t) (r; t) d dr h(r; t) dr: (30)
The guideline to establishing the adhesive viscoelastic contact equations is thus to express these equations at time t as a function of known quantities at time ¡ t. However, for a given ¡ t, we directly know Â(r; ) for r ¡ a( ) only (from the datum of ( ) and 0 (r)) and g(r; ) for a( ) ¡ r only. This deÿnes two time-space domains, the g-[resp. Â-] domain where g [resp. Â] is known (Fig. 1). The calculation of the extensions of each function (g in the Â-domain and conversely) is now more complex, because the simple proportionality relation Eq. ( 9) is replaced by Eq. ( 28) or its inverse Eq. ( 29), which sometimes must be used recursively, leading to cumbersome nested temporal integrals.

We will therefore restrict ourselves to the two cases which are most signiÿcant: for (a) increasing and (b) subsequently decreasing contact radius. Indeed, as in the adhesionless case, for more complex loading histories, the equations which describe the response of the system become increasingly unwieldy, although they can still be derived with the methods developed in the present paper.

The next two sections will brie y present the general results for the viscoelastic adhesive contact, while the rest of the paper is devoted to the explicit development of a special case as an example.

Increasing contact radius

This is the most straightforward case, since the history of the system is simple enough to easily be taken into account: for a given radius r, we meet at most once with the boundary B between the two domains (Fig. 1). Except for the stress distribution inside the contact zone, and thus the force (Section 3.2.4), all relevant quantities can be expressed with one single temporal integral.

Displacement

From Eq. ( 29), the continuity of  at the edge of the contact zone, which determines , now yields

(t) = 0 (a(t)) + t 0 d (t -) 9 9 g(a(t); ); (31)
which is the form of Eq. ( 13) valid for the viscoelastic case in the increasing contact radius regime.

The gap

Now in the gap equation Eq. ( 15), we have

Â(r; t) = t 0 d (t -) 9 9 g(r; ) (32)
and thus

h(r; t) = h Hertz (r; t) + 2 t 0 d (t -) 9 9 r a(t) ds g(s; ) -g(a(t); ) √ r 2 -s 2 (33)
from which the self-consistency equation is derived through Eq. (30).

The stress distribution inside the contact zone

Inside the contact zone (r ¡ a(t)), the stress distribution can be calculated from the successive use of Eqs. ( 28) and (11). However, in Eq. ( 28), the time domain extends from 0 to t, while Â(r; t) is known only inside the Â-domain, that is between t a-(r) and t, where t a-(r) is deÿned in the increasing contact radius zone (Fig. 1) by a(t a-(r)) = r:

(34)

Therefore, outside the Â-domain, we have to express Â(r; ) as a function of g(r; ) through Eq. ( 29), from which the following nested integral expression results:

g(r; t) = t ta -(r) d (t -) 9 9 { ( ) -0 (r)} + ta -(r) 0 d (t -) 9 9 0 d ( -) 9 9 g(r; ): (35)
Note that this expression is actually also valid in the decreasing contact radius regime, where it lies at the root of the equation giving the penetration.

The force

The total force

F(t) = F ext (t) + F int (t) (36)
is the sum of the outer term

F ext (t) = 4 ∞ a(t)
dr g(r) (37)

and the inner terms

F int (t) = F int;a (t) + F int;b (t) (38)
with

F int;a (t) = 4 t 0 d (t -) d d min(a(t);a( )) 0 dr ( ( ) -0 (r)) (39)
and, in a manner similar to the previous paragraph,

F int;b (t) = 4 ta -(a(t)) 0 d (t -) d d 0 d ( -) d d a(t) a( ) dr g(r; ): (40)
Here again, the expression for the force is valid for both in-and outward runs.

Decreasing contact radius

In the decreasing contact radius regime, we may meet up to twice with the boundary B (Fig. 1), from which increasingly complex expressions result.

Displacement

Using Eq. ( 35), continuity of g at the periphery of the contact now gives

g(a(t); t) = t ta -(a(t)) d (t -) 9 9 { ( ) -0 (a(t))} + ta -(a(t)) 0 d (t -) 9 9 0 d ( -) 9 9 g(a(t); ) ; (41)
which is the form of Eq. ( 13) valid in the decreasing contact radius regime. We note that this equation is simply an extension of the integral equation characteristic of the decreasing contact radius case in the adhesionless contact [START_REF] Ting | The contact stresses between a rigid indenter and a viscoelastic half-space[END_REF], which is recovered here for g = 0. Thus, as in the non-adhesive case, the penetration is to be calculated from this integral equation.

In addition, we note that at t = t a-(a(t)), that is to say when the contact radius reaches its maximum a max = a(t max ), the second term on the right-hand side of Eq. ( 41) is shown through Eq. ( 27) to equal g(a(t); t). Thus the convolution integral of ( ) is zero, with the simple but prominent implication that (t) is continuous at t max .

The gap and the self-consistency relation

In the calculation of the gap, for a given radius we meet twice with the boundary B between domains, leading to an increasing number of nested temporal integrals. As an example of the recursive use of Eqs. ( 28) and ( 29) already applied in Section 3.2.3, we here explicitly give the suitable expression for  to calculate the gap, in which each term is expressed as a function of known quantities (the displacement in the Â-domain, the stress distribution in the g-domain):

Â(s; t) = ta -(s) 0 d (t -) 9 9 g(s; ) + t ta + (s) d (t -) 9 9 g(s; ) + ta + (s) ta -(s) d (t -) 9 9 ta -(s) d ( -) 9 9 Â(s; ) + ta -(s) 0 d ( -) 9 9 0 d ( -) 9 9 g(s; ) ; ( 42
)
where t a+ (r) is deÿned in the decreasing contact radius zone (Fig. 1) by

a(t a+ (r)) = r: (43)
Inverting the order of the temporal and the spatial integrals, we obtain an expression for u where all terms can be calculated as a function of the history of the system.

Here again, from the gap expression, the self-consistency equation is derived through Eq. ( 30). These expressions are given explicitly in Appendix B (Eq. (B.42)).

The resolution of the adhesive viscoelastic contact problem

In this case, the solution proceeds as follows: one has to simultaneously solve the penetration equation (Eqs. ( 31) or (41) -the latter one being an integral equation) [resp. the force equation (Eq. ( 36))] and the self-consistency equation Eq. ( 30) for a ÿxed grip [resp. ÿxed load] experiment. We now show that this approach is feasible in practice with the following examples.

Double-Hertz model

To model experiments with a sphere, we assumed a paraboloidal shape of the indenter. Other shapes can naturally be dealt with, and nano-indentation addicts may here assume a sphere-terminated cone, and modify 0 and 0 accordingly. As for the interaction zone, we at ÿrst implemented a Dugdale model. The resulting model was thus the viscoelastic Dugdale-Maugis model. In practice however, we found it much more advantageous to avail ourselves of the so called 'double-Hertz' model, put forward by [START_REF] Greenwood | An alternative to the Maugis model of adhesion between viscoelastic spheres[END_REF], where the attractive stress distribution inside the interaction zone, bounded by the contact radius a and the interaction radius c, is assumed to be ellipsoidal. Indeed, in this model, we beneÿt both from the absence of stress discontinuity at the periphery of the interaction zone and from the analyticity of all the spatial integrals. As a result, the numerical calculations are signiÿcantly easier. All the equations used are given explicitly in Appendix B.

A program has been implemented to work with any pair of viscoelastic functions, provided they are inverse for the convolution integrals Eqs. ( 28) and ( 29). In the 44)) and relaxation (Eq. ( 45)) functions used for the results shown below. The characteristic creep time is the time unit, and the characteristic creep parameter is k = 0:09. following, we assume a normalized standard viscoelastic solid model with creep and relaxation functions given by

(t) = 1 + k 1 -k (1 -exp(-t)); ( 44
) (t) = k + (1 -k) exp(-t=k); ( 45
)
where k is the characteristic parameter of the normalized creep process. The relaxation and creep functions for k = 0:09, as assumed below, are illustrated in Fig. 2. Note that for this value of k, the relaxed modulus is about ten times smaller than the instantaneous modulus.

As mentioned in the Introduction, in normalized form, in the elastic case (paraboloidal indenter with radius R, adhesion energy w: see Eqs. (B.19)-(B.22) for the details of the normalization), the adhesive interaction is characterized by a single parameter . We here use the same normalization, where the instantaneous modulus (0) is now substituted for the elastic modulus.

As an example, let us assume (0)=100 MPa, R=2 cm and w=50 mJ m -2 ; we then have the following normalizing factors: force, P 0 3 mN; contact radius, A 0 90 m; penetration: 0 0:4 m and stress: 0 0:4 MPa.

In the following section, we discuss numerical results obtained for an adhesion parameter = 4. Thus, we present examples for a rather extreme (and computationally more demanding) case of notably viscous solid with sizeable adhesion.

We assume a ÿxed grip experiment, and simultaneously solve Eqs. (B.30) and (B.37) (inward run), or (B.38) and (B.56) (outward run), for the normalized radius of the contact zone A(t) and the normalized radius of the interaction zone C(t). We suppose zero penetration and an instantaneous switch-on of the interactions at t=0. The temporal integrals were calculated with a Simpson rule method, with a convergence criterion of 10 -4 . The three-integrals term has been neglected throughout, due to long computation times. We checked in a few cases that inclusion of this term leads to minor changes in the contact radius history at the turning point (maximum contact radius), and is not signiÿcant at later stages.

Numerical results

Example 1: adhesionless contact-adhesive contact

We here contrast the adhesive and adhesionless systems. The indenter is driven with a simple penetration history (Fig. 3) with a penetration velocity of +10 (inward part of the curve-zone 1) or -10 (outward-zone 3), normalized to the characteristic creep time. An intermediate period of reduced penetration velocity (+2-zone 2) is provided for, so as to better illustrate the response of the system.

Between these two cases, the main di erence in the inward part of the run (zones 1 and 2) lies in the contact radius (Fig. 4), which, for an identical penetration, is increased through the interactions. However, in the inward run, the force curves for the two cases are quite similar, featuring in particular a sharp relaxation process in zone 2 (Fig. 5). The one noteworthy di erence is to be found at the very beginning of the curve, where the initially adhesive (i.e. negative) total force, in the adhesive case, swiftly decays into repulsive.

In the third zone, the overall di erences between the two cases are more striking. When the sphere is retracted, a sharp decrease of the contact radius (Fig. 4) to zero occurs for the adhesionless contact, while the contact radius goes on increasing and then stays roughly constant for one unit of time when adhesion is present, before it ÿnally bends down. In the constant contact radius region, the force vs. penetration curve (Fig. 6) is found to be linear. Finally, the maximum tensile force is -13:4, considerably larger than the bare -1:5 expected in an elastic model for such a large value of . Thus in the viscoelastic case, the contact rupture is characterized by two parameters: the "stick time" between the inception of retraction and the decrease of the contact radius, and the (large) adhesive force. In zones 1 and 2, the contact radius increases under the e ect of the interactions, in a manner similar to the classical elastic contact. In zone 3 (outward), in the adhesive contact, the contact radius stays roughly constant for about one time unit before curving down, in sharp contrast to the adhesionless case. We will call this constant contact radius domain the "stick" domain.

- In the adhesive and adhesionless cases, the forces are similar in zones 1 and 2 (inward), due to fast stress relaxation, most evident in zone 2. In the adhesive case however, the contact zone starts receding (in zone 3) only if su cient tensile stresses have built up in the contact zone, which results in a large adhesive force.

Example 2: in uence of the retraction velocity

In the second series of calculations, we changed the retraction velocity, for identical inward runs (Fig. 7). The stick time increases markedly when the retraction velocity decreases (Fig. 8). The adhesive force (Fig. 9), however, is little a ected (an apparent dependence as the power 0.1 of the contact radius velocity at maximum tensile force). This behaviour, which cannot be explained by crack tip e ects only, will be qualitatively explained in the discussion.

Example 3: shallower indents

Everyday experience or laboratory experiments evidence an increase of the adhesive force with dwell time or maximum penetration depth. The present system features such a behaviour as well. Using the penetration velocities displayed in Fig. 10, decreasing the maximum penetration, which occurs at t = 2, from 20 to 6, the maximum contact radius (Fig. 11) drops from about 5 to about 3, and with identical retraction velocities, we calculate a decrease of the adherence force (Fig. 12) from 13.1 to 8.6. A weak decrease of the adherence force with decreasing retraction velocity (-15:8 for v = 30, -14:9 for v = 15, -13:7 for v = 7:5) is calculated, while the penetration decreases markedly, to negative values for the two slower velocities, because of the competition between creep under mainly compressive (zones 1 and 2) or mainly tensile (zone 3) stresses.

Discussion

Crack tip e ects-creep

The typical size of the interaction zone is of the order 5 × 10 -2 , to be compared with the 8 × 10 -2 expected for an elastic solid with the instantaneous modulus. Typical contact radius velocities are in the range 1-4, except in the stick zone, where it goes to zero. Thus, in the crucial zones-zone 1, which determines the maximum contact radius, and the rupture zone (later stage of zone 3)-, the characteristic time of the interaction stress is of the order of a few 10 -2 , well below the creep time (unity). As a result, the interaction zone behaves essentially elastically, with an e ective modulus close to the instantaneous value. It is only in the stick zone (constant contact radius) that, due to reduced crack tip velocities, viscoelastic crack tip e ects become apparent mainly under the form of a slow increase of the contact radius, due to creep, while the interaction zone contracts by a factor of the order of unity (Fig. 13), because the solid is e ectively softer at long time scales.

Contact zone e ects-stress relaxation

On the other hand, most of the e ects described above may be ascribed to stress relaxation within the contact zone. This relaxation is particularly noticeable in zone 2, where the force decay has a ÿtted characteristic time of 0.089 (Fig. 5), to be compared -13:1) is due to the smaller maximum contact radius: for small penetrations, the contact radius relevant for the calculation of the adherence force is mainly controlled by the maximum contact radius. with the theoretical k =0:09 (Eq. ( 45)). Similarly, in the region where the contact radius is about constant, the stresses relax rapidly and the system actually behaves as a at punch of radius equal to the contact radius a with a modulus equal to the inÿnite-time modulus K ∞ . Indeed, in this case, the predicted slope (or contact sti ness) is

Ä = 3 2 K ∞ a; (46)
which in the present case is equal to 0.69, identical to the slope in the calculated solution. Similarly, the stick zone, where the contact radius is constant, results from the weak elastic energy build-up in the contact zone. It is only for a sizeable backward displacement that tensile stresses at the crack tip ÿnally meet up with the adhesive interactions. The stick time can thus be estimated from the following argument: in Eq. ( 41), the integral

I = t ta -(a(t)) d (t -) 9 9 { ( ) -0 (a(t))} (47)
is an expression of (the Abel transform of) the stress at a(t) induced by the retraction motion ( ), where the last (double-integral) term has been neglected. It has to equal the (Abel transform of the) stress generated by the interactions, which is g(a(t)). Now, taking rather arbitrarily an interaction zone size of 4 × 10 -2 -not quite the elastic limit, since we are still in a low crack velocity region (Fig. 13), a contact radius of 8, Eq. (B.23) yields a value g(a(t)) 3. Now, during the "stick" phase, the penetration is monotonously decreasing. Then, taking into account the present form for the relaxation function Eq. ( 45), for k1,

I k d dt t -t max + 1 -exp - t -t max k : (48)
From this equation and g(a(t)) = 3, we estimate that t -t max is, respectively, 0.2, 1 and 3 for penetration velocities of 30, 15 and 7.5, in reasonable agreement with the stick times observed in Fig. 8. Along the same line of thought, in the outward run, reduction of the crack size requires the build-up of su cient tensile stresses inside the contact zone. However, in contrast to an elastic behaviour, where compressive and tensile stresses evolve simultaneously, in the viscoelastic case, the compressive contact zone stresses have relaxed and tensile stresses dominate in the outward leg. Thus, the balance between compressive and tensile stresses found in the elastic case is o set for a viscoelastic contact, which entails an adhesive force much larger than expected in a completely elastic problem, along with a reduced release of elastic energy. Thus, stress relaxation inside the contact zone is the main source of energy dissipation.

More precisely, in a JKR-like regime, where the interaction zone is small, the force is dominated by the inner term in

P 3 2 t 0 d ˜ (t -) d d A(t) ( ) - A(t) 3 3 ; (49) while G(A) t 0 d ˜ (t -) d d ( ( ) -A(t) 2 ); (50)
which combine to give an expression quite similar to the one proposed by [START_REF] Johnson | Contact mechanics and the adhesion of viscoelastic spheres[END_REF]:

P ˜ (t)A(t) 3 + 3 2 G(A(t))A(t) (51)
with

G(A) = -3 C 2 -A 2 : (52)
decrease. This "sticking" phenomenon is due to the fast relaxation of the stresses (both compressive and tensile) inside the contact zone, so that crack back-propagation requires the renewal of tensile stresses inside the contact zone, as described by Eq. ( 48). As a result, the adherence force is magniÿed because the tensile stresses needed for crack retraction are no longer balanced by equivalent compressive stresses, since the latter have already relaxed (Eq. ( 55)). Thus, we ÿnd that the ampliÿcation of the adherence force occurs because strain and stress inside the contact zone are out of phase, as typical of viscoelasticity. Crack tip e ects play a more subdued role.

Ongoing work deals with the development of approximate models and the incorporation of roughness therein.

Appendix A. General relations

A.1. Inversion of the auxiliary functions g and Â

The inverse relations to Eqs. ( 6) and ( 7

A.2. General expression for the gap

We derive Eq. ( 15) from Eq. ( 14)

h(r) ≡ u(r) -+ f(r) (A.4) = 2 a 0 ds -0 (s) √ r 2 -s 2 -2 ( -f(r)) + r a ds Â(s) √ r 2 -s 2 (A.5) = 2 r a ds 0 (s) - √ r 2 -s 2 + r a ds Â(s) √ r 2 -s 2 (A.6) = 2 r a ds 0 (s) -0 (a) √ r 2 -s 2 + r a ds Â(s) -Â(a) √ r 2 -s 2 ; (A.7)
where the equality of Eq. (A.5) and (A.6) results from the inversion of the deÿnition of 0 (cf. Eq. ( 6) and its inverse Eq. (A.2)). In the case of the paraboloidal indenter, the ÿrst term is the Hertz gap Eq. ( 16).

In a double-Hertz model,

F ext = 0 √ c 2 -a 2 - 2 3 c 3 - 1 3 a 3 + c 2 a : (B.6) Altogether, F(a) = 8Ka 3 3R - 2 0 3 c 3 -a 3 √ c 2 -a 2 : (B.7) B.1.4.
The gap With Eqs. ( 15) and ( 9), and using the relation

Y (c -s)(s 2 -c 2 ) -(a 2 -c 2 ) = Y (s -c)(c 2 -s 2 ) + (s 2 -a 2 ); (B.8)
where Y is the Heaviside step function, we have

h(r) = h Hertz (r; a) + 0 2K √ c 2 -a 2 r a (s 2 -a 2 ) -Y (s -c)(s 2 -c 2 ) √ r 2 -s 2 ds (B.9) = h Hertz (r; a) + 0 2K √ c 2 -a 2 r a (s 2 -a 2 ) √ r 2 -s 2 -Y (r -c) r c (s 2 -c 2 ) √ r 2 -s 2 ds (B.10) = h Hertz (r; a) + 0 2K √ c 2 -a 2 {f H (r; a) -Y (r -c)f H (r; c)}; (B.11)
where f H is deÿned by Eq. ( 17). Combining these expressions

h(r) = 2 R + 0 2K √ c 2 -a 2 f H (r; a) - 0 2K √ c 2 -a 2 Y (r -c)f H (r; c): (B.12) B.1.5. Self-consistent approach
The ÿrst term only in Eq. (B.12) contributes to the integral. We have to calculate

= 2 R + 0 2K √ c 2 -a 2 0 √ c 2 -a 2 I (B.16) and w = 1 3 1 R + 0 4K √ c 2 -a 2 0 √ c 2 -a 2 I H (c; a): (B.17) B.1.6. Normalization
In the case of the sphere, we deÿne

K = 8K 3 : (B.18)
Following [START_REF] Maugis | Adhesion of spheres: the JKR-DMT transition using a dugdale model[END_REF], we normalize F by wR, and introduce

P = F wR ; (B.19) A = a ( wR 2 =K) 1=3 ; (B.20) = ( 2 w 2 R=K 2 ) 1=3 ; (B.21) = 2 0 ( wK 2 =R) 1=3 : (B.22) Then G(S) = 3 S 2 -C 2 √ C 2 -A 2 ; (B.23) = A 2 -3 C 2 -A 2 ; (B.24) P = A 3 -3 C 3 -A 3 √ C 2 -A 2 (B.25)
and

H (R) = 2 1 + 3 1 √ C 2 -A 2 f H (R; A) -3 1 √ C 2 -A 2 Y (R -C)f H (R; C) ; (B.26) 1 = 6 1 + 3 1 √ C 2 -A 2 I H (C; A) √ C 2 -A 2 : (B.
(r; t) =        -0 c(t) 2 -r 2 c(t) 2 -a(t) 2 if a(t) ¡ r ¡ c(t); 0 i f c(t) ¡ r: (B.28)
From Eq. ( 31), the penetration is given by

(t) = 0 (a(t)) -4 0 t tc -(a(t)) d (t -) 9 9 c 2 ( ) -a 2 (t) c 2 ( ) -a 2 ( ) (B.29) or, in normalized form, (t) = A 2 (t) -3 t tc -(A(t)) d (t -) 9 9 C 2 ( ) -A 2 (t) C 2 ( ) -A 2 ( ) : (B.30) B.2.2. Gap
In this case, Eq. ( 33) becomes h(r; t) = h Hertz (r; t)

+ 0 2 t tc -(a(t)) d (t -) 9 9 [f H (r; a(t)) -Y (r -c( ))f H (r; c( ))] c 2 ( ) -a 2 ( ) (B.31)
or, in normalized form, where M a stands for min(A(t); A( )) and M for min(A(t); C( )).

H (R; t) = 2 f H (R; A(t)) + 3 t tc -(A(t)) d (t -) 9 9 × [f H (R; A(t)) -Y (R -C( ))f H (R; C( ))] C 2 ( ) -A 2 ( ) : (B.

Fig. 1 .

 1 Fig.1. Schematic plot of contact radius vs. time (line B). The r-t plane is divided into two parts: the g-domain (a(t) ¡ r), in which the (attractive) surface stress distribution is known, and the Â-domain (r ¡ a(t)), where the surface displacement is known. A solution is a time dependent-stress distribution which veriÿes the penetration equation (Eq. (31)) [resp. Eq. (4)] for t ¡ tmax [resp. tmax ¡ t]) and the self-consistency equation (Eq. (30)). The shaded rectangles exemplify the r-t zones relevant for the calculation of the gap at r 1 (inward) and r 2 (outward). Like the penetration, calculation of the stress inside the contact zone ( point) is done along a horizontal line. The functions ta -(r) and ta + (r) are deÿned in the text (Eqs. (34) and (43)).

Fig. 2 .

 2 Fig. 2. Exponential creep (Eq. (44)) and relaxation (Eq. (45)) functions used for the results shown below. The characteristic creep time is the time unit, and the characteristic creep parameter is k = 0:09.

Fig. 3 .

 3 Fig. 3. Example 1-Penetration history imposed on the system. Inward (fast (+10): zone 1-slow (+2): zone 2) and outward (fast (-10): zone 3) normalized penetration velocities.

Fig. 4 .

 4 Fig.4. Example 1-Contact radius as a function of time. In zones 1 and 2, the contact radius increases under the e ect of the interactions, in a manner similar to the classical elastic contact. In zone 3 (outward), in the adhesive contact, the contact radius stays roughly constant for about one time unit before curving down, in sharp contrast to the adhesionless case. We will call this constant contact radius domain the "stick" domain.

Fig. 5 .

 5 Fig.5. Example 1-Force as a function of time. In the adhesive and adhesionless cases, the forces are similar in zones 1 and 2 (inward), due to fast stress relaxation, most evident in zone 2. In the adhesive case however, the contact zone starts receding (in zone 3) only if su cient tensile stresses have built up in the contact zone, which results in a large adhesive force.

Fig. 6 .Fig. 7 .

 67 Fig.6. Example 1-Force as a function of penetration: note the linear regime in the "stick" domain. The theoretical slope for a at punch with the relaxed modulus (Eq. (46)) is indicated by the line segment under the curve.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Example 2-Contact radius as a function of time. The points at which the maximum tensile (i.e. adherence) force is reached are indicated. Note the variation of the "stick time" with retraction velocity (see text for a qualitative analysis).

Fig. 10 .Fig. 11 .

 1011 Fig.10. Example 3-Penetration history imposed on the system: identical retraction velocity, but two di erent maximum penetrations.

Fig. 12 .

 12 Fig.12. Example 3-Force as a function of penetration. The weaker adherence force (-8:6 vs. -13:1) is due to the smaller maximum contact radius: for small penetrations, the contact radius relevant for the calculation of the adherence force is mainly controlled by the maximum contact radius.

Fig. 13 .

 13 Fig.13. Evidence for crack tip e ects (taken from Example 1). The contact zone radius (dashed-left-hand scale, note the high magniÿcation) continues increasing well into the decreasing displacement zone (zone 3) as a result of creep at the crack tip. The interaction zone extension c -a (line-right-hand scale) decreases-weakly-with decreasing crack tip velocity da=dt.

F

  gap and the self-consistency relation Due to the intertwining of time and space dependances, we introduce an equivalent form for the gap h(r) = h Hertz (r; a) + u JKR (r; a) + u ext (r; a); (B.39) where in Eq. (15), we single out the JKR term in the gap equation u JKR (r; a) = -s 2 ds = f d (r; d; e) -Y (r -f)f d (r; d; f) the case of the double-Hertz model, the three spatial integrals in Eqs. (B.43)-(B.45) are, respectively,I ext1 = -4 0 f t (r; c 2 ( ); max(a(t); a( )); c( )) c 2 ( ) -a 2 ( ) ; (B.48) I int1 = f t (r; ( ); a(t); min(a( ); a( ) 2 -r 2 9 9r h(r; t) dr: (B.51) The core of the calculation is I d (a; c; d; e) = I H (C(t); A(t)) -2 (C(t) -A(t))( (t) -0 (A(t))) + A int1 (B.57) andT = A ext1 + A int2 : ) -A 2 ( ) {I d (A(t); C(t); C 2 ( ); max(A(t); A( ))) -Y (C(t) -C( ))I d (A(t); C(t); C 2 ( ); C( ))} ; ×{I d (A(t); C(t); ( ); A(t))) -Y (c(t) -min(A( ); A( ))I d (A(t); C(t); ( ); min(A( ); A( )))} (B.60) ) -A 2 ( ) ×{I d (A(t); C(t); C 2 ( ); max(A(t); A( ))) -Y (C(t) -A( ))I d (A(t); C(t); C 2 ( ); A( ))}: (B.61) B.4. Viscoelastic bodies-the force As in Section B.1.3, the total force is the sum of the outer term F ext (t) = 0 c 2 (t) -a 2 (t) int (t) = F int;a (t) + F int;b t) + A(t)C 2 (t)

If the interaction zone may be modelled as elastic, then 1 2 12 2 (C(t) -A(t)):

Assuming fast stress relaxation, (t) (+∞) (54) and

which is the normalized JKR equation with relaxed Hertzian stresses. The maximum force is obtained for a radius

and is

In our example, this approximation results in an adherence force P max equal to -16:7, and a contact radius at the maximum force A Pm of 5.7, consistent with the values computed for large penetrations (Example 2). What happens, then, for penetrations smaller than A Pm (Example 3)? The contact radius at maximum tensile force is of course no longer given by Eq. (56). It is not exactly given by the maximum contact radius either because the interaction zone extension, and |g(a)|, increase when the contact radius velocity increases. However, at the same time, the contact radius decreases. It is a balance of these two competing e ects which determines the contact radius at maximum tensile force. The approximation (Eq. ( 55)) may still be useful, however. Indeed, substituting the calculated contact radius values in Eq. ( 55), we obtain approximate adherence forces of, respectively, -9:3 and -14:7 for computed values of -8:6 and -13:1.

Conclusion

In the ÿrst part of the paper, we have shown that the restricted self-consistent model for the adhesive contact of linear elastic bodies lends itself to generalization to linear viscoelasticity. In the relatively simple model developed in the present paper, two scalar equations (expressing the continuity of the stress distribution at the contact zone boundary and the self-consistency of the description of the adhesive stresses) have to be solved simultaneously for the contact zone radius and the interaction zone extension as a function of time. Complex histories of the contact radius lead to increasing numbers of nested time integrals in these equations, with decreasing practical utility for numerical calculations.

In the second part of the paper, a case study, we chose a highly viscoelastic response. We then evidence a "pinning" of the contact zone in the initial stages of the penetration Let us investigate the behaviour of the gap h(r) for r → a(t). For r = a(1 + ), one can show that, for a di erentiable function j(s), r a

Thus, continuity of  at a entails that the second term in Eq. ( 15) behaves as 3=2 at a. The Hertz term h Hertz (r; a) behaves similarly, since it as the same form (Eq. ( 16)), and therefore the O( 3=2 ) behaviour applies to the total gap.

Appendix B. Adhesive contact models

We derive the main results of the double-Hertz model in the elastic case before extending the results to the viscoelastic case.

B.1. Elastic case B.1.1. Stress

In the double-Hertz model, the normalized stress distribution outside the contact zone is ellipsoidal:

We calculate g outside the contact zone: from Eq. (A.2) we infer that for a 6 r 6 c