On modeling and measuring viscoelasticity with dynamic Atomic Force Microscopy

The interaction between a rapidly oscillating atomic force microscope tip and a soft material surface is described using both elastic and viscous forces with a moving surface model. We derive the simplest form of this model, motivating it as a way to capture the impact dynamics of the tip and sample with an interaction consisting of two components: interfacial or surface force, and bulk or volumetric force. Analytic solutions to the piece-wise linear model identify characteristic time constants, providing a physical explanation of the hysteresis observed in the measured dynamic force quadrature curves. Numerical simulation is used to fit the model to experimental data and excellent agreement is found with a variety of different samples. The model parameters form a dimensionless impact-rheology factor, giving a quantitative physical number to characterize a viscoelastic surface that does not depend on the tip shape or cantilever frequency.


I. INTRODUCTION
An increasingly important application of the Atomic Force Microscope (AFM) is the characterization of viscoelastic materials and interfaces, such as cell membranes and tendons [1][2][3][4], polymer blends and composites [5][6][7][8][9][10][11], the liquid-gas and liquid-solid interfaces [12,13], and suspended membranes [14]. The present trend toward higher scanning speeds [15,16] and higher resolution mapping of mechanical properties [17,18] necessitates more rapid tip motion and therefore a larger viscous contribution to the tip-sample force. Viscous force may even dominate over elastic force with soft materials and a proper characterization of the material must therefore rely on a dynamic measurement that distinguishes viscous from elastic force. Here we employ such a method to demonstrate how traditional models for extracting material properties from AFM data fail to explain measurements on soft materials. We introduce the most simple form of an alternative model that does explain the data. We show how our model parameters relate to physically meaningful numbers that characterize mechanical response of the viscoelastic surface.
The AFM tip oscillating in and out of contact with the sample is a nanometer-scale example of a broader class of dynamical systems known as impact oscillators. These are often modeled with piece-wise non-smooth impact forces that produce interesting bifurcations [19]. The oscillation in dynamic AFM typically has a frequency close to a cantilever resonance with high quality factor Q, where the inertial and linear restoring forces of the cantilever body dominate the system dynamics. Nev- * haviland@kth.se ertheless, much smaller impact forces that weakly perturb the harmonic motion of the cantilever can be obtained with the help of lockin measurement techniques. We use a well established multifrequency lockin method where many Fourier components of this nonlinear perturbation are detected as intermodulation products near resonance [20], providing a great deal of information about the impact forces.
When the rigid AFM tip impacts the soft sample it experiences different types of force: elastic force associated with strain in the contact volume and curvature of the contacting interface, and viscous force associated with the rate of change of strain and curvature. These sample deformations occur at the nanometer scale, where the surface-to-volume ratio is much larger than that of macroscopic contact mechanics. We therefore expect interfacial forces to play an important role in AFM. A viscoelastic model of AFM contact forces must look beyond bulk rheology to also include the rheology of the interface [21].
The traditional approach to quantitative AFM is based on contact mechanics, where tip-sample force is considered to result from bulk elasticity. The Hertz model [22] gives a repulsive contact force. Assuming an axially symmetric rigid tip, one can parameterize the tip profile in terms of a power m ∈ [1,2] and a length scale [23], to express the contact force as a function of tip indentation z 0 − z, where z is the position of the tip and z 0 is the equilibrium position of the surface (see fig. 1a).
Here E eff is the effective modulus. Typically one does not know , m, or z 0 and they should therefore be treated as free parameters when fitting the measured force curve F T S (z) to extract E eff . With stiff samples one should not assume that the tip is rigid, in which case several more parameters are required to actually get the material modulus from E eff . Whatever the case, this class of models gives only constant forces that are not explained by dissipation. Such models can not say anything about the sample viscosity, as it gives rise to a velocity-dependent force.
A viscoelastic version of the Hertz model has been studied in the context of two-body collisions, where impact forces can be calculated in the center-of-mass reference frame [24]. With the center-of-mass in the laboratory frame, defined as the inertial reference frame where the entire sample is at rest (i.e. sample mass cantilever mass), the viscoelastic Hertzian model reduces to a Kelvin-Voigt expression for the tip-sample force where dissipation is introduced via a viscous damping coefficient η(z) that depends on tip position. One can in principle extract η(z) from a dynamic AFM measurement [25], but relating it to sample viscosity would involve a complicated model requiring knowledge of the tip geometry. Other models of contact viscoelasticity use a creep-compliance picture, where the elastic modulus is time-dependent and force is determined by integrating over deformation history [1][2][3]. One can also use finite element methods with linear [26] or nonlinear force-displacement relations that account for an attractive tip-sample force [27]. Independent of these bulk viscoelastic models, tipsample force in AFM may also result from interfacial energy or surface tension γ. The work of adhesion is the change in total surface energy after contact, where the subscripts T and S refer to tip and sample respectively. Typically W > 0 giving attractive force upon contact. The role of surface energy in small mechanical contacts was originally discussed by Johnson, Kendal and Roberts (JKR) [28]. The JKR model has been shown to break down when the contact radius of curvature R is small in comparison with the elastocapillary length, L = Υ/E, where Υ is the surface stress and E is the Young's modulus of the sample bulk [29,30]. For a soft material E ∼ 3 MPa forming a contact with relatively low interfacial energy γ TS ∼ 30 mN/m and no additional surface stress (in which case Υ = γ TS ), we find L ∼ 10 nm, the typical radius of an AFM tip. Thus, an AFM tip contacting a soft material should resemble a liquid-like sample wetting and forming a meniscus around the tip, as opposed to the tip compressing an elastic solid. Indeed, W is called the spreading parameter in the context of wetting phenomena if the sample corresponds to the liquid and the tip to the solid substrate, and W > 0 is called total wetting.
The discussion above makes clear that a threedimensional continuum model of AFM impact forces is a) Coordinates and piece-wise-linear (PWL) interaction. The cantilever deflection d = z − h is measured by the AFM detector, where z(t) is the instantaneous tip position and the constant h is the equilibrium (zero force) tip position. z0 is the equilibrium position of the sample surface.
(a) The traditional view of contact forces in AFM has the tip and surface moving together when they are in contact, z = zs, and interaction force is considered to be a function of the surface indentation (z0 − z). (b) The moving surface model treats the surface position zs(t) as an independent dynamic variable. The model introduces elastic and viscous forces that depend on surface deflection, ds = zs − z0 and velocityḋs, and an interaction force that depends on the separation s = z − zs andṡ . Forces are balanced in the inertial reference frame where the cantilever has a fixed working distance to the sample w = h − z0. (c) The piece-wise linear (PWL) interaction force plotted together with a calculated DMT model frequently used in AFM [31]. The parameters for the DMT model are: reduced modulus E * = 100 MPa, Hamaker constant H = 8 × 10 −20 J, tip radius R = 10 nm, inter-molecular spacing a0 = 0.16 nm. Parameters for the PWL model are: F ad = 5 nN, kv = 0.9 N/m. The range of separation shown is that typically covered by the cantilever oscillation for dynamic AFM on soft materials.
difficult to formulate if one includes both bulk and interfacial phenomena, and both elastic and viscous contributions. Knowledge of the tip and sample geometry is required, something that is difficult to determine either before or after an AFM measurement. Even if a realistic model could be formulated, determining its many free parameters from AFM data appears to be a hopeless task because the AFM measurement gives the dynamics of only one degree of freedom, namely the vertical position of the tip z(t). It is therefore well-motivated to formulate reduced models which approximate the tip-sample collision dynamics. The basic assumption of nearly all models in AFM is that the tip-sample interaction force can be expressed as a function of two dynamic variables, We argue that this basic assumption is incorrect for soft materials as it neglects the fact that the sample itself also has dynamics, which clearly influence the tip-sample interaction. Our approach is to reduce the complicated three-dimensional dynamics of the sample to a simplified one-dimensional dynamics in terms of a single degree of freedom z s (t), a generalized or spatially-averaged vertical position of the surface in the laboratory frame [see Fig. 1b)]. If we neglect the inertia of the sample contained in the very small interaction volume, our interaction force depends on this one additional 'hidden' dynamic variable This type of model was introduced by Cantrell and Cantrell to account for externally forced oscillations of the sample in the context of ultrasonic AFM. [32] Below we describe a simplified version of the moving surface model. A more complicated variation was presented and compared to experiments on soft materials in our previous publication [33]. The simpler model has analytic solutions in special cases which give physical intuition. We simulate the simplified model and introduce numerical optimization to find the model parameters that best fit experimental data. We argue that the traditional approach to quantitative AFM, where model parameters are bulk material properties such as elastic modulus and viscosity, is not appropriate for soft materials because a physically correct and complete model would involve far too many free parameters and uncontrolled assumptions. Rather, we demonstrate that a dimensionless ratio of our model parameters, called the impact-rheology factor R, gives a physically relevant and useful quantity for nanoscale mechanical characterization of viscoelastic surfaces. Our measurements show that the impact-rheology factor is independent of oscillation frequency and details of the tip geometry.

