A fast first-principles approach to model atomic force microscopy on soft, adhesive, and viscoelastic surfaces

Quantitative atomic force microscopy (AFM) on soft polymers remains challenging due to the lack of easy-to-use computational models that accurately capture the physics of the interaction between the tip and sticky, viscoelastic samples. In this work, we enhance Attard’s continuum mechanics-based model, arguably the most rigorous contact model for adhesive viscoelastic samples, via three key enabling strategies. First, the original model’s formalism is rearranged to enable a fast and explicit solution of the model’s ordinary differential equations (ODEs). Second, the deformed surface is reconstructed using a complete set of optimized orthogonal basis functions as opposed to Attard’s original, computationally expensive radial discretization. Third, the model’s governing ODEs are solved using a multi-step numerical method to further stabilize the solution when using for soft and sticky samples. Implementing these enhancements, enhanced Attard’s model (EAM) is more stable, 3+ orders of magnitude faster, and equally accurate when compared to the original model. These facilitate EAM’s inclusion into simulations of various AFM operating modes. We demonstrate EAM based simulations of quasi-static force spectroscopy and amplitude modulation AFM approach curves on soft sticky polymer surfaces. On a typical desktop computer, simulation of an amplitude modulation approach curve with EAM takes less than a minute as compared to ≈15 h by the original Attard’s model. We expect EAM to be of interest to the AFM community because it facilitates the inclusion of rigorous models of tip-sample contact in simulations on polymer samples. EAM is available as part of the VEDA set of simulation tools deployed on nanoHUB.org cyber-infrastructure.


Introduction
Successful quantitative estimation of local sample properties from AFM data, whether in the form of forcedistance curves or maps of dynamic AFM observables, hinges on how well the utilized tip-sample contact model reflects the actual physics of the interaction. Soft, adhesive, and viscoelastic materials represent an exciting area of growth of new materials whether for biomedical, critical node photolithography, flexible electronics, or in a wide range of consumer products [1][2][3][4]. Quantitative AFM on such materials requires accurate, versatile, and fast contact mechanics models to link experimentally observed hysteresis, frequency dependence, and creep to the material's viscoelastic and adhesive properties [5][6][7][8][9][10][11].
However, the currently utilized models to calculate the tip-surface interactions on soft viscoelastic and adhesive materials are limited in various ways. Either they do not utilize first-principles methods, make ad hoc self-inconsistent assumptions, neglect relevant physical parameters of the contact phenomenon, employ computationally inefficient and slow procedures, and/or lack a simple and direct path to integrate with AFM operational modes [5][6][7][11][12][13][14][15]. Furthermore, almost all viscoelastic contact mechanics models used by the AFM community calculate the tip-surface interaction forces assuming that the force is only a function of position and velocity of the tip: surface particularly on soft adhesive materials [17]. The long-range attractive forces between interacting surfaces cause these experimentally observed surface instabilities to occur prior to and after the actual tip-surface contact and lead to energy dissipation in the form of heat [18,19]. In many viscoelastic contact models commonly used by AFM community, the creep compliance function is simply inserted into Hertz [20] or Johnson-Kendall-Roberts (JKR) [21] theories formalism. However, the static free energy minimization procedure implemented in Hertz or JKR theories does not hold in the case of these ad hoc viscoelastic models [22]. Interested readers are referred to [11] for further details on the deficiencies associated with the ad hoc viscoelastic models and their non-reliable predictions. On the other side, the continuum mechanics-based viscoelasticity models which are soundly based on the established elasticity and viscoelasticity equations, do not usually result in analytical closed-form solutions, and consequently lead to higher computational expenses. Ting [23] proposed one of the oldest viscoelastic models in this category which despite its high reliability, does not account for the surface forces. Therefore, a fast first-principles contact mechanics model which enables rigorous prediction and interpretation of AFM images on soft, adhesive, and viscoelastic polymers is required to address these shortcomings. Attard et al [12,13,24] introduced a novel approach that provides versatility in modeling tipsample interactions by including different linear viscoelastic constitutive models as well as arbitrary surface forces within a continuum mechanics framework. The model can predict the surface deformation as the tip interacts with adhesive, viscoelastic samples by applying the correspondence principle on the Boussinesq solution for a semi-infinite half-space [25] and calculates the interaction force by implementing radial and temporal discretization. Attard's original implementation, which is akin to a boundary element method, computes the distributed pressure and the associated surface deformation over the contact region at each radial node and each timestep. However, Attard's implementation does not lead to a closed-form solution to calculate the tip-surface interaction forces. The implementation is computationally expensive due to the intrinsic temporal and radial discretization in addition to the iterative loops needed to solve the model's governing differential equations. Moreover, Attard's original implementation leads to computational instabilities when the surface undergoes abrupt bifurcations/instabilities under the action of the tip, especially on soft and sticky surfaces.
In this work, we use three enabling strategies to significantly enhance the computational part of Attard's model and render it into a faster and more robust method. First, we rearranged the model formalism so that its governing differential equations become explicitly solvable. Next, we approximated the deformed surface profile using a complete set of optimized orthogonal basis functions to replace the radial discretization in the original model. Finally, instead of the forward Euler method employed by Attard, we used the two-step Adams-Bashforth method to solve the model's ordinary differential equations. Implementation of these strategies accelerates the model's computations by at least three orders of magnitude. Furthermore, the enhancements improve the computational robustness of the model when predicting tip-surface interaction on highly soft, viscoelastic, adhesive surfaces. We refer to the combination of Attard's model and the three identified improvements as the enhanced Attard's model (EAM). Leveraging EAM's capabilities within an AFM framework leads to a better understanding of the complex phenomena that occur during tip-surface interaction, especially on soft sticky samples.
To demonstrate the utility of EAM in the AFM context, we performed simulations of a quasi-static AFM force curve and a dynamic approach/retract curve on a polymeric surface. The resulting simulation tools based on EAM are now included as part of the VEDA suite of tools [26] on the cloud computing cyber-infrastructure of nanoHUB.org and are thus easily accessible to the AFM community worldwide. We expect that the utilization of EAM simulations on the nanoHUB to provide a fast and accessible path to gain further insight into the complex tip-surface interaction phenomena on soft viscoelastic sticky polymers samples.

