Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 1, 2021, pp. 133 - 153 https://doi.org/10.22190/FUME200920009F © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper GUIDELINES TO SIMULATE LINEAR VISCOELASTIC MATERIALS WITH AN ARBITRARY NUMBER OF CHARACTERISTIC TIMES IN THE CONTEXT OF ATOMIC FORCE MISCROSCOPY Maximilian Forstenhäusler1, Enrique A. López-Guerra1,2, Santiago D. Solares1 1The George Washington University, Department of Mechanical and Aerospace Engineering, Washington DC, USA 2Park Systems Inc., Santa Clara, CA, USA Abstract. We provide guidelines for modeling linear viscoelastic materials containing an arbitrary number of characteristic times, under atomic force microscopy (AFM) characterization. Instructions are provided to set up the governing equations that rule the deformation of the material by the AFM tip. Procedures are described in detail in the spirit of providing a simple handbook, which is accompanied by open-access code and workbook (Excel) sheets. These guidelines seek to complement the existing literature and reach out to a larger audience in the awareness of the interdisciplinary nature of science. Examples are given in the context of force-distance curves characterization within AFM, but they can be easily extrapolated to other types of contact characterization techniques at different length scales. Despite the simplified approach of this document, the algorithms described herein are built upon rigorous classical linear viscoelastic theory. Key Words: Modeling Viscoelasticity, Generalized Maxwell Model, Generalized Kelvin-Voigt Model, Multiple Characteristic Times, Numerical Simulation, Atomic Force Microscopy Received September 20, 2020 / Accepted December 22, 2020 Corresponding author: Santiago D. Solares The George Washington University, Department of Mechanical and Aerospace Engineering, 800 22nd Street NW, Suite 3000 Washington, DC 20052, USA E-mail:ssolares@gwu.edu mailto:ssolares@gwu.edu 134 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES 1. INTRODUCTION Viscoelastic materials are those that can simultaneously store and dissipate energy when deformed [1, 2]. This occurs because they have structural mechanisms to relax some of the stresses that build up when deformed [3, 4]. These relaxations occur at different characteristic times which are often referred to as relaxation times [2, 5]. Viscoelastic materials are very common in nature and diverse scientific fields focus on them. For example, biofilms are known to be viscoelastic and knowledge of their mechanical properties is crucial to understanding how they disseminate and how they could be eradicated [6-9]. Similarly, human cells are viscoelastic and knowledge of their mechanical properties can, for example, help to discriminate cancerous cells from healthy ones in the early stages of the disease [10-12]. In the field of engineering, polymeric fuel cells and organic solar cells have viscoelastic components whose mechanical response in the operation of the devices is known to be ultimately linked to their performance [13-15]. All of these applications can benefit from accurate material modeling approaches to complement the experimental procedures and thus ensure reliable characterization. Unfortunately, due to mathematical complexity, the viscoelastic nature of such materials is often overlooked or oversimplified in scientific studies. Recent investigations [16-26] have explored more rigorous approaches to take into consideration the intricacies of real viscoelastic behaviors by combining the classical theory of viscoelasticity [1, 2, 27-29] with modern characterization techniques, such as atomic force microscopy (AFM) [10, 30-39]. These efforts are especially important to close the gap between the thorough and rigorous mechanical modeling approach of early theories and the interests of scientists that deal with modern techniques. Despite these efforts, we believe that there is still a long way ahead to make this viscoelastic modeling truly accessible to a broader scientific community, regardless of their mathematical background. Therefore, in this study we focus on a more ‘hands-on approach’ where detailed straightforward descriptions to implement viscoelastic modeling are provided. The practical aspects of the modeling implementations are prioritized over mathematical rigor, such that only a very basic knowledge of ordinary differential equations is required to follow our presentation. The reader will be walked through the algorithms of how to solve numerically the constitutive equations that govern the deformation of viscoelastic materials, for the case of an indenter probe that deforms a viscoelastic material with multiple characteristic times. This is of special interest for AFM and for nano- and micro-indentation techniques, although with few adaptations the routines can also be used for other types of characterization techniques. For the convenience of the reader, the algorithm descriptions are accompanied by an Excel spreadsheet available online that contains examples with the numerical calculation methods [40]. Further, users with some background in programming will benefit from the numerical implementations given in Python code in an open-access repository [40]. The repository provided alongside this paper contains two libraries, afmsim and viscoelasticity, as well was implementations and examples of the functions available in these libraries. The contact-mode implementation of AFM will be discussed in Section 3. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...135 2. GUIDELINES FOR THE SIMULATION OF LINEAR VISCOELASTIC BEHAVIOR 2.1. Model Selection: Generalized Maxwell Model & Generalized Kelvin-Voigt Model In this practical guide, generalized models are used to simulate the response of real viscoelastic materials. This provides the necessary flexibility as real viscoelastic materials have multiple characteristic times at which they accommodate and relax stresses when deformed. Specifically, this document will be referring to the Generalized Maxwell Model and the Generalized Kelvin-Voigt Model which are shown in Fig. 1. These two models are mechanical analogs [2] which means that they have the same mechanical response when the appropriate equivalent parameters are chosen. The selection of one over the other obeys pure algebraic convenience, depending on whether force/stress or deformation/strain is regarded as the input. For example, if the user wishes to perform a simulation where the deformation/strain history is given or calculated ‘a priori,’ then the Generalized Maxwell Model is the most convenient. On the other hand, if the force/stress history is known or given, then the Generalized Voigt Model is the most convenient. Fig. 1 a) Generalized Maxwell Model and b) Generalized Kelvin-Voigt Model; Mechanical model diagrams representing the linear viscoelastic relationship between stress and strain in the complex plane with multiple (N) characteristic times when strain or stress is regarded as the excitation, respectively. In a) the Laplace transformed strain, (s), is regarded as the excitation and the Laplace transformed stress, (s), as the response; in (b) the opposite occurs. Jn and Gn refer to the compliance and modulus of the n th spring, respectively. Jg and Ge refer to the glassy compliance and rubbery modulus, respectively. n and n refer to the fluidity and viscosity of the n th dashpot (damper), respectively. Regardless of the chosen model, one must be aware that the model is only a visual representation of the material behavior described by ordinary differential equations relating the stress and strain. These are the governing equations that describe the linear viscoelastic behavior, and their most general form is as follows [20]: ∑ 𝑢𝑛 𝑑𝑛𝜎(𝑡) 𝑑𝑡 𝑛 ∞ 𝑛=0 = ∑ 𝑞𝑚 𝑑𝑚𝜀(𝑡) 𝑑𝑡 𝑚 ∞ 𝑛=0 (1) 136 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES where un and qm are differential coefficients, and (t)and (t)are stress and strain, respectively. This equation is model independent, where the nth and mth time derivatives are acting on the stress and strain tensors, respectively. For mathematical convenience, the differential equation can be transformed into an algebraic equation by applying the Laplace transform. Doing so will transform Eq. (1) to, �̅�(𝑠)𝜎(𝑠) = �̅�(𝑠)𝜀(̅𝑠) (2) where �̅�(𝑠) = ∑ 𝑢𝑛𝑠 𝑛 𝑛 , �̅�(𝑠) = ∑ 𝑞𝑚 𝑠 𝑚 𝑚 and 𝜎(𝑠), 𝜀(̅𝑠) are the transformed stress and strain, respectively. Here, it is assumed zero initial conditions which should not affect generality as they can be incorporated when necessary [2]. Note that �̅�(𝑠) and �̅�(𝑠) are now simply polynomials in the complex variable ‘s’. The interested reader can learn more about the Laplace transform by consulting appropriate references [41-43], although this is not a strict requirement to understand the practical implementation of the algorithms described in this manuscript. We can rearrange Eq. (2) in two ways, depending on whether stress or strain is regarded as the input. For example, when regarding strain as the input, Eq. (2) is rearranged in the following manner: 𝜎(𝑠) = �̅�(𝑠)𝜀(̅𝑠) (3) where �̅�(𝑠) = �̅�(𝑠) 𝑢(𝑠) is the relaxance, which is the mathematical transform carrying the viscoelastic information that will rule the stress response to a given strain input. On the other hand, if the stress is regarded as the input, Eq. (2) can be rearranged as: 𝜀(̅𝑠) = 𝑈(𝑠)𝜎(𝑠) (4) where the retardance, 𝑈(𝑠) = 𝑢(𝑠) �̅�(𝑠) , is introduced. Here, the material transform 𝑈(𝑠), containing the mechanical information of a given material, governs the strain response when an excitation stress is imposed. As can be observed in Eq. (3) and Eq. (4), the retardance is the inverse of the relaxance: �̅�(𝑠) = 1/𝑈(𝑠) (5) In the case of the generalized models in Fig. 1, the material transforms are ratios of polynomials in the complex variable ‘s’. These polynomial expressions, �̅�(𝑠) and �̅�(𝑠), can be expressed in terms of the spring and dashpot (damper) values in the generalized models in Fig. 1. The process to find them is analogous to formulating mesh equations in electric circuit theory [2, 41]. An example for this calculation is given later in Section 2.3. For now, the derivation is skipped and the relaxance and retardance for the generalized Maxwell and Kelvin-Voigt Models are given here: �̅�(𝑠) = 𝐺𝑒 + ∑ 𝐺𝑛𝜏𝑛𝑠 1+𝜏𝑛𝑠 𝑁 𝑛=0 (6) 𝑈(𝑠) = 𝐽𝑔 + ∑ 𝐽𝑛 1+𝜏𝑛 𝑠 𝑁 𝑛=0 (7) where Ge is the equilibrium (rubbery) modulus, which describes the elastic response of the material at long timescales, and Jg is the classy compliance, which describes the elastic response of the material at short timescales [2]. Both quantities are visually Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...137 defined in Fig. 1. N corresponds to the total number of characteristic times in the material. n refers to the n th characteristic time. For the generalized Kelvin-Voigt Model, 𝜏𝑛 = 𝐽𝑛/𝜙𝑛 is the retardation time related to the nth Voigt unit, given by the ratio of the compliance of the n-th spring over the fluidity of the nth dashpot. For the generalized Maxwell Model, n = N / GN is the ratio of the viscosity of the n th dashpot over the modulus of the nth spring. Parameters of this type can be found in the literature for certain materials [17, 44, 45]. The next section describes how to set up the constitutive equation to be solved in our viscoelastic simulation. Here, it will become evident the importance of the mathematical concepts that were just presented in this section. 2.2. Defining the Governing Equation For the specific case of a tip indenting a surface, a customized equation resembling Eq. (1), including geometrical aspects of the physical problem under consideration, is needed. This is because there needs to be a relationship between force and indentation (instead of a relationship between stress and strain as in Eq. (1)), since the AFM instrument measures forces and distances, not stresses. To construct this relationship, geometrical aspects involving a boundary value problem need to be considered for which the correspondence principle may be conveniently invoked. It states that if the elastic solution for a contact problem is known, the viscoelastic solution can be directly formulated as they are analogous in the Laplace domain [46]. The mathematical details are beyond the scope of this manuscript, but curious readers are directed to the literature that discusses its validity [27, 29, 46, 47]. For the case of a spherical indenter with a radius of curvature R, penetrating a viscoelastic half-space (e.g., an AFM tip penetrating a flat viscoelastic surface), the relationship between force Fts(t) and indentation h(t) is [27]: ∑ 𝑞𝑚 𝑑𝑚[{ℎ(𝑡)}3/2] 𝑑𝑡 𝑚 ∞ 𝑚=0 = 3 16√𝑅 ∑ 𝑢𝑛 𝑑𝑛𝐹𝑡𝑠(𝑡) 𝑑𝑡 𝑛 ∞ 𝑛=0 (8) Notice that this equation is very similar to Eq. (1), but the stress has been substituted with the force, and the strain has been substituted with the indentation. Furthermore, to account for the geometry of the system, there is a coefficient on the right-hand side that depends on the radius of curvature of the indenter, and the derivatives on the left-hand side of the equation are taken on the indentation raised to the power 3/2. Eq. (8) assumes that the material is incompressible with a time-independent Poisson’s ratio  = 0.5) [27, 44]. With a finite number of N characteristic times, this equation can also be expressed as: 𝛼 [𝑞0𝑝 + 𝑞1 𝑑𝑝 dt + 𝑞2 𝑑2𝑝 𝑑𝑡 2 +. . . +𝑞𝑁 𝑑𝑁𝑝 𝑑𝑡 𝑁 ] = 𝑢0𝐹+𝑢1 d𝐹 dt +𝑢2 𝑑2𝐹 𝑑𝑡 2 +. . . +𝑢𝑁 𝑑𝑁𝐹 𝑑𝑡 𝑁 (9) where  = 16√𝑅 3 , and p = h(t)3/2. Eq. (9) is our target governing equation that describes the force-indentation relationship of a parabolic tip penetrating a viscoelastic surface. The differential coefficients (u0, u1…, uN, q0, q1…, qN) will depend on the numerical values of the spring moduli and dashpot coefficients in the chosen model (Fig. 1). The next two sections will cover the process of finding these coefficients. Once this is done, the governing equation is completely set up to be solved numerically. Fig. 2 lays out the general process of viscoelastic modeling that we have discussed so far, as well as the remaining steps detailed in next sections. 138 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES As a final note, the correspondence principle is formally invalidated for the retract portion [28, 47], i.e., when the AFM tip is retracting from the surface but still in contact with it. We observed that for the timescales associated with AFM experiments disregarding this fact generally introduces only very small inaccuracies, so we will omit it within this manuscript. To include it, more complex contact-mechanics schemes need to be considered [28, 48], as it has been done on some recent AFM studies [16, 21]. Fig. 2 Summary of the main steps involved in the modeling of a physical problem involving the deformation of a viscoelastic material. These steps are detailed through Sections 2.1 to 2.5 of this manuscript. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...139 2.3. Obtaining the material transform for the viscoelastic model chosen In the last section we targeted the governing equation for the contact mechanics problem of a parabolic tip penetrating a viscoelastic sample, which can be solved with Eq. (9). For simplicity in our illustration of the method, we will now assume that the viscoelastic model chosen has only one characteristic time: The Standard Linear Solid (SLS) model. Also, it will be assumed that deformation can be regarded as input, such that the Maxwell-SLS Model will be the most convenient model to use. Specifically, the Maxwell-SLS Model is a three-parameter model consisting of one spring arm with the equilibrium modulus Ge in parallel with one Maxwell arm (an arm with a spring and a dashpot in series). The Maxwell-SLS Model can be visualized by truncating Fig. 1a to contain only one Maxwell arm (that is, setting N = 1). After selecting an adequate model (e.g., Maxwell-SLS) and selecting the appropriate constitutive equation (e.g., Eq. (9)) describing the force/deformation relationship, the next step is to obtain the material transform for the selected viscoelastic model. The concept of material transform was introduced in Section 2.1. As discussed there, obtaining the material transform expression for a given model involves an analog process to formulating mesh equations in electric circuit theory. Details of this process can be found in the literature [2, 41]. Here, for simplicity, general rules are given as a recipe, without discussing the mathematical or physical justification. These rules for obtaining the material transform are as follows: 1) Relaxances are added up in parallel. For example, for the Maxwell-SLS Model, to obtain the overall relaxance of the model we need to add up the individual relaxances of the individual elements. In this case, we need to add up the relaxance of the isolated spring (Ge) with the relaxance of the Maxwell arm that contains spring G1 and dashpot1. However, before we can perform this addition, we must first determine the relaxance of the Maxwell arm itself, for which the 2nd rule below will be needed. 2) Retardances are added up in series. For example, to obtain the relaxance of the Maxwell arm (a spring in series with a dashpot), we need to add up the retardance of the spring (1/G1) with the retardance of the dashpot (1/1s). Recall from Section 2.1 that relaxances and retardances hold an inverse relationship. Thus, for the Maxwell arm the retardance is 1/1s + 1/G1 and its relaxance is the inverse of this expression, equal to 1 1/η1s+1/G1 . Thus, according to the above two rules, for the example of the Maxwell-SLS Model, the total relaxance is the addition of the isolated spring’s relaxance (Ge) with the Maxwell arm’s relaxance ( 1 1/η1s+1/G1 , which yields 𝐺𝑒 + 1 1/η1s+1/G1 . After some algebraic arrangements and with the substitution n = 1/G1this expression can be transformed into: 𝑄𝑆𝐿𝑆̅̅ ̅̅ ̅̅ (𝑠) = 𝐺𝑒 + 𝐺1𝜏1𝑠 1+𝜏1𝑠 (10) where the characteristic time 1is the ratio of dashpot coefficient to the modulus of the spring in the Maxwell arm. For the generalization of a Maxwell Model with an arbitrary number of characteristic times (Fig. 1a), the relaxances of the rest of the Maxwell arms should be added. It can be easily inferred that following this process would lead us to the generalized expression in Eq. (6), which governs the response of the generalized model in Fig. 1a. 140 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES 2.4 Calculating the Coefficients for the Governing Equation Once we have calculated the material transform (either �̅�(𝑠) or �̅�(𝑠), depending on the chosen viscoelastic model), the next step is to find the differential coefficients (u0, u1…, uN, q0, q1…, qN) from this material transform. To do that, we need to perform algebraic manipulation of the material transforms �̅�(𝑠) and 𝑈(𝑠) such that they can be expressed as a ratio of polynomials in the complex variable ‘s’. For illustrative purposes we continue with the Maxwell-SLS Model example and assign the following arbitrary values for the parameters: Ge = 1.0x10 6 Pa, G1 = 1.0x10 8 Pa, 1 = 1.0x10 -2 s. The first step is to perform the algebraic addition in Eq. (10): 𝑄𝑆𝐿𝑆̅̅ ̅̅ ̅̅ (𝑠) = 𝐺𝑒(1+𝜏1𝑠)+𝐺1𝜏1𝑠 1+𝜏1𝑠 (11) Then, we expand all terms in the numerator and denominator (the denominator will have more than one factor when the model has more than one characteristic time) and gather all the terms with common ‘s’ exponent: 𝑄𝑆𝐿𝑆̅̅ ̅̅ ̅̅ (𝑠) = 𝐺𝑒+(𝐺𝑒+𝐺1)𝜏1𝑠 1+𝜏1𝑠 (12) All the terms that are multiplied by s0 in the numerator are combined and assigned to q0. Similarly, all the terms that are multiplied by s 0 in the denominator are combined and assigned to u0. Then, all terms multiplied by s 1 in the numerator are combined and assigned to q1, while all terms multiplied by s 1 in the denominator are combined and assigned to u1, and so on. Then, it follows that all terms multiplied by s n are grouped and assigned to qn in the numerator and to un in the denominator. Thus, for our specific example we have for the numerator: 𝑞0 = 𝐺𝑒 = 1.0 × 10 6 𝑃𝑎 (13) 𝑞1 = (𝐺𝑒 + 𝐺1)𝜏1 = 1.01 × 10 6 𝑃𝑎 𝑠 (14) And for the denominator we have: 𝑢0 = 1 (15) 𝑢1 = 𝜏1 = 1.0 × 10 −2s (16) Once we have calculated all these differential coefficients, we can plug them into Eq. (9), which concludes the process of explicitly setting up the governing equation. The rest of the manuscript will now focus on explaining the steps to solve this equation numerically (i.e., through computational simulation). As a note of caution, it is worth mentioning that the algebra used to find the differential coefficients becomes quite involved when a large number of characteristic times is chosen. For convenience, the reader is invited to explore the Jupyter Notebook “Calculation of the coefficients, (u0, u1…, uN, q0, q1…, qN), for the differential equation for a 3 Arm Gen. Maxwell, Kelvin-Voigt Model” in an open-access GitHub repository [40]. This Notebook conveniently performs the algebra by means of the SymPy library [49]. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...141 2.5. Numerical Solution of the Governing Equation In the next two subsections, we describe the algorithms to solve the governing equation (Eq. (9)) numerically. The solution is divided into two sections: the first one relates to the portion where tip and sample are in contact during the probe approach and the initial portion of the tip’s retract trajectory, before losing tip-sample contact. The second portion corresponds to the recovery of the surface after tip-sample contact is lost. In this second portion the surface remains temporarily depressed after the tip retracts, followed by surface recovery in a stress-free fashion (the rebound portion [50]). 2.5.1. Numerical Solution for the Indentation into a Viscoelastic Material: Contact Portion of the Force Spectroscopy Curve For the numerical procedure of Eq. (9), the left-hand-side (LHS) and right-hand-side (RHS) of Eq. (9) are separated according to the variable that one wishes to solve for (either force or indentation). In the case of AFM simulations, the known input parameter for Eq. (9) at a given time step is the surface indentation. This is because solution of the AFM tip’s equation of motion (e.g., the harmonic oscillator model [30]) yields the trajectory of the tip, from which we can calculate the indentation (details of this will be provided in Section 3. Thus, the input for the viscoelastic governing equation (Eq. (9)) is the indentation, and the output is the force (the tip-sample force). The first step in solving numerically Eq. (9) is to calculate all higher order derivatives on the LHS using the indentation as the initial input (calculated from the AFM tip trajectory). Specifically, the tip position at the ith time step, hi, provides the information needed to calculate the zero-order derivative on the LHS of the equation, namely the term p = h(t)3/2. However, since we know that depression of the surface will lead a positive (upward) force exerted by the sample on the AFM tip, we must change the sign of hi when calculating the zero-order derivative of the deformation term: 𝑝𝑖 = (−ℎ𝑖 ) 3/2 (17) We can now calculate the rest of the derivatives on the LHS of Eq. (9) using the finite difference method [42, 43]. For example, for the first order derivative ( 𝑑𝑝 dt = �̇�𝑖 ) the numerical approximation using backward difference [42] becomes: �̇�𝑖 = 𝑝𝑖 − 𝑝𝑖−1 Δ𝑡 (18) where pi=(-hi) 3/2 corresponds to the ith time step, as indicated in Eq. (17), and pi-1=(-hi-1) 3/2 is the analogous term corresponding to the previous time step, (i–1). The finite time step is denoted by Δ𝑡. Also note that at t = 0, all higher order derivatives on the LHS are initially set to zero: �̇�𝑖=0 = 0, �̈�𝑖=0 = 0, 𝑝𝑖=0 = 0, etc., and only p0=(-h0) 3/2 has a non-zero value, with h0 being the initial sample deformation at time zero. The higher order derivatives at each time step are determined similarly, using Eq. (19) and Eq. (20): �̈�𝑖 = �̇�𝑖− �̇�𝑖−1 Δ𝑡 (19) … 𝑝𝑖 𝑁 = 𝑝𝑖 𝑁−1 − 𝑝𝑖−1 𝑁−1 Δ𝑡 (20) 142 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES where N is the number of characteristic times in the viscoelastic model. The key for this numerical approach is to respect the order of calculation of derivatives. When the input parameter is the deformation, as in this example, the lowest-order derivative of the deformation term must be calculated first, followed by the calculation of the higher order derivatives in ascending order. Once all the derivatives of indentation have been calculated, it is time to focus on the RHS of Eq. (9). When calculating the RHS of Eq. (9), we use the reverse strategy described in the previous subsection. In this case, we solve first for the highest-order force derivative, and then from it we calculate the lower-order derivatives of the force. As a first step, we rearrange Eq. (9) and solve for the highest force derivative at the ith time step ( 𝑑𝑁𝐹𝑖 𝑑𝑡 𝑁 ) introducing two intermediate variables: 𝑑𝑁𝐹𝑖 𝑑𝑡 𝑁 = 𝑠𝑢𝑚𝑎𝑄𝑖 – 𝑠𝑢𝑚𝑎𝑈𝑖 (21) where 𝑠𝑢𝑚𝑎𝑄𝑖 = 𝛼 𝑢𝑁 [𝑞0𝑝𝑖 + 𝑞1 𝑑𝑝𝑖 dt + 𝑞2 𝑑2𝑝𝑖 𝑑𝑡 2 +. . . +𝑞𝑁 𝑑𝑁𝑝𝑖 𝑑𝑡 𝑁 ] (22) is the LHS of Eq. (9), divided by the coefficient (differential constant) of the highest- order force derivative, and where 𝑠𝑢𝑚𝑎𝑄𝑖 = 𝛼 𝑢𝑁 [𝑞0𝑝𝑖 + 𝑞1 𝑑𝑝𝑖 dt + 𝑞2 𝑑2𝑝𝑖 𝑑𝑡 2 +. . . +𝑞𝑁 𝑑𝑁𝑝𝑖 𝑑𝑡 𝑁 ] (23) is the RHS of Eq. (9), with the highest-order derivative term removed, divided by the differential constant of the highest-order force derivative. Note that in calculating the variable 𝑠𝑢𝑚𝑎𝑈𝑖 we use the force derivatives of the previous time step, (𝑖 − 1), , instead of the current time step. Recall also that  = 16√𝑅 3 for an incompressible material, and p(t) = h(t)3/2 (ensuring that h has a positive sign, as previously discussed – see Eq. (17). Once we have calculated the highest-order derivative, the lower-order force derivatives can be determined by using Euler integration method [42, 43]. 𝐹𝑖 𝑁−1 = 𝐹𝑖−1 𝑁−1 + 𝐹𝑖 𝑁 Δ𝑡 (24) … and we can continue the process until arriving at the lowest force derivative, namely the force itself (Eq. (26)): �̇�𝑖 = �̇�𝑖−1 + �̈�𝑖Δ𝑡 (25) 𝐹𝑖 = 𝐹𝑖−1 + �̇�𝑖Δ𝑡 (26) Eq. (26) is the resulting viscoelastic force response due to indentation, which is transmitted to the AFM tip. This procedure is repeated for each time step dt until the end of the AFM simulation is reached, and in this manner, we obtain the force history. Note that if there has been no previous indentation history (e.g., when the AFM tip first approaches a pristine, undeformed sample during the simulation), the initial conditions (t = 0) for the force and for all its derivatives up to the (N–1)th derivative are equal to zero, 𝐹𝑖=0 𝑁−1 = . . . = �̇�𝑖=0 = 𝐹𝑖=0 = 0. Therefore, at time zero only the highest-order derivative can have a non-zero value. An example calculation for a 3-Maxwell-arm Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...143 model is provided in the next subsection. Fig. 3 summarizes in a flow chart the procedure to calculate the viscoelastic surface deformation that has been discussed so far in this section. Fig. 3 Flow chart describing the calculation steps to obtain the force response for the contact portion of an indenter penetrating a viscoelastic surface. For each timestep of the algorithm, the indentation is regarded as the input for Eq. (9). The calculation steps are detailed in Section 2.5.1. This algorithm is conveniently implemented with Python code in an open-access repository [40] under the afmsim library with the contact_mode function. It is also available in the excel spreadsheet within the same repository. 144 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES 2.5.2. Example – 3 Arm Generalized Maxwell Model To illustrate the process described above, consider a Generalized Maxwell Model with three characteristic times (3 Maxwell arms, refer to Fig. 1a). Again, the procedure begins with the calculation of the deformation derivatives, since the deformation is known: 𝑝𝑖 = (−ℎ𝑖 ) 3/2 (27) �̇�𝑖 = 𝑝𝑖 − 𝑝𝑖−1 Δ𝑡 (28) �̈�𝑖 = �̇�𝑖− �̇�𝑖−1 Δ𝑡 (29) 𝑝𝑖 = �̈�𝑖− �̈�𝑖−1 Δ𝑡 (30) Once all deformation derivatives are determined, the force response can be calculated. 𝑠𝑢𝑚𝑎𝑄𝑖 = 𝛼 𝑢3 [𝑞0𝑝𝑖 + 𝑞1�̇�𝑖 + 𝑞2�̈�𝑖 + 𝑞3𝑝𝑖 ] (31) 𝑠𝑢𝑚𝑎𝑈𝑖 = 1 𝑢3 [𝑢0𝐹𝑖−1 + 𝑢1�̇�𝑖−1 + 𝑢2�̈�𝑖−1] (32) Eq. (21) for the 3-Arm Generalized Maxwell Model will therefore equal: 𝐹𝑖 (𝑡) = 𝑠𝑢𝑚𝑎𝑄𝑖 − 𝑠𝑢𝑚𝑎𝑈𝑖 (33) Using Euler integration, we find the lower-order derivatives and the force itself (zero- order derivative, Eq. (36)): �̈�𝑖 = �̈�𝑖−1 + 𝐹𝑖Δ𝑡 (34) �̇�𝑖 = �̇�𝑖−1 + �̈�𝑖Δ𝑡 (35) 𝐹𝑖 = 𝐹𝑖−1 + �̇�𝑖Δ𝑡 (36) where the initial force and force derivatives at time zero are all equal to zero, �̈�(0) = �̇�(0) = 𝐹(0) = 0. The corresponding implementation in Python code is illustrated in the contact_mode function. This algorithm is also implemented in the Excel spreadsheet provide in the online repository [40]. 2.5.3. Numerical Solution for the Viscoelastic Recovery after the AFM Tip Loses Contact with the Surface during the Retract Portion of the Spectroscopy Curve: The Rebound Problem This portion of the simulation applies only for the retraction portion of the indenter trajectory, after the indenter loses contact with the viscoelastic sample. At this point the temporarily indented viscoelastic surface becomes ‘stress-free’ (since there is no longer an indenter exerting forces on it) and recovers in time according to a deformation profile that depends on its viscoelastic properties and on the previous indentation history. This process is also known as the rebound indentation portion [50]. This portion is seemingly irrelevant for the characterization of viscoelastic materials using AFM force-distance curve methods because there is no way to observe from experimental observables the surface recovery. However, in the simulations it is relevant as we need to track surface Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...145 position for each time step even if the tip temporarily loses contact with the sample. This allows us to have continuity of tip-sample distance knowledge, which is relevant in the calculation of long-range noncontact forces (Section 3.2). Focusing now on the technical details, the onset of this portion in the simulation will be indicated by a change of sign in the Force term. Since forces cannot become negative (assuming the absence of van der Waals forces), as this would indicate that the AFM tip is grabbing and pulling the surface upwards, a flag variable in the simulation should be established to determine when the force changes from positive to negative. From this point on, the stress-free condition begins, and the calculations described in the previous section (2.5.1) are no longer applicable. Specifically, the input parameter for the rebound portion will not be the indentation as in the previous section. Instead, the stress-free condition, which governs the recovery, dictates that in our calculations all force values and force derivatives become zero. Thus, Eq. (9) reduces simply to: 𝛼 [𝑞0𝑝 + 𝑞1 𝑑𝑝 dt + 𝑞2 𝑑2𝑝 𝑑𝑡 2 +. . . +𝑞𝑁 𝑑𝑁𝑝 𝑑𝑡 𝑁 ] = 0 (37) and we solve for the highest-order derivative of the deformation (recall that we are now solving for the deformation profile): 𝑑𝑁𝑝 𝑑𝑡 𝑁 = − 1 𝑞𝑁 (𝑞0𝑝 + 𝑞1 𝑑𝑝 dt + 𝑞2 𝑑2𝑝 𝑑𝑡 2 +. . . +𝑞𝑁−1 𝑑𝑁−1𝑝 𝑑𝑡 𝑁−1 ) (38) The lower-order deformation derivatives can now be determined by using Euler integration: 𝑝𝑖 𝑁−1 = 𝑝𝑖−1 𝑁−1 + 𝑝𝑖 𝑁 Δ𝑡 (39) … �̇�𝑖 = �̇�𝑖−1 + �̈�𝑖 Δ𝑡 (40) 𝑝𝑖 = 𝑝𝑖−1 + �̇�𝑖 Δ𝑡 (41) And finally, using Eq. (27), the sample surface position can be calculated: ℎ𝑖 = 𝑝𝑖 2/3 (42) It is noteworthy that in this procedure the indenter shape parameter (α) has not played a role. This is in accordance to the literature, where it has been noted that for the rebound viscoelastic problem the viscoelastic surface recovery is independent of the indenter shape [50]. Further, as all the viscoelastic forces are zero, due to the stress-free condition, the only forces still present between the tip and the sample are the long-range noncontact forces between the tip and the sample surface, which correspond to dispersion forces (London or van der Waals forces) and are discussed in Section 3.2. The flow chart in Fig. 4 illustrates the procedure to calculate the surface deformation profile as a function of time for the stress-free condition. 146 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES Fig. 4 Flow chart describing the calculation steps to obtain the deformation response for the viscoelastic recovery portion (rebound [50]) when the indenter loses contact with the viscoelastic surface during retract of the indenter. For each timestep the stress-free condition dictates that the force and all its derivatives in Eq. (9) are zero. This allows the calculation of higher order derivatives of indentation in RHS of Eq. (9) to finally obtain the resulting time-dependent surface recovery profile. The calculation steps are detailed in Section 2.5.3. The algorithm has been implemented in Python code and is available in an open-access repository [40] under the afmsim library and the contact_mode function. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...147 3. CASE STUDY OF VISCOELASTIC INDENTATION WITH AN ATOMIC FORCE MICROSCOPE IN ACQUISITION OF FORCE DISTANCE CURVES 3.1. Viscoelastic Implementation in Force-Distance Curve Methods During an off-resonance force-distance curve acquisition, the base of the AFM cantilever is brought towards the surface at a constant velocity. The governing equation that captures the dynamics of the cantilever tip interacting with the surface via the tip-sample forces is [30]: 𝑚�̈�(𝑡) + 𝑐�̇�(𝑡) + 𝑘𝑧(𝑡) = 𝑘𝑧𝑏 (𝑡) + 𝐹𝑡𝑠(𝑡) (43) where t is time, k is the stiffness, c = 2𝜋𝑓𝑚 𝑄1 is the damping coefficient, f the fundamental eigenfrequency (natural frequency), m = 𝑘 (2𝜋𝑓)2 is the equivalent mass, Q1the quality factor of the fundamental eigenmode, �̈�(𝑡)is the tip acceleration, �̇�(𝑡) is the tip velocity, �̇�(𝑡) is the tip position, zb(t) is the cantilever base position and Fts(t) is the tip-sample force. As indicated, the dynamical variables are time-dependent. The initial position of the cantilever base, zb,initial serves as the starting point for the approach process towards the sample surface. During this experiment, the user specifies the cantilever-base approach velocity, �̇�𝑏. 𝑧𝑏 (𝑡) = �̇�𝑏 𝑡 (44) Numerically, the new base position, zb, for each time step, i, is obtained by multiplying the approach velocity by the total simulation time, and adding the result to the initial cantilever base position (the total simulation time is simply the index i multiplied by the time step dt): 𝑧𝑏,𝑖 = 𝑧𝑏,𝑖𝑛𝑖𝑡𝑖𝑎𝑙 + �̇�𝑏,𝑖 (𝑑𝑡) (45) We now calculate the tip acceleration, �̈�(𝑡), by solving the equation of motion (EOM) of the cantilever tip (Eq. (46)) in introducing a subindex corresponding to the time step: �̈�𝑖 = −𝑘𝑧𝑖−1 − 2𝜋𝑓𝑚�̇�𝑖−1 𝑄1 + 𝑘𝑧𝑏,𝑖+𝐹𝑡𝑠,𝑖−1 𝑚 (46) At time zero, Fts = 0, as there is no tip-sample contact yet (the experiment begins with the cantilever placed far away from the sample). z and �̇� are given initial values by the user, z0 = zb,initial and �̇�𝑖 = 0, since we assume that the cantilever was initially at rest and in equilibrium. Using the tip acceleration we calculate the tip position z (Eq. (47)) and velocity �̇� (Eq. (48)) using the Verlet Integration process (a numerical method used to integrate Newton’s Equations of motion) and the Central Difference, respectively [42, 51, 52]. 𝑧𝑖+1 = 2𝑧𝑖 − 𝑧𝑖−1 + �̈�𝑖Δ𝑡 2 (47) �̇�𝑖 = 𝑧𝑖+1− 𝑧𝑖−1 2Δ𝑡 (48) Note that in order to calculate the position for time step (i + 1) with Eq. (47) and the velocity for time step i with Eq. (48), it is necessary to know the position at the previous two time steps i and (i - 1). Therefore, at the beginning of the simulation, it is necessary to define one additional initial condition of position for time step i = -1. One generally 148 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES sets this initial condition to be equal to the initial cantilever base position, namely zi=-1 = zb,initial. In our code we rename the tip position as TipPos, which is given the value zi+1. If TipPosis greater than the current surface position, there is no tip-sample contact. As the cantilever approaches the sample, this variable will become negative, indicating that tip- sample contact has been established and that viscoelastic forces can now be calculated using Eq. (9) as explained in Section 2.5. Fig. 5 provides examples of tip-sample interaction force curves calculated using the above procedures for different sets of parameters, using a 3-Arm Generalized Kelvin- Voigt Model. Note here that the spring parameters for this model are given in terms of compliances (J) instead of moduli (G). As explained in Section 2.1, retardances and relaxances hold an inverse relationship (Eq. (5)). Specifically, the relaxance of a spring is its modulus (a measure of its stiffness) while its retardance measures how soft the spring is: its compliance. Thus, the larger the value of compliance, the softer the spring is. Generally, viscoelastic materials have low values of glassy compliance (Jg), which means that they behave in a stiff-elastic manner when quickly deformed (i.e., at short timescales). On the other hand, viscoelastic materials have larger equilibrium (rubbery) compliances (Je) when probed with very slow excitations (Je = J1 + J2 + … + JN) [1, 2]. Fig. 5 Tip-sample interaction force curves calculated using the contact_mode function contained in the afmsim library [40] for a 3-Arm Generalized Kelvin-Voigt Model. The solid line corresponds to the tip approach and the dash-dotted line corresponds to the tip retract. Hysteresis is evident in all plots, as expected for a (dissipative) viscoelastic material. Fig. 5a provides results for different characteristic time 1, 2x10-2 s, 1.5x10-1 s and 4.5x10-1s. The remaining viscoelastic model parameters are Jg = 2.0x10 -10 Pa-1, J1 =9.0x10 -9 Pa-1, J2=7.0x10 -9 Pa-1, J3 = 1.0x10 -10 Pa-1, 2 = 0.5x10-3 s, 3 = 0.5x10 -2 s. Fig.5b displays results for different cantilever stiffness. The simulation parameters are Jg = 2.0x10 -10 Pa-1, J1 = 5.0x10 -9 Pa-1, J2 = 7.0x10 -9 Pa-1, J3 = 1.0x10-10 Pa-1,1 = 1.5x10 -1 s, 2 = 0.5x10 -3 s, 3 = 0.5x10 -2 s and k = 0.5, 1 and 2 N/m. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...149 3.2 Long-range Noncontact Forces To accurately simulate an AFM experiment, noncontact forces (van der Waals or London dispersion forces) must be accounted for. As the simulation begins, if the tip- position is greater than zero, meaning there is no tip-sample contact yet, the tip and sample will experience an attractive force, which is generally modeled using the Hamaker equation [30]: 𝐹𝑣𝑑𝑊 = − 𝐻∗𝑅 6 (𝑇𝑖𝑝𝑃𝑜𝑠)2 (49) where H is the Hamaker constant, R is the radius of curvature of the AFM tip (assumed to be nearly spherical) and TipPos gives the tip-sample distance (the distance between the tip and the unperturbed surface). The negative sign indicates that the forces are attractive (i.e., the tip experiences a downward force). When there is no tip-sample contact, Eq.(49) describes the only force present, so Fts the tip-sample force in the cantilever tip equation of motion, Eq. (43)) must be calculated using Eq. (49), and there is no viscoelastic force component. Once tip-sample contact is established, the tip-sample force, Fts, has two contributions: the first one is the viscoelastic force described above, Fi, and the second one is the attractive tip-sample force from Eq. (49). When the tip is pushing onto the surface, the tip-sample distance is set to≈0.2 nm, which corresponds approximately to the diameter of a single atom of average size (that is, we assume that the atoms of the tip and the sample are touching one another, so the centers of the atoms in contact are one atomic diameter away). The total tip-sample force is then: 𝐹𝑡𝑠 = 𝐹𝑖 − 𝐻∗𝑅 6 𝑎2 (50) For simplicity, we assume that the van der Waals forces are unable to deform the surface. Strictly speaking this is not true, since attractive tip-sample forces are in some cases able to pull the surface upwards and cause it deform above the initial unperturbed position. This can be important when modeling the imaging of materials having very low stiffness. 3.3 Time Step The time step, or iteration step, is of central importance, as it has a strong effect on the accuracy of the results and significantly affects the stability of numerical methods. Numerical methods are prone to instabilities when the time step is too large and can diverge very rapidly. On the other hand, extremely small time steps bring about high memory demands and lead to unnecessarily large computation time. Therefore, the right balance between too small and too large a time step needs to be established. In selecting the time step, one should consider the total simulation time required, the smallest characteristic time in the material and the period of the AFM cantilever (the inverse of its natural frequency). The smaller the total simulation time, the smaller the time step must be in order to ensure stability (for the case of high deformation rates). Instead of recommending a specific time step that provides stability, a guideline for ranges of time steps seems more applicable. It is recommended that a time step smaller than the smallest characteristic time divided by 10 or the fundamental period T of the cantilever divided by 10,000 be used for the simulation to obtain stable results. The lower limit depends on the computational resources available. Since every simulation has its own peculiarities, these general suggestions may not always be valid, so the optimum time step must in the end be determined by trial and error. 150 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES 3.4 Verification Upon successfully running a numerical simulation, for example using the contact_mode function we have provided [40], the question remains as to whether it is actually accurate and correct. One way to verify the results is by comparing the results of the left-hand-side (LHS) with the right-hand-side (RHS) of the following equation, the Boltzmann Superposition Integral, adjusted for a spherical tip according to Lee and Radok [17]. 16√𝑅 3 ℎ(𝑡)3/2 = ∫ 𝑈(𝑡 − 𝜉)𝐹(𝜉)𝑑𝜉 𝑡 0 (51) where the LHS equals 𝐿𝐻𝑆 = 𝛼 ∗ ℎ(𝑡)3/2 (52) and the RHS, the convolution of force with retardance (U*F), equals 𝑅𝐻𝑆 = 𝐽𝑔 𝐹(𝑡) + ∑ ∫ 𝐽𝑛 𝜏𝑛 𝑒 −(𝑡−𝜉)/𝜏𝑛 𝑑𝜉 𝑡 0𝑛 (53) where 𝑈(𝑡) = 𝐽𝑔 + ∑ 𝐽𝑛 𝜏𝑛 𝑒 −𝑡/𝜏𝑛𝑛 is the compliance of the material. Recall that= 16√𝑅 3 , h(t) and F(t) are the indentation and the force, respectively, as before. 𝜉is the integration variable. Once the simulation data is collected, one can use the deformation data to calculate the LHS, and the force data to calculate the RHS. The RHS can be calculated with the conv function in the viscoelasticity library in our GitHub repository [40]. Calculating the RHS requires passing the force and time arrays from the simulation to the numerical convolution described in Eq. (53). This can also be done using the conv function in the viscoelasticity library in our GitHub repository [40]. When doing this, it is sometimes computationally very expensive to pass the entire force and time arrays, which could be very large if the time step of the simulation is very short. Instead, one could pass a shorter (scattered) version of the arrays containing data every defined number of time steps. If the latter is decided for computational convenience, a word of caution should be given as to how scattered the version of the arrays should be. Specifically, we strongly advice following the time step boundaries provided in Section 3.3. An example of the consequence of passing too scattered force and time arrays is illustrated in Fig. 6. This plot verifies the simulation as correct, if the convolution integral (RHS) and the results from the simulation (tip position multiplied by α constant: LHS) overlay each other. As can be seen, when the arrays are too scattered (the time step is too large), LHS and RHS do not overlay (see the blue line labeled with a time step dtconv = 9.9x10 -5 s). In contrast, the yellow line (with a time step dtconv = 9.9x10 -7) perfectly overlays the scattered purple star points (LHS). Note that for this verification we have not considered van der Waals tip-sample forces. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...151 Fig. 6 Impact of computing the convolution of force with retardance (RHS, Eq. (53)) with arrays that have different time step values. The convolutions were computed with the following time step values dtconv = 9.9x10 -5 s; 9.9x10-6 s and 9.9x10-7 s and plotted with different colors as indicated in the figure’s legend. These results are compared with Eq. (53) (LHS), which is plotted with purple stars. It can be observed that a time step of approximately 9.9x10-7 s has to be passed to the convolution integral (RHS, Eq. (53)) in order to match the results from the simulations (LHS, Eq. (52)). Note how the yellow line matches the purple star markers. 𝑈(𝑡) = 𝐽𝑔 + ∑ 𝐽𝑛 𝜏𝑛 𝑒 −𝑡/𝜏𝑛𝑛 is the retardance of the material, h is the indentation, F(t) the force, the integration variable, and the coefficient is = 16√𝑅 3 . 4. CONCLUSION We have outlined detailed guidelines for implementing linear viscoelastic simulations in the context of AFM force-distance curve methods, where an indenter deforms a linear viscoelastic material. The manuscript has been written in a practical manner in an effort to make it accessible to a broad audience of researchers and students of different disciplines. In this spirit, procedural details have been prioritized over mathematical details, for which the user has been referred to relevant literature. Specifically, we have provided detailed explanations for setting up the governing equation ruling the behavior of the linear viscoelastic material being indented, and the subsequent numerical solution of this equation. For the numerical solutions, we have provided conceptual explanations and given specific examples, along with ready-to-use Python codes and Excel spreadsheets. Despite the simplicity of the presentation we have not sacrificed rigor, as the modeling has been placed in the context of complex viscoelastic materials possessing multiple characteristic times. This manuscript targets a large audience desiring to perform analysis and characterization of viscoelastic materials in a rigorous manner, but having only a moderate background in the field of linear viscoelasticity. 152 M. FORSTENHAEUSLER, E. A. LÓPEZ-GUERRA, S. D. SOLARES REFERENCES 1. Ferry, J.D., 1980, Viscoelastic Properties of Polymers, John Wiley & Sons. 2. Tschoegl, N.W., 2012, The Phenomenological Theory of Linear Viscoelastic Behaviour: an Introduction, Springer Science & Business Media. 3. Rouse Jr, P.E., 1953, A theory of the linear viscoelastic properties of dilute solutions of coiling polymers, The Journal of Chemical Physics, 21(7), pp. 1272-1280. 4. Plazek, D., 1993, Breakdown of the Rouse model for polymers near the glass transition temperature, The Journal of chemical physics, 98(8), pp. 6488-6491. 5. Roylance, D., 2001, Engineering viscoelasticity, Department of Materials Science and Engineering– Massachusetts Institute of Technology, Cambridge MA, 2139, pp. 1-37. 6. Gordon, V.D., 2017, Biofilms and mechanics: a review of experimental techniques and findings, Journal of Physics D: Applied Physics, 50(22), 223002. 7. Kovach, K., 2017, Evolutionary adaptations of biofilms infecting cystic fibrosis lungs promote mechanical toughness by adjusting polysaccharide production, npj Biofilms and Microbiomes, 3(1), 1. 8. López-Guerra, E.A., Shen H., Solares, S.D., Shuai, D., 2019, Acquisition of time–frequency localized mechanical properties of biofilms and single cells with high spatial resolution, Nanoscale, 11(18), pp. 8918- 8929. 9. Shen, H., López-Guerra, E.A., Zhu, R., Diba, T., Zheng, Q., Solares, S.D., Zara, J.M., Shuai, D., Shen, Y., 2018, Visible-light-responsive photocatalyst of graphitic carbon nitride for pathogenic biofilm control, ACS applied materials & interfaces, 11(1), pp. 373-384. 10. Krisenko, M.O., Cartagena, A., Raman, A., Geahlen, R.L., 2015, Nanomechanical property maps of breast cancer cells as determined by multiharmonic atomic force microscopy reveal Syk-dependent changes in microtubule stability mediated by MAP1B, Biochemistry, 54(1), pp. 60-68. 11. Lekka, M., Laidler, P., 2016, Applicability of AFM in cancer detection, Nature nanotechnology, 4(2), pp. 72-72. 12. Plodinec, M., Loparic, M., Monnier, C.A., Obermann, E.C., Zanetti-Dallenbach, R., Oertle, P., Hyotyla, J.T., Aebi, U., Bentires-Alj, M., Lim, R.Y.,Schoenenberger, C.A., 2012, The nanomechanical signature of breast cancer, Nature nanotechnology,7(11), pp. 757-765. 13. Arrechea, S., Aljarilla, A., de la Cruz, P., Palomares, E., Sharma, G.D, Langa, F., 2016, Efficiency improvement using bis (trifluoromethane) sulfonamide lithium salt as a chemical additive in porphyrin based organic solar cells, Nanoscale, 8(41), pp.17953-17962. 14. Bruner, C., Dauskardt, R., 2014, Role of molecular weight on the mechanical device properties of organic polymer solar cells, Macromolecules, 47(3), pp. 1117-1121. 15. Noh, H., Diaz, A.J., Solares, S.D., 2017, Analysis and modification of defective surface aggregates on PCDTBT: PCBM solar cell blends using combined Kelvin probe, conductive and bimodal atomic force microscopy, Beilstein Journal of Nanotechnology, 8, pp. 579-589. 16. Efremov, Y.M., Wang, W.H., Hardy, S.D., Geahlen, R.L., Raman, A., 2017, Measuring nanoscale viscoelastic parameters of cells directly from AFM force-displacement curves, Scientific reports, 7(1), pp. 1-14. 17. Zhai, M., McKenna, G.B., Viscoelastic modeling of nanoindentation experiments: A multicurve method, Journal of Polymer Science Part B: Polymer Physics, 52(9), pp. 633-639. 18. López‐Guerra, E.A., Eslami, B., Solares,S.D., 2017,Calculation of standard viscoelastic responses with multiple retardation times through analysis of static force spectroscopy AFM data, Journal of Polymer Science Part B: Polymer Physics, 55(10), p. 804-813. 19. López‐Guerra, E.A., Solares, S.D., 2017, Material property analytical relations for the case of an AFM probe tapping a viscoelastic surface containing multiple characteristic times, Beilstein Journal of Nanotechnology, 8(1), pp. 2230-2244. 20. Garcia, P.D., Guerrero, C.R., Garcia, R., 2017, Time-resolved nanomechanics of a single cell under the depolymerization of the cytoskeleton, Nanoscale, 9(33), pp. 12051-12059. 21. Garcia, P.D., Guerrero, C.R., Garcia, R., 2020, Nanorheology of living cells measured by AFM-based force– distance curves, Nanoscale, 12(16), pp. 9133-9143. 22. Parvini, C.H., Saadi, M.A.S.R., Solares, S.D., Extracting viscoelastic material parameters using an atomic force microscope and static force spectroscopy, Beilstein Journal of Nanotechnology, 11(1), pp. 922-937. 23. Rajabifar, B., Jadhav, J.M., Kiracofe, D., Meyers, G.F. Raman, A., 2018, Dynamic AFM on viscoelastic polymer samples with surface forces, Macromolecules, 51(23), pp. 9649-9661. 24. Chyasnavichyus, M., Young, S.L., Tsukruk, V.V., 2015, Recent advances in micromechanical characterization of polymer, biomaterial, and cell surfaces with atomic force microscopy, Japanese Journal of Applied Physics, 54(8S2), 08LA02. 25. Radmacher, M.,Tillmann, R., Gaub, H., 1993, Imaging viscoelasticity by force modulation with the atomic force microscope, Biophysical journal, 64(3), pp. 735-742. Guidelines to Simulate Linear Viscoelastic Materials with an Arbitrary Number of Characteristic Times...153 26. Garcia, R., 2020, Nanomechanical mapping of soft materials with the atomic force microscope: methods, theory and applications, Chemical Society Reviews, 49(16), pp. 5850-5884. 27. Lee, E., Radok, J.R.M., 1960, The contact problem for viscoelastic bodies, Journal of Applied Mechanics, 27(3), pp. 438-444. 28. Ting, T., 1966, The contact stresses between a rigid indenter and a viscoelastic half-space, Journal of Applied Mechanics, 33(4), pp. 845-854. 29. Graham, G.A., 1965, The contact problem in the linear theory of viscoelasticity, International Journal of Engineering Science, 3(1), pp. 27-46. 30. Garcia, R., Perez, R, 2002, Dynamic atomic force microscopy methods, Surface science reports, 47(6), pp. 197-301. 31. Melcher, J., Hu, S., Raman, A., 2008, Invited Article: VEDA: A web-based virtual environment for dynamic atomic force microscopy, Review of Scientific Instruments, 79(6), 061301. 32. Guzman, H.V., Garcia, P.D., Garcia,R., 2015, Dynamic force microscopy simulator (dForce): A tool for planning and understanding tapping and bimodal AFM experiments, Beilstein Journal of Nanotechnology, 6(1), pp. 369-379. 33. Attard, P., 2007, Measurement and interpretation of elastic and viscoelastic properties with the atomic force microscope, Journal of Physics: Condensed Matter, 19(47), pp. 473201. 34. Amo, C.A., Garcia, R., 2016, Fundamental high-speed limits in single-molecule, single-cell, and nanoscale force spectroscopies, ACS nano, 10(7), pp. 7117-7124. 35. Cartagena, A., Raman, A., 2014, Local viscoelastic properties of live cells investigated using dynamic and quasi-static atomic force microscopy methods, Biophys J, 106(5), pp. 1033-43. 36. Garcia, P.D., Garcia, R., 2018, Determination of the elastic moduli of a single cell cultured on a rigid support by force microscopy, Biophysical Journal, 114(12), pp. 2923-2932. 37. Garcia, R., 2006, Identification of nanoscale dissipation processes by dynamic atomic force microscopy, Physical review letters, 97(1), 016103. 38. Herruzo, E.T., Perrino, A.P., Garcia, R., 2014, Fast nanomechanical spectroscopy of soft matter, Nature Communications, 5(1), pp. 1-8. 39. Hu, S., Raman, A., 2008, Inverting amplitude and phase to reconstruct tip–sample interaction forces in tapping mode atomic force microscopy, Nanotechnology, 19(37), 375704. 40. Forstenhäusler, M., López‐Guerra, E.A., 2020, J.L. AFMviscoelastic Github Repository, Available from: https://github.com/mforstenhaeusler/AFMviscoelastic. 41. Gardner, M.F., Barnes, J.L., 1956, Transients in Linear Systems Studied by the Laplace Transformation, J. Wiley & Sons. 42. Kreyszig, E., 2007, Advanced Engineering Mathematics, John Wiley & Sons. 43. Stroud, K.A., Booth, D.J., 2011, Advanced Engineering Mathematics, Palgrave Macmillan. 44. Brinson, H.F., Brinson, L.C., 2008, Polymer Engineering Science and Viscoelasticity, Springer. 45. Simon, S.L., Mckenna, G.B, Sindt, O., 2000, Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation, Journal of Applied Polymer Science, 76(4), pp. 495-508. 46. Lee, E., 1995, Stress analysis in visco-elastic bodies, Quarterly of Applied Mathematics, pp. 183-190. 47. Graham, G., 1968, The correspondence principle of linear viscoelasticity theory for mixed boundary value problems involving time-dependent boundary regions, Quarterly of Applied Mathematics, 26(2), pp. 167-174. 48. Popov, V.L., Willert, E., Heß, M., 2018, Method of dimensionality reduction in contact mechanics and friction: A user’s handbook. III, Viscoelastic contacts, Facta Universitatis-Series Mechanical Engineering, 16(2), pp. 99- 113. 49. Meurer, A., Smith, C.P., Paprocki, M., Čertík, O., Kirpichev, S.B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J.K., Singh, S. Rathnayake, T., 2017, SymPy: symbolic computing in Python, PeerJ Computer Science, 3, 103. 50. Argatov, I.I., Popov, V.L., 2016, Rebound indentation problem for a viscoelastic half‐space and axisymmetric indenter—Solution by the method of dimensionality reduction, ZAMM‐Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 96(8), pp. 956-967. 51. Khan, I.R., Ohba, R., 2003, Taylor series based finite difference approximations of higher-degree derivatives, Journal of Computational and Applied Mathematics, 154(1), pp. 115-124. 52. Grubmüller, H., 1991, Generalized Verlet algorithm for efficient molecular dynamics simulations with long- range interactions, Molecular Simulation, 6(1-3), pp. 121-142.