II. THE MOVING SURFACE MODEL
All coordinates of the moving surface model are defined in Fig. 1b. The essential difference to the traditional approach (Fig. 1a) is that we introduce the instantaneous position of the surface z s (t) as an independent dynamic variable and we express the interaction force as a function of the tip-surface separation s = z−z s . Unlike traditional models, we do not require that the tip and surface move together (i.e. z = z s ) when they are in contact.
One could imagine many different nonlinear models for the interaction force, but we make a sweeping simplification by linearizing the interaction when in contact. Figure 1c shows a piece-wise linear (PWL) force model described by the following equations, The interaction force is zero when the tip is out of contact, s > 0. When in contact, s ≤ 0, the adhesion force F ad turns on, correspond to a lowering of the total interfacial energy. Adhesion is counteracted by a repulsive force linear in the penetration −s, with force constant k v . We also include the possibility of a viscous force, linear inṡ =ḋ −ḋ s (h and z 0 are constant) associated with material flow as the tip penetrates the sample. The model preserves one essential feature of the interaction which is well known in AFM -the large force gradient d dz F TS (s,ṡ) localized near the point of contact s = 0. This rapid change of force is responsible for the jump-to-contact and pull-off hysteresis seen in nearly all quasi-static force curves, or measurements of d(h) for smallḣ, traditionally analyzed in quasi-static AFM. In dynamic AFM the amplitude of oscillation is typically much larger than the range of this localized interaction. We may therefore approximate this region of large interaction gradient as an adhesion force which instantly turns on and off when crossing the point of contact s = 0.
The interaction force couples the dynamics of the cantilever's flexural eigenmode to the dynamics of the viscoelastic surface, described by the following set of equations Here d s = z s − z 0 is the deflection of the surface from its equilibrium position z 0 and F drive (t) is the drive force. The free motion of the tip is given by eq. (7a) with F TS = 0, describing a driven damped harmonic oscillator. This model is valid in a narrow frequency band around a high quality factor resonance of the cantilever. The three mode parameters: stiffness k, resonant frequency ω 0 = 2πf 0 = k/m and dimensionless quality factor Q = ω 0 τ = η 2 /mk, are independently determined using a calibration procedure (see Methods).
The free motion of the surface is described by eq. (7b) with F TS = 0. Note that eq. (7b) does not have an inertial term corresponding to force arising from acceleration of the sample mass. Neglecting sample mass is valid when deformation occurs only in a local volume around the tip, however inertial forces may arise if the cantilever excites standing surface waves [33]. Thus the model effectively puts the center-of-mass at the equilibrium position of the tip, which we treat as fixed in the laboratory frame. We therefore neglect the very small amplitude base motion needed to inertially actuate the high Q resonance [25], as well as any forces arising from rapid changes of the probe height due to overly active surface-tracking feedback.
We can develop intuition for the sample dynamics by considering what happens when the tip is held rigidly fixed in the laboratory frame (i.e. not connected to a flexible cantilever) at the height of the unperturbed surface. As the tip just touches the surface from above [see Fig. 2(a)], the case s = 0 in eq. (6) together with the condition z = z 0 , or equivalently s = −d s in eq. (7b), gives Solving this equation we find that upon contact, the adhesion force lifts the surface forming a meniscus with asymptotic height, in a characteristic contact formation time, Similarly, when the tip just separates from the lifted surface, eq. (6) for the case s > 0 and eq. (7b) give, describing free relaxation of the surface to its equilibrium position in a characteristic time The contact formation and free relaxation dynamics are depicted in Fig. 2. We may also define the time constant associated with tip penetration into the sample. The model behaves as we would intuitively expect for an interaction consisting of two opposing forces: attractive surface force resulting from minimization of interfacial energy, and repulsive volumetric force resulting from compressive stress in the bulk. The stiffness parameters k s and k v , include both surface and bulk forces with the relative contribution depending on the size of the deformation in relation to the elastocapillary length. Nevertheless, a liquid-like interaction is described by K ≡ k s /k v 1, in which case the meniscus lifts to a maximum height, A solid-like interaction is described by the opposite limit K 1, in which case δ δ 0 as capillary adhesion is counteracted by compressive stress in the contact volume.
The time constants τ c , τ s and τ v represent ratios of viscous to elastic force constants of the model. The dynamics of both the tip and the sample depend on these characteristic times scales and their relation to the time spent in and out of contact during a single oscillation cycle, where the latter is set by the frequency of cantilever oscillation, the amplitude of motion and working distance to the surface w = h − z 0 . When the contact formation time τ c is larger than the time spent in contact, the surface can not fully deform to achieve contact equilibrium. When the free relaxation time τ s is large compared to time spent out of contact, the surface can not relax to its equilibrium position before the next tap of the tip. Repeated taps result in a steady-state dynamics characterized by a time-average up-lifted or indented position of the surface. In the following section we demonstrate this behavior of the model and correlate it with data from experiments. We show that this simple model explains the experimental data remarkably well for a variety of different soft samples. Fitting the model to experimental data we extract the parameters which characterize the material and its surface. We then comment on how these parameters may be useful for a general, quantitative characterization of viscoelastic samples with AFM.