Attard's original model and implementation
Attard's model and computational implementation provide a framework to compute tip-sample force and surface deformation history induced by a prescribed time-dependent motion of a rigid, axisymmetric AFM tip on an adhesive, viscoelastic surface [12,13]. in this section, we briefly summarize key elements of Attard's original model and its implementation [12,13,22].
The elastic equation for a semi-infinite half-space predicts the surface deformation due to a distributed applied force as follows [27]: where r is the position vector to the point of interest on the undeformed flat surface, and 'r' is the radial coordinate of the position vector such that = r 0 corresponds to the position of a point on the undeformed sample surface that lies directly beneath the axisymmetric AFM tip apex. u r t , ( ) is the instantaneous out-ofplane surface deformation at r, and p r t , ( )defines the distributed pressure applied on the surface with elastic modulus E and Poisson ratio n. When an axisymmetric rigid tip contacts a surface, p r t , ( )arises from the local instantaneous tip-surface interaction forces which their magnitude depends on h r t , , ( ) the tip-deformed surface distance at time t. These parameters for an interacting tip-surface ensemble are shown in figure 1. With a parabolic approximation of the rigid tip profile, the geometry parameters can be related as: where, h 0 is the distance between the tip apex and the undeformed surface at time t (figure 1). h t 0 ( ) defines the prescribed tip trajectory as a function of time.
To extend equation (1) to tip-surface contact on linear viscoelastic materials, the surface deformation history is accounted for by implementing the correspondence principle on the linear elastic solution as follows: where, = p h r t dp h r t dt , , ( ( )) ( ( ))  and t 0 is the initial time instant of the computations. Tip position at t 0 is assumed to be far enough from the sample so that the tip-surface interaction is negligible and the surface is flat and stationary » u r t , 0. 0 ( ( ) ) As seen in equation (3), the deformation of viscoelastic surfaces in a continuum mechanics perspective requires a time integral over its preceding deformation history.
To model the surface viscoelasticity behavior, Attard employs the standard linear solid (SLS) constitutive relationship whose creep compliance function is defined as: where, E 0 and ¥ E are the short-and long-term modulus of the surface, E t ( ) is the effective instantaneous viscoelastic modulus, and τ is the characteristic creep(retardation) time. When τ is very short or long comparing with the tip-sample interaction time, the SLS model reduces to the Kelvin-Voigt or Maxwell models, respectively. However, when τ is chosen to be comparable with the tip-sample interaction time, the SLS model can capture more complex material behavior that cannot be predicted by the Kelvin-Voigt or Maxwell models. The tip-surface pressure between any point on the surface with the tip is computed using a Lennard-Jones (LJ) pressure model: where H is the Hamaker constant and z 0 is the equilibrium distance. These two physical parameters, H and z 0 define the work of adhesion which is the reversible thermodynamic work to separate two interfaces from equilibrium state to distance of infinity:  where, g 1 and g 2 are surface energies of the two bodies and g 12 is the interfacial energy. Attard's model calculates the interaction between each of the introduced radial elements on the surface with their associated counterparts on the tip using the same method invoked in the derivation of the well-known Derjaguin approximation [28,29]. The model can accommodate alternative and/or extra surface force models or more sophisticated linear viscoelasticity constitutive models such as the generalized Maxwell model (also known as Prony series) or power-law rheology to capture more complex sample behavior.
Substitution of the SLS model's creep compliance function, equation (4), into equation (3) and its subsequent differentiation with respect to time cast the equations in the form of differential equations and remove the time-convolution integral, as follows: Due to the axisymmetry of the problem, Attard further simplifies the equation by expressing the kernel of the integrals in equations (7) and (8) in terms of the complete elliptic integral of the first kind with r and s representing the radial coordinates of a point located at r and s respectively: where r d is the computational domain beyond which the tip-sample interaction is considered negligible. k r s , ( ) is defined as: where, K is the complete elliptical integral of the first kind. Since = k r s k s r , , ( ) ( )the k square matrix is symmetric.
To solve equation (10), Attard discretizes the computational domain < < r r 0 d ( ) into n A equal-in-size segments and assumes that all points in each segment have the same u r t , , ( )  u r t , , ( ) and h r t , .
( ) The parameters associated with each of these discretized segments are calculated based on its central point, hereinafter is called a node. Then, Attard rewrites equation (10) for each of these n A nodes and solves the resultant set of dependent ODEs for the deformation rate of each of the nodes [12,13,17]. Due to implicit u r t , term on the right side of equation (10), this set of ODEs cannot be treated and solved explicitly. Therefore, Attard solves them through iteratively adjusting a guessed u r t , ( )  until equality between both sides of equation (10) is established. Attard suggests using the forward Euler's method to predict u r t , ( ) for each timestep based on the calculated u r t , ( )  from the preceding timestep. Attard numerically calculates the integrals in equations (10) using introduced radial nodes as integration partitions. Simultaneously these nodes are used to reconstruct the deformed surface profile at each timestep and act like degrees of freedom (DOF) of the surface. This bifunctionality of the radial nodes is not computationally efficient since the number of radial discretization for sufficiently accurate integral calculations might not be necessarily the same as the required DOF to render the deformed surface. We will refer to the above-described computational method proposed by Attard as the 'iterative' approach.