III. COMPARISON WITH EXPERIMENT
We have validated the model by fitting it to data collected on several different samples. Figure 3 shows scanned images of four samples discussed below. To display the results of both experiment and theory, we show dynamic force quadrature curves using a technique called Intermodulation AFM (ImAFM) [20]. Force quadratures have been introduced in previous publications [34,35] but they are rather unfamiliar to the majority of AFM scien- tists so we briefly describe them here. Force quadrature curves do not display the instantaneous force on the tip as a function of tip position F TS (z) (i.e. 'force curves'). Rather, force quadratures represent integrals of force over a single oscillation cycle, making them the natural force curves of dynamic AFM.
The reader familiar with dynamic mechanical analysis (DMA) will recognize a similarity between force quadratures and typical DMA experiment. In DMA, the strain in a bulk material is measured while a sinusoidal stress of fixed amplitude is applied. The in-phase and quadrature components of the strain are then monitored as a function of the frequency of excitation (or the temperature of the sample) to characterize the complex modulus E * = E + iE of the material under test. To determine the AFM force quadratures we instead monitor the in-phase and quadrature response of the cantilever, subject to a sinusoidal excitation at fixed frequency but with slowly varying amplitude. This allows us to study how the nonlinear features associated with the impact of the tip on the viscoelastic sample change with amplitude. To explain these features we must look beyond the bulk concepts of storage and loss moduli to also include surface forces and adhesion.
The high quality factor of the cantilever oscillations means that the stored energy in the cantilever oscillation is much larger than the tip-sample interaction potential or energy lost to dissipation in the cycle. The tip motion is therefore well-described by harmonic oscillation at a fixed 'carrier' frequencyω ω 0 . The second drive tone of ImAFM gives rise to a slowly-modulated amplitude and phase of the carrier. From the intermodulation spectrum of these two drive tones, which is concentrated near resonance, we extract the Fourier coefficients of the tip motion and tip-sample force, at this carrier frequency. Rotating the force coefficients by the motion phase, we project out the Fourier coefficients of force which are in phase with the motion F I , and that which are quadrature to the motion F Q [34,35]. Both are presented as functions of oscillation amplitude A.
F I represents a conservative force (i.e. energy recovered in the oscillation cycle), and F Q represents a dissipative force (i.e. energy lost in the cycle).
When the tip is tapping on the sample, large Fourier components of the tip-sample force also exist at frequencies far above resonance but well below the detection noise floor. Thus, we can not determine the instantaneous force at any given time in the cycle. However, using the multifrequency lockin technique embodied in ImAFM, we can determine the integrated force over single cycles with very good signal-to-noise ratio. The two curves F I (A) and F Q (A) are found by direct transformation of the intermodulation spectrum, requiring only the calibrated linear response function of the cantilever eigenmode. No assumptions are made regarding the specific nature of the tip-sample interaction. The transformation is quite efficient computationally, allowing for immediate examination of the force quadrature curves at individual pixels, while scanning. Figure 4 shows force quadrature curves taken at two different points on a blend of polystyrene (PS) and polycaprolactone (PCL), shown in Fig. 3(a). The force quadratures on the softer PCL (nominal bulk E ∼300 MPa) show larger magnitude of the dissipative force F Q and a conservative force which is dominantly attractive, F I > 0, even at the highest amplitude. To fit the theory to the experimental data, we simulate the model dynamics by numerical integration of the equations of motion eq. (7a) and eq. (7b), adjusting the parameters of the model to find a best fit. Details of this procedure are given in the Methods section.
On the soft PCL we find an excellent fit. The simulation accurately reproduces the magnitude, complex shape and hysteresis of both F I (A) and F Q (A). The simulated   We compare the behavior on the softer PCL with that on the stiffer PS domain (expected E ∼3 GPa). On PS we find that the conservative force quadrature F I becomes dominantly repulsive (F I < 0) at large amplitude, and we see a smaller magnitude of the dissipative force F Q . Simulations show that the tip penetrates as much as 10 nm below the surface, but the surface lifts only slightly and quickly relaxes to its equilibrium position before the next impact. For increasing A the simula-tion captures the shape and magnitude of both F I (A) and F Q (A). However for decreasing A we find that the simulation does not reproduce the hysteresis observed in the experiment. A different interaction function F TS (s) could improve the quality of the fit, but at the expense of introducing additional model parameters. In our experience, the simplified model presented herein often does not capture hysteresis observed together with F I < 0, but it nearly always captures hysteresis when F I > 0. Figure 5 shows an additional example of model fitting at two locations on a sample consisting of PDMS mixed with hydrophobic silica nanoparticles, shown in Fig. 3b. This soft sample has faster free relaxation (smaller τ s ) such that no hysteresis is seen in the force quadratures. The model reveals large amplitude surface motion when tapping on the PDMS matrix. When tapping on a region with dispersed nanoparticles, the model shows a much slower surface and very low amplitude surface motion, but equally deep penetration.
We also investigate how the fitted model parameters depend on a change of tip and cantilever oscillation frequency. Using three standard cantilevers with calibrated parameters given in Table I, we study a wellknown soft material consisting of micron size domains of LDPE (nominal bulk E ∼100 MPa) in a matrix of PS [37], as shown in Fig. 3c . The cantilever resonant   frequencies and mode stiffnesses span one order of magnitude. In order to make a reasonable comparison we adjust the excitation to keep the amplitude of free motion A free such that the stored energy in the free oscillation E free = 1 2 kA 2 free is approximately the same in each measurement. Table I also shows the best-fit parameters of the moving surface model.
For each cantilever we also study how the force quadrature curves change as we vary the working distance, w = h − z 0 , or static probe height above the relaxed surface. We control w by simply changing the scanning feedback set-point at regular intervals along the slow scan axis [see Fig. 3(c)]. Figure 6 shows an example of the measured and simulated force quadratures for the Tap525 cantilever at three different probe heights. All simulation curves shown in the figure use the same parameters given in Table I, changing only the working distance w which is given in Fig. 6. Equivalent plots for the NSC15 and AC55 cantilevers are shown in the supplemental material Figs. S1 and S2 respectively. We see that the model reproduces the shape and offset of the force quadrature curves rather well.
To further test the model we also analyze a nanoparticle membrane sample shown in Fig. 3(d). This sample consists of a monolayer of Au nanoparticles, bound together by organic ligands. The membrane is suspended over a circular hole 2.5µm in diameter, forming a ultrathin drum head [38,39]. Previous studies have shown large apparent stiffness, corresponding to a large effective bulk modulus of the membrane material [40]. Interpretation of the model parameter k v in terms of bulk compression is not obvious with these samples that have thickness 7nm. However, adhesion gives rise to very large curvature of upper and lower surfaces of the membrane as it covers the tip with very small radius. The Laplace pressure change across both of these interfaces adds, giving the force described by k v which acts in opposition to adhesion. The stiffness constant k s then describes elastic tension in the membrane as the surface deviates from equilibrium. It is therefore not surprising that our simple model fits the data quite well. Figure 7 shows a comparison of experiment and simulation of the model for two different cantilevers. The qualitative shape of the force quadratures and the simulated dynamics of the surface are the same for both cantilevers, despite the difference in resonance frequency of AC55 (∼ 2 MHz) and Tap300 (∼ 300 kHz). Low-density polyethelene (LDPE) with Tap525 cantilever. Experimental (left) and simulated (right) force quadrature curves, each offset vertically for clarity with dashed lines corresponding to zero force. The experimental curves result from analysis of data at pixels marked with an × of the corresponding color in Fig. 3(c). The working distance is changed in the experiment by adjusting the amplitude-feedback setpoint to the value given in left panels (% of free amplitude). For simulation the working distance w shown in the right panels is found by numerical optimization with all other parameters of the model fixed at the values given in Table I ). Ratios of the model parameters (bottom box) are: the relaxation times τs = ηs/ks and τv = ηv/kv, contact formation time τc = (ηs + ηv)/(ks + kv), stiffness ratio K = ks/kv and impact-rheology factor R = τs/τv. Although the stiffness and damping parameters of the model show considerable variation between probes with different frequency and tip geometry, the time constants show much less variation. The impact-rheology factor R shows little variation between probes.