Enhanced Attard's Model (EAM)
To improve the computational implementation of Attard's model, we first remove the slow and computationally intensive iterative solution of the model's differential equations. To do so, we expand the rate of pressure change in equation (10) as follows: Then, we substitute equation (13) into equation (10) and rearrange it so that all terms involving u r t , ( )  are on the left side of the equation, as below: The result is a set of nonlinear coupled ordinary differential equations with explicit time-dependent forcing through the prescribed tip motion term, h t , 0 ( )  on the right-hand side of equation (14). Hereafter, this method is referred to as the 'explicit' approach in this article. Equation (14) can be solved by introducing radial and temporal discretization like the iterative approach. Using the explicit approach solution significantly reduces the computational time in comparison with the iterative approach.
Furthermore, we propose a method to optimize the computations speed by separating the dual functionality of the radial nodes in the original contact model. To do so, we apply separation of temporal and radial variables on u r t , ( ) as follows: where a t i ( ) represent the time-dependent coefficients of the radial basis functions and f r i ( ) are a complete set of functions that approximate the surface deformation profile at each timestep. Substitution of equation (15) into equation (14) results into a set of differential equations as follows: Thus, instead of solving a set of ODEs in equation (14) each associated with a radial element, we solve equation (16) for the time-dependent coefficients, a t , i ( ) that are each associated with a selected radial basis function. The resultant advantage is that the number of basis functions can now be adjusted and optimized independently of the number of radial integral partitions.
While almost any complete set of functions can be chosen for f r , i ( ) an accurate reconstruction of the deformed surface is achieved with smaller n b if the shape of selected functions resembles the expected deformed surface profiles during an interaction cycle. To further simplify equation (16), we choose f r i ( ) to be a complete orthogonal set of functions, for which we have: Since the problem is axisymmetric, we only consider the even terms of selected basis functions. Then, we multiply both sides of equation (16) with f r j ( ) and integrate over (0, r d ) to utilize the terms orthogonality and simplify the resultant equation: where, i and j are integer numbers between 0 to n . b r s j ( ) and z j are time-independent variables and hence, are not required to be calculated for each timestep. We rearrange equation (18) into a matrix representation and introduce temporal discretization as follows: To evaluate and optimize the performance of the proposed computational approach, we studied an alternative way of implementing the approximating basis functions. In this scenario, we express the tip-surface distance h r t , ( ) instead of u(r,t) in terms of a complete set of orthogonal radial basis functions, as follows: The implementation process for equation (25) is similar to the previously explained procedure.
The third important improvement made in EAM relates to computational stability. In cases where the selected model to calculate p h r t , ( ( )) switches between attractive and repulsive forces depending on the gap between interacting bodies, Attard's model computations may become unstable. The computational instability occurs when the surface undergoes rapid non-equilibrium movements during the interaction cycle, specifically during snap-in or -off surface processes for which the model predicts large u r t , .
where 'dt' is the infinitesimal time interval between timesteps. Using forward Euler's method, the projected of the proceeding timestep. Therefore, when u r t , q ( )  is very large, it may lead to discontinuities in the reconstructed deformed surface profile and resultant computational instabilities. This may happen due to surface deformation instabilities during interaction time such as snap-in and/or snap-off instances when the tip approaches or retracts off the surface, respectively, especially on highly soft and sticky surfaces. To alleviate this issue, we employed the general form of multistep Adams-Bashforth method to establish a more smoothly controlled link between u r t , ( ) at each timestep and u r t , ( )  of the proceeding timesteps: Figure 2. AFM force spectroscopy schematics and used parameters. In this mode, the input signal to Z-piezo prescribes the defined Z motion. The microcantilever moves toward the surface with a constant velocity, indents, and then retracts as it reaches the defined force setpoint.
where, as per definition the rate coefficients must satisfy this condition: å [30]. Moreover, we successfully investigated the use of other methods such as Runge-Kutta methods to further improve the computational efficiency. However, for simplicity, we use the above-described method in the simulations of this study. We name the described method as formulated by equations (17)- (19) combined with equation (27) as 'EAM'.