IV. DISCUSSION
An obvious criticism of the model presented herein is that it is ad hoc, or not based on first principles. Nonlinear forces arising from sample deformation, the analysis of which forms the basis of traditional AFM nanomechanics, are avoided by linear approximation. The model keeps only one nonlinearity to capture the sudden impact and release of the tip and sample. Reduction of the sample dynamics to one effective degree of freedom d s (t) is a sweeping simplification which neglects the fact that the surface deformation has a transverse profile.
We emphasize again that our intent is to describe the measured dynamics of the cantilever with a minimum number of free parameters and additional degrees of freedom. A physically more realistic short-range interaction force which smoothly interpolates between the two linear regions could be constructed, but at the expense of additional parameters which are not necessary to explain the data. We consider this to be the minimal model necessary for explaining dynamic AFM on viscoelastic materials, and in our experience it is often sufficient.
Starting from first principles, even very idealized models would involve many parameters: tip radius and sidewall angle -assuming an axial-symmetric tip; elastic modulus, Poisson ratio and viscosity -assuming homogeneous half-space; surface energies and local curvature of the unperturbed surface. All parameters not independently determined should be treated as free when fitting AFM data. The fit procedure here uses the multifrequency ImAFM data at some 40 frequencies near resonance. While the quadrature data (80 values) do have different weight (they are measured with different signalto-noise ratio) each represents an independent observable with information about the tip-sample interaction. In stark contrast to other quantitative methods in AFM, the large number of observables in relation to the number of free model parameters makes physical interpretation of the fitted parameter values meaningful.
Examining Table I, we recall that the cantilever parameters f 0 , Q and k are independently calibrated (see Methods), and the working distance w can be determined by inspection (amplitude at the initial sharp rise in F I ). Thus the 5 parameters F ad , k s , η s , k v and η v , are free parameters. We expect that all of these will depend on the detailed shape of the tip. It is therefore not so surprising that the table shows considerable variation of these parameters between the different probes. We expect that the relaxation time τ s = ηs ks is independent of tip shape as it involves the free relaxation of the sample. Furthermore τ v = ηv kv might, to first order, be sensitive to tip shape as a blunt tip would result in both larger viscous and larger elastic force, in comparison with a sharper tip. Indeed, we do observe little variation in the values of τ s(v) for the two cantilevers with lower frequency, but there is a reduction of both τ s and τ v for the higher frequency cantilever.
It is interesting to note in is apparently independent of the AFM probe used in the experiment. We call R the impact-rheology factor [41], a dimensionless number formed from all four force constants of our PWL viscoelastic model. Our measurements suggests that this factor may be a good quantity for physical characterization of surface viscoelasticity. In order to quantify measurement uncertainty in R, one should study how various sources of noise in the experiment propagate to the model parameters in the fitting procedure. Further experimental studies and more detailed model analysis are required to study uncertainties, and verify if the impact-rheology factor is indeed a reliable measurement of an intrinsic property of the sample interface.
While this moving surface model is minimal with regard to the number of parameters, it nevertheless has very complicated dynamics. General statements about its behavior in different parameter regimes are therefore difficult to formulate. The impact-rheology factor represents the relative strength of two viscoelastic processes: large R means that the sample's free relaxation is slow in comparison to the time needed for the penetrating tip to locally deform the sample. Conversely, small R means that tip penetration is slow in comparison with free relaxation. However, the magnitude of R alone is insufficient to predict the shape of the force quadrature curves, which also depend on the adhesion force, working distance, and frequency of oscillation. Hysteresis is associated with ω 0 τ s 1, large surface lifting with small k s or large δ 0 , and deep penetration with small k v .