Embedding EAM into quasi-static and dynamic AFM force spectroscopy
EAM and Attard's original implementation both assume that the AFM tip motion relative to the undeformed sample surface is a known, prescribed function of time. In most AFM modes, however, while the AFM microcantilever's excitation or Z oscillation range may be prescribed, the tip's trajectory depends on the ensuing tip-sample interactions as well as on the microcantilever's mechanics and dynamics. Here we demonstrate how EAM can be embedded within a simple microcantilever mechanics model and within the Amplitude Modulation AFM (AM AFM) amplitude reduction and phase lag formulas to predict the force spectroscopy response on soft, adhesive, viscoelastic surfaces.
First, we embed EAM computations into quasi-static AFM force spectroscopy simulations. The quasi-static AFM, in which a constant velocity Z-piezo expansion-retraction above a sample is prescribed, is commonly used to study the surface mechanical properties [31][32][33]. During a single cycle of this mode, the microcantilever approaches the surface, snaps into the surface, indents into the sample until it reaches the defined setpoint, and then withdraws back to its initial state. The Z-piezo periodic motion, Z t , ( ) has a frequency that is much smaller than the fundamental resonance frequency of the microcantilever hence this method is called 'quasi-static' force spectroscopy. A schematic of an AFM force spectroscopy experiment is shown in figure 2.
The tip trajectory with respect to the unperturbed surface level, = + h t Z t q t , 0 ( ) ( ) ( ) is not known a priori in quasi-static AFM force spectroscopy even if Z(t) is a known prescribed motion. This is because the tip deflection, q t , ( ) does not follow the triangular trajectory of Z(t), but rather depends on the tip-sample interactions. Therefore, we need an algorithm that calculates q t ( ) for each timestep during the cycle based on the associated tip-sample interaction calculated by EAM. Since the Z frequency is far below the fundamental frequency of the microcantilever, the q t ( ) output of the designed algorithm for each time step needs to satisfy the relevant quasi-static solution: where k is the quasi-static bending stiffness of the microcantilever. The algorithm determines q t ( ) by iteratively adjusting an initial guess until the difference between two sides of equation (28) falls below the defined tolerance. Note that the quasi-static assumption in the algorithm implies that transient ringing of the microcantilever upon detachment from the surface cannot be captured by this approach.
Next, we summarize how the EAM can be embedded inside AM-AFM's amplitude and phase equations to predict amplitude and phase response and force curves during a Tapping mode approach curve on an adhesive viscoelastic surface. The algorithm to do so with Attard's model has been discussed in detail elsewhere [11], so here we mention that the approach involves using an energy balance based algorithm employing the amplitude where, A ratio and A are the amplitude setpoint and tip oscillation amplitude, respectively. Furthermore, the virial V ts ( ) and energy dissipation E ts ( )in a single oscillating cycle of AM-AFM are defined as: where q t ( )  is the tip velocity and T is the period of a single oscillation cycle. Therefore, we need an algorithm to find the specific Z distance for which the computed V ts and E ts using equations (31) and (32) satisfies equation (29) within a prescribed tolerance [11].

EAM verification
To verify EAM and assure the reliability of the developed code, we compare EAM predicted F-d curves for a prescribed triangular tip motion with the ones presented in Attard's original work [12]. In this set of simulations, the paraboloid tip has a radius of 10 μm which oscillates through a triangular wave between h 0 = 10 nm and h 0 = −10 nm. The employed material properties and interaction parameters are E , and τ are Hamaker constant, equilibrium distance, Poisson's ratio, short-term modulus, Long-term modulus, and creep time, respectively. The tip radius in these simulations is set to be 10 nm.

EAM performance
We simulate tip-surface interaction in amplitude-modulation AFM (AM-AFM) on a low-density polyethylene (LDPE) sample using various computational setups of EAM to evaluate, compare, and optimize their output computational accuracy. We assume tip oscillation is steady-state and sinusoidal in AM-AFM which is legitimate when the operation is done in air/vacuum [34]: where ω and f are the excitation frequency and phase lag relative to excitation force, respectively.
To demonstrate the achieved computational stability enhancement, we used the Adams-Bashforth method with = n 0 E (forward Euler's method) and = n 1 E to predict a single force-distance (F-d) history when a rigid tip oscillating through a sinusoidal wave interacts with an LDPE sample as characterized in table 1. We regulated the adhesion of the surface by decreasing z 0 in the LJ pressure equation (equation (5)) to envisage when computational instabilities arise. As illustrated in figure 4, the Adams-Bashforth method implementation stabilizes the computations in comparison with the forward Euler's method.
To quantify the convergence of EAM solutions, we use virial V ts ( ) and energy dissipation E ts ( )values of a single oscillation cycle of the microcantilever on an LDPE sample. V ts and E ts are commonly used parameters to quantify AM-AFM and represent the conservative and non-conservative part of the interaction energy during tip-surface contact [11]. Since the output computational error of 'iterative' and 'explicit' approaches can be minimized, we consider them as 'exact' solutions. EAM simulations converge to exact predictions if a sufficient number of basis functions (n b in equation (15)) are employed.
The LDPE and EAM parameters used for this analysis are listed in table 1 and table 2. The free amplitude, A , 0 is 62 nm and the amplitude setpoint is 0.8 for these simulations. EAM computational setups 'C' and 'E'/'G' are computed using equation (16) and equation (18), respectively, and computational setups 'D' and 'F' refer to employing Fourier expansion basis functions to approximate h r t , ( ) instead of u r t , ( ) as described in equation (25). All solutions are checked to be independent of the selected domain and temporal discretization values. Since the problem is assumed to be axisymmetric, only the even terms of the employed basis functions are considered. The number of the employed basis functions is in the range of 0 to 56. The convergence of the  EAM computational setups is evaluated by comparing to the exact solutions ('A' and 'B') [11][12][13] and the associated energy errors are expressed as a percentage. We considered the 1% error on both V ts and E ts parameters as the convergence threshold and we determined the convergence when this criterion is met constantly beyond a specific n .
b The number of the timesteps is identical for all simulations 10 4 ( ) and the number of radial discretization terms used for the 'exact' computational solution is 70.
The results illustrated in figures 5(a)-(b) show that V ts and E ts computed using EAM computational setups A-G in table 2 gradually converge to the associated exact values by increasing n .
b Figures 5(c)-(d) depicts the percentage error of each of these computational setups and horizontal black dashed lines show the 1% error criteria with respect to exact solutions. Figure 5(c) shows that the Fourier basis function's implementation leads to a more optimized V ts convergence when used to approximate u r t , ( ) (computational setups 'C' and 'E') rather  The repulsive phase of the force history converges faster than the attractive phase.
The time elapsed to compute a single tip-sample interaction cycle when using the computational setups listed in table 2 is compared in figure 7. The significant difference between iterative ('A') and other computational approaches depicts the considerable enhancement achieved through the applied strategies. Furthermore, figure 7 shows that the basis functions orthogonality implementation (computational setups 'E', 'F', and 'G') leads to an optimized computational time due to the more simplified formalism. The Bessel set of basis functions ('G') despite the acceptable convergence rate, is slower than Fourier ones due to the intrinsic computational complexity of the Bessel functions and their orthogonality implementation. Figure 7 illustrates that in terms of the computational time, EAM can outperform the 'explicit' approach if the number of employed basis functions in EAM is less than the number of the implemented radial discretization in the 'explicit' approach.
Considering both results in figures 5 and 7, amongst the evaluated EAM computational setups, the tipsample interaction phenomenon is most efficiently computed by EAM as formulated in equations (21)-(23), when we use: • Equation (15) to approximate the surface deformation, and, • The cosine terms of Fourier expansion as the basis functions (setup 'E'):  Table 3. The values employed for LJ pressure and surface viscoelasticity properties for polystyrene (PS) in the conducted simulations. The tip radius in these simulations is set to be 50 nm.

AFM Simulations using EAM
We applied the algorithm for the quasi-static solution of force spectroscopy (Eqn. (28)) to simulate an AFM force spectroscopy on an LDPE sample with viscoelasticity and LJ pressure parameters as specified in table 1 except t = 0.01 ms. In this simulation, the microcantilever equivalent stiffness is 28.1 N m −1 and the microcantilever base travels through a triangular wave with 1 Hz frequency and 50 nm amplitude.
The results shown in figure 8 illustrate the model's ability to capture the surface instabilities both when approaching and retracting the surface. During retraction, the adhered surface to the tip gradually peels away until the adhesive forces again dominate and the microcantilever deflects downward (negative force). This will lead to a hysteretic loop in which the long-range surface forces may significantly contribute to the resultant energy dissipation. The whole cycle computation took less than a minute to complete using a typical desktop computer as compared to »12 h of the original Attard's model.
Next, an AM-AFM approach/retract curve on a polystyrene (PS) surface with viscoelasticity and adhesive properties as specified in table 3 is simulated using the algorithm described in section 2 for AM-AFM approach curves. We calculated E , ts V , ts and f as a function of A ratio as shown in figure 9. Computing this AM-AFM approach curve in figure 9 using EAM with = n 25, b which meets the 1% V ts and E ts error threshold, takes about a minute on a typical desktop computer as compared to ≈15 h taken by the original Attard's model. The results are checked to be independent of the computational domain selection and temporal discretization.
The hysteretic approach-retraction branches in figure 9 demonstrate the co-existence of both repulsive and attractive regimes in a range of A ratio 's. The tip-sample energy dissipation predicted in the repulsive regime follows the same trend previously predicted [35]. However, the method is uniquely capable of calculating the viscoelastic and adhesive dissipation in the attractive regime of AM-AFM operation. This is particularly helpful for sticky and/or highly delicate samples for which imaging in the repulsive regime is challenging and/or sample is susceptible to damage with forces in the repulsive regime.
Finally, we use an algorithm to visualize the tip-surface interaction in a single cycle of AM-AFM on an LDPE sample as characterized in table 1. A series of figures illustrating the tip position and surface geometry at different timesteps are shown in figure 10. The finite-range attractive tip-surface forces cause the surface to deform slightly toward the approaching tip ( figure 10(c)). When retracting, the soft, relaxed surface forms a meniscus around the tip profile and progressively detaches until it completely peels off and continues to gradually return to its initial state (figures 10(g) to (i)). Surface attractive forces cause the contact radius during retraction to become larger than the one during the approach. The surface deformation profiles shown in figure 10 and described above match qualitatively well with in situ confocal microscopy observations of soft polymer surfaces being indented with colloidal probes [36].

Conclusions
This work features an approach to enhance the computational part of Attard's continuum mechanics viscoelasticity contact model. Three enhancements are implemented: (1) the model's formalism is optimized to enable the explicit solution of the governing differential equations, (2) instead of using radial discretization, the deformed surface is reconstructed using a complete set of basis functions, and (3) instead of the forward Euler's method, higher-order numerical procedures are employed to solve the model's ordinary differential equations. By implementing the enhancements, EAM is more than three orders of magnitude faster than Attard's originally proposed computational model. Furthermore, the enhancements improve the computational stability of EAM to better tolerate unstable and abrupt movements of the surface. EAM is a fast first-principles viscoelasticity model that is versatile in terms of the inclusion of various tip-surface interaction forces and surface linear viscoelasticity models. EAM was implemented within the AFM framework to predict force spectroscopy observables and dynamic approach/retract curves by AM-AFM. Moreover, we used EAM to calculate the timeresolved surface dynamics during a single tip-surface interaction cycle of AM-AFM. The excellent agreement between EAM simulation results and the ones by exact approaches verifies EAM's accuracy. The accuracy, speed, and robustness of EAM facilitate simulations for AFM operations on soft, adhesive, and viscoelastic surfaces.