V. CONCLUSIONS
We presented a minimal model for dynamic AFM on viscoelastic materials, accounting for both surface and bulk forces. The model is simple in that it is linear where possible, but its dynamics is complex due to the nonlinearity describing the sudden impact and release of the oscillating tip and surface. Unlike traditional nanomechanical analysis of AFM, our model takes into account the viscoelastic dynamics of both the penetrating tip and free sample. The sample dynamics was shown to be quite significant on soft materials, where simulations revealed large amplitude surface motion. We validated the model by showing excellent agreement with experimental data on a variety of samples. By fitting the model to the data we extracted viscous and elastic force constants for the surface and bulk. Our analysis indicated that the impact-rheology factor, formed from the dimensionless ratio of these constants, is independent of tip shape and cantilever resonance frequency. This simple model describes numerous measured force quadrature curves with complex and differing shape, instillings confidence it captures the essential physics of dynamic AFM on soft materials.
Appendix: Methods

Measurement details
Each measurement (one pixel of the scan) requires 2 ms, corresponding to a frequency spacing in the comb of intermodulation products of 500 Hz, also the measurement bandwidth. To enhance the signal-to-noise ratio, we average the measured intermodulation spectra over several neighboring pixels which show the same type of response, thereby using spatial correlation to reduce the noise. This averaging results in smoother experimental force quadrature curves without reducing the scan speed.

Equations of motion
For numerical simulation of the tip and surface motion it is convenient to introduce a non-dimensional time u = ω 0 t, where f 0 = ω 0 /2π is the resonant frequency of the cantilever in Hz. With this scaling the relaxation times become dimensionless parameters, u s = ω 0 τ s and u v = ω 0 τ v . Equations (7a) and (7b) then describe a threedimensional dynamical system with state variables d , d and d s . For the case s < 0 this system reads, Since the right hand side of eq. (A.1b) depends on s = d − d s at each time step, we must first solve eq. (A.1c) for d s and then substitute it in to the right hand side of eq. (A.1b). In terms of the impact-rheology factor R = us uv and the stiffness ratio K = k s /k v we can write eq. (A.1b) as: 2) The prefactor 1 1+RK in front of eq. (A.2) describes the strength of the dynamic coupling between tip and surface, giving a measure of how much surface motion one can expect.

Numerical integration and optimization
Numerical integration of the model is performed with the package CVODE, part of the SUNDIALS suite of nonlinear solvers [42]. This integrator has adaptive timestepping. Because the dynamical system is coded in C, simulation time is reduced by a factor of 100 in comparison with coding in higher-level languages such as MAT-LAB or Python. This speed up is of great importance as many integrations are required when iterating to find the optimal parameters. When performing the numerical optimization, the initial conditions for d and d are reset to their measured values for each iteration, whereas the initial condition for d s is set to zero on the first iteration, and then estimated from the previous integration for all subsequent iterations.
We use the Python library scipy.optimize.leastsq, which is a wrapper for MINPACK's implementation of the Levenberg-Marquardt algorithm. The algorithm performs a least-square minimization of the array of residuals r = {real(F exp −F sim ), imag(F exp −F sim )} where theF (ω) are experimental and simulated complex force amplitudes at 40 frequencies. The method finds a local minimum, starting from good initial parameter values determined by trial and error.

Background force compensation
When the cantilever is oscillating above a surface we often observe significant background forces, not resulting from tip-sample interaction but rather acting over the entire cantilever body. The origin of these background forces might be e.g. long-range electrostatic force or hydrodynamic squeeze-film damping. In order to deduce the tip-sample force we need to remove this background interaction. The procedure we use for removing any linear background force is described in a previous publication [43].

Calibration
Cantilever parameters are determined by the noninvasive thermal calibration method described by Higgins et al. [44], which combines the fluctuation dissipation theorem with Sader's method based on analysis of hydrodynamic damping [45]. With this approach one can extract the three parameters of the cantilever transfer function, k, f 0 and Q, as well as the inverse responsivity of the optical detector used to measure cantilever deflection α −1 [nm/V], all from one measurement of the thermal Brownian motion of the cantilever near resonance. The method is encapsulated in the recently launched Global Calibration Initiative (GCI) [46], where a thermal noise measurement of f 0 and Q can be used to get k, based on a single hydrodynamic constant determined from averaging over the measurements of may users on the same type of cantilever. We take our measured f 0 and Q and use the GCI to determine k Sader , allowing us to then determine α −1 . Note that error in the calibration of k and α −1 result in a re-scaling of the force and amplitude axes respectively, which does not change the general shape of the force quadrature curves.