Abstract
Atomic force microscopy (AFM) techniques have provided and continue to provide increasingly important insights into surface morphology, mechanics, and other critical material characteristics at the nanoscale. One attractive implementation involves extracting meaningful material properties, which demands physically accurate models specifically designed for AFM experimentation and simulation. The AFM community has pursued the precise quantification and extraction of ratedependent material properties, in particular, for a significant period of time, attempting to describe the standard viscoelastic response of materials. AFM static force spectroscopy (SFS) is one approach commonly used in pursuit of this goal. It is capable of acquiring rich temporal insight into the behavior of a sample. During AFMSFS experiments the cantilever base approaches samples with a nearly constant velocity, which is manipulated to investigate different timescales of the mechanical response. This manuscript seeks to build upon our previous work and presents an approach to extracting useful linear viscoelastic information from AFMSFS experiments. In addition, the basis for selecting and restricting the model parameters for fitting is discussed from the perspective of applying this technique on a practical level. This work begins with a guided discussion that develops a fit function from fundamental laws, continues with conditioning a raw SFS experimental dataset, and concludes with the fit and prediction of viscoelastic response parameters such as storage modulus, loss modulus, loss angle, and compliance. These steps constitute a complete guide to leveraging AFMSFS data to estimate key material parameters, with a series of detailed insights into both the methodology and supporting analytical choices.
Introduction
Modern AFM applications commonly involve testing samples that are soft, biological, or polymeric in nature. Understanding the dissipative nature of these materials at the nanoscale is especially important to their use in many applications [16]. However, performing such measurements using AFM has been difficult due to the complexity of material phenomena that govern AFM observables. For example, the history and ratedependence of tip–sample interactions can produce large errors when estimating the surface stiffness, leading to poor and inconsistent data quality. To avoid complex analytical derivations, it can be convenient to ignore wellestablished viscoelastic models in favor of elastic relationships. Unfortunately, these equations are sometimes oversimplified to such an extent that they no longer properly represent the behavior of real materials [7,8]. Such approaches critically lack the capability of reproducing the complex temporal behavior of samples with one or more characteristic times. It follows that more detailed viscoelastic models are necessary to describe and predict AFM tip–sample behavior accurately [912]. To this end, an appropriate viscoelastic model must be chosen (or derived) and subsequently specified for a geometry specific to AFM tip–sample interactions.
The process of developing and fitting a viscoelastic model is somewhat organic. One approach is to use theoretical arguments in conjunction with continuum mechanics fundamentals, but it is most common to design more intuitive spring–dashpot linear viscoelastic mechanical models [1315]. The spring–dashpot, or “mechanical equivalent” approach to linear viscoelasticity involves connecting different combinations of springs and dashpots in series and/or parallel to mimic the action of a material. These physical systems are then described algebraically, transformed into Laplace space, and rearranged to create transfer functions that describe the material response for a given stress or strain excitation. The mechanicalequivalent approach is simple to explain, but can require more assumptions and some additional knowledge of the Laplace transform to derive analytical stress–strain relationships. Alternatively, continuum mechanics can be used to create more general models to describe the material under study. The continuum method involves the derivation and restriction of constitutive relations that apply to practical problems in mechanics. These relations are then subjected to balance laws in order to develop a pipeline from problem statement to final configuration [16]. The continuum approach is more general at the cost of requiring significant expertise and complex explanations that are more difficult to follow.
Either method is capable of generating a model that provides significantly more insight into the material response than elastic models, and the choice between approaches usually depends upon the material and complexity of interactions being represented. In the case of AFM, many common problems involve nanoscale viscoelastic systems exposed either to intermittent probe contact (tapping mode) or a similar type of probe–sample interaction where the indentations are usually small and fast. In such cases mechanicalequivalent models usually suffice, and the continuum approach is more complicated than necessary to reproduce the material response. When the contact occurs for longer periods, such as during contactmode AFM or when tip–sample nonlinearities are clearly visible in the dataset, continuum models may be necessary to capture all of the deformation complexity. Using either technique, the goal remains to create a physically accurate model of the sample, and use experimental data to quantify meaningful material properties.
Lopez et al. [17] recently formulated a generalized solution to the physical contact problem in AFMSFS experiments, particularly for cases where the sample is viscoelastic. The approach is also capable of being adapted to any linear viscoelastic stress–strain model. Their method utilized a theoretical basis for viscoelastic indentation problems proposed by Lee and Radok [18] in 1960. Lee and Radok presented a framework that used Hertzian contact relationships and the “viscoelastic correspondence principle”, in which elastic parameters are replaced by their viscoelastic counterparts to account for differences in the mechanical response of the material. Approaches utilizing viscoelastic correspondence were then cited frequently in the literature of that time [1921].
The aforementioned work of Lopez et al. focused upon materials with an arbitrary number of retardation times and included analytical descriptions both with and without the assumption of a linear forcevstime curve. This generalized approach neatly ties together the theory of Lee and Radok with a modern understanding of linear viscoelastic mechanicalequivalent (i.e., spring–dashpot) models, and the present article seeks to provide a more detailed discussion on practical applications of the methodology of Lopez and coworkers. Topics discussed here include the basis in theory, how to specify various material models for fitting, how to calculate relevant viscoelastic quantities, and importantly a discussion of where the extracted parameters are valid. The extension of this methodology to nonlinear viscoelastic materials is the subject of a future manuscript in development. In order to begin extracting viscoelastic parameters, the solution outlined by Lee and Radok [18] must be revisited and subsequently specified for an AFMSFS experiment as described by Lopez and coworkers [17]. Afterwards, the relationships are generalized for an arbitrary load history and assumptions are made about the fundamental characteristics of the material being studied. The “Theoretical Background” section focuses on the analytical approach, with the detailed data handling outlined in the section immediately afterward. For clarity, a specific example is then provided in the “Results and Discussion” section, showing how to extract viscoelastic parameters using AFMSFS data collected from hollow nylon 6,6 tubing.
Theoretical Background
To understand the analytical choices available for viscoelastic models in AFM indentation problems, it is helpful to begin with the viscoelastic correspondence principle framework mentioned above. An AFM tip–sample interaction involves measuring force and deformation, which can then be compared to the solution of a contact problem where a spherical probe tip indents a flat viscoelastic surface. Depending upon the material being tested, different models can be inserted into the resulting force–indentation relationship and then compared to AFMSFS data. The parameter values within the model are modified until a close approximation is made between the experimental force and indentation datasets, resulting in a set of “bestfit” viscoelastic parameters for the model. The following sections explain this workflow in a similar manner: First, the framework is described briefly, then the force–indentation relationship is modified to be more easily compared with AFM data, and last, the features of several viscoelastic models are discussed. For a detailed derivation of the Lee and Radok framework, and a review of the extension to arbitrary load history (previously presented by Lopez et al. [17]), the reader is directed to Supporting Information File 3 included with this manuscript.
Spherical indentation of a viscoelastic halfspace
In an ideal case, the exact tip geometry is known a priori. In lieu of such information, it is expedient to make simplifying assumptions about the tip shape at the point of contact. One common assumption is that the local surface deformations are small, such that the tip has a roughly spherical contact geometry during indentation. This assumption is best made for largerdiameter tips or specialized colloidal probes, although it can be argued that smaller tips can still be approximated provided the indentation depth is small. As mentioned above, Lee and Radok [18] proposed a solution to the rigid spherical contact of a viscoelastic halfspace nearly sixty years ago. The indentation configuration is visualized in Figure 1, which is based upon Figure 1 in their original manuscript.
Here, the deepest indentation occurs at the center of the sphere (r = 0) and is labeled h(t). The indenter has a radius of curvature R, and the distance from the center axis to the edge of contact, known as the contact radius, is labeled l(t).
Using the geometry presented in Figure 1, a series of algebraic manipulations, Laplace transformations, and viscoelasticequivalence arguments, Lee and Radok presented the following relationship between applied load and spherical indentation depth for a viscoelastic material having characteristic operators u and q
Equation 2 is the relationship upon which Lopez et al. build a solution for the AFM experiment. It is originally presented as Equation 3 in their paper [17].
Extending the solution to an arbitrary load history
Traditionally, when using creeprecovery experiments to parameterize the viscoelastic models under study, a constant stress is first applied to a sample and then subsequently removed. For example, one method involves hanging weights from a material for some time, then removing the weights and allowing the material to recover. A strain gage or displacement sensor captures the deformation occurring during both phases, and the data is used for fitting. The boundary conditions for such an approach dictate a step function in stress, and that one end of the sample is fixed. In that case, the constant load history would suggest using the “creep compliance” J(t), an engineering quantity that represents the change in strain as a function of the time for a medium exposed to an instantaneous constant stress [17,18]. While applicable for this type of study, the tip–sample interaction during AFM experiments does not apply a constant force (i.e., stress) as a function of the time. The load history is more reminiscent of a discrete impulse function in which the contact time is short. Therefore, while the solution form could use the creep compliance for fitting, it is more direct to use the material retardance. The creep compliance can be defined in terms of the applied stress σ(t), the strain ε(t), and the material retardance as:
By taking the Laplace transform of Equation 2, rearranging, substituting the creep compliance relationship above, and taking the inverse Laplace transform, one arrives at the following result:
Equation 5 allows for the straightforward definition of the retardance U(t) according to an appropriate material model. The convolution integral replaces the multiplication of and in the Laplace domain during translation back to the time domain. As mentioned previously, the approach outlined here is preferred for AFM as it both removes the requirement of a step function in stress and of a measurement of the force application rate. Since an AFMSFS experiment provides force and indentation (deformation), the data streams can be utilized directly without requiring a discrete derivative of force in time. Taking a discrete derivative can introduce undesirable oscillatory errors into the resulting dataset, and thus obscure some of the information contained within. Equation 5 represents the strain (or deformation) response of a material to the unit stress impulse.
Specifying a viscoelastic material model for fitting
Selecting an appropriate material model is often an iterative process, unless one can leverage available datasets or a priori knowledge about the nature of a given sample. There is a broad array of literature that covers spring–dashpot models and viscoelastic rheology, many of which provide helpful visualizations of the predicted material response [14,15,22]. Provided a relationship between model parameters and the material retardance can be derived, one can replace U(t − ζ) in the above equations to make an implicit assumption of how the material reacts to applied stress. Acquiring these relationships is generally done by taking the derivative of the creep compliance J(t) of the model, as the two quantities are related in Laplace space by Equation 4. Table 1 summarizes several commonly used viscoelastic models and their reported applicability.
Table 1: Summary of viscoelastic models and the corresponding applications [14,15,22].
type  name  shape  stress–strain equation  parameters 
linear  Kelvin–Voigt^{a}  2  
standard linear solid (Voigt form)^{b}  3  
generalized Kelvin–Voigt^{c}  1 + 2n  
Maxwell^{d}  2  
Wiechert^{e}  1 + 2n  
nonlinear  Schapery^{f}  —  7+ 
^{a}Model includes exponential (and reversible) strain creep under stress, but no stress relaxation with constant strain. The response will decay to zero residual stress provided sufficient time.
^{b}Model includes both creep and relaxation phenomena, but is not accurate when literally compared to common material responses.
^{c}If including the fluidity term, an additional dashpot is added in series with the mechanical model. This model accounts for creep strain occurring at multiple retardation times.
^{d}Model includes exponential (and reversible) stress relaxation under constant strain, but linear (nonreversible) strain creep under constant stress.
^{e}Most easily applied to relaxation experiments (constant strain). This model accounts for relaxation occuring at multiple time scales.
^{f}By including additional terms in ϕ and ϕ^{′}, more environmental effects can be accounted for and the number of parameters increases.
Often these linear spring–dashpot descriptions utilize an expanded series of compliances and characteristic time constants, which introduces a common source of error during analysis of the viscoelastic data. Each term in the series is intended to represent the reaction of the material to applied stress on a particular time scale, represented by τ in the governing equations shown in Table 1. As stated previously, provided the model has a derived retardance, it can be inserted into Equation 5 directly without an inherent restriction on the number of characteristic terms included in the series. Often additional terms are inserted with the intent of more accurately fitting the dataset in question. However, if the data does not include information at those timescales, whether due to a lack of resolution (for small time scales) or limitations on the experiment length (long time scales) these additional terms will introduce significant error in the parameters extracted. Furthermore if the parameter set obtained is extrapolated for new loading conditions, it will incorrectly predict the material response as a result of an incorrect number of terms (i.e., branches in the generalized mechanicalequivalent model).
It is critical to evaluate the timescales contained within a dataset before performing a parameter extraction, especially, before predicting new material responses. For example, since modern AFM instruments are capable of acquiring data on the scale of 2 kHz for quasistatic characterization, the smallest characteristic timescales that can be resolved are of the order of 10^{−4} s. If parameters are fit to a dataset of the order of 10^{−4} to 10^{0} s (e.g., data collected at 2 kHz for a total experiment length of 1–3 s), it is improper to utilize that parameter set to mimic the action of a material for an experiment lasting more than 10 s. Issues of this kind can lead to significant problems with the performance of a parameterized model. As technology continues to advance and AFM setups can sustain faster data acquisition speeds, the resulting datasets will continue to become richer in timedependent material information and provide a more complete understanding of how the sample material responds to stress.
To maintain consistency between the approach outlined here and that of Lopez et al., the generalized Kelvin–Voigt model has been selected for analysis. The retardance is usually found by taking the derivative of the compliance of a model [17]:
Several new parameters have been introduced. The first is the “glassy compliance” (J_{g}), representing the elastic response of the material at short timescales; the second is the characteristic time (τ_{n}), which scales the third parameter, the characteristic compliance (J_{n}), such that it contracts on a specific timescale corresponding to the nth branch of the generalized Voigt model. By inserting Equation 7 into Equation 5, we obtain:
Note that the final term containing the steadystate flow ϕ_{f} has been removed as a higherorder effect with relatively small magnitude for mostly solid viscoelastic materials. Materials that exhibit steadystate flow under stress, such as plastic or glass at high temperatures, are thus excluded from the analysis shown here. If this term was included under special circumstances, the set of fit parameters would increase by a number determined by the steadystate flow model utilized for ϕ_{f} and an additional numerical convolution would become necessary.
Equation 9 allows for extracting parameters for n characteristic terms without any assumptions about the load history, provided the penetration grows monotonically, the indentation is not “sufficiently large”, and the correspondence principle applies. It represents the objective fit function described by Lopez et al. [17] and requires a numerical convolution of U(t) and F(t) in terms of the material parameter set {J_{g}, J_{n}, τ_{n}}.
Useful viscoelastic quantities
When characterizing the response of viscoelastic materials to external stress, especially during cyclic loading (such as with dynamic mechanical analysis (DMA) machines), it is common to calculate harmonic quantities such as the storage modulus E′, loss modulus E′′, and loss angle δ [2325]. Each has a physical interpretation, representing the elastic and viscous motion of the material and their magnitudes relative to one another. For example, a material that is very stiff will have a high storage modulus and a low loss modulus. Such a sample will tend to store a majority of the applied load within its molecular structure and elastically return most or all of that energy when unloaded. Alternately, a medium that is susceptible to large shear forces (such as fluids) will have a low storage modulus and high loss modulus. In this case, most of that input energy will be lost to friction and heat, and therefore the material will return far less energy than the stiff elastic material when unloaded. To acquire storage and loss modulus as functions of the frequency, a Fourier transform of the compliance equation is applied to a particular model. Once translated into the frequency domain, the real portion of the expression corresponds to the storage modulus, and the imaginary portion to the loss modulus. In the case of a generalizedVoigt material, a definition in terms of modulus is analytically inconvenient. The model is most clearly written using material compliances (J_{n}), and as such the storage compliance and the loss compliance (J′ and J′′, respectively) are used. The corresponding moduli can be written as [26]:
For ideally elastic behavior, the storage modulus or compliance dominates and the loss angle tends towards zero. In contrast, when the dissipative capacity of the material dominates, the sample action is more fluidlike and the loss angle increases. Understanding these viscoelastic quantities in the context of material characterization can be helpful when describing the relationship between external loads and the viscoelastic response, especially when the excitations are periodic. Within the context of AFM, both the storage modulus and the loss modulus are critical to evaluating the dissipated energy during tappingmode analysis [24]. The key equations required for viscoelastic parameter inversion are as follows:
Equation 2, general Lee and Radok solution for spherical indentation [18]:
Equation 5, integral form of Lee and Radok spherical indentation solution [17]:
Equation 9, Lopez solution for spherical indentation of generalizedVoigt materials with multiple retardation times [17]:
Equation 10, storage compliance for generalizedVoigt viscoelastic materials [26]:
Equation 11, loss compliance for generalizedVoigt viscoelastic materials [26]:
Equation 12, loss angle for generalizedVoigt viscoelastic materials [26]:
Results and Discussion
Extracting viscoelastic parameters
The following sections discuss how to leverage AFMSFS data for viscoelastic model parameterization, including how to condition raw datasets for fitting, extract material model parameters, and calculate the viscoelastic quantities mentioned above. While the current approach is primarily geared towards MATLAB implementation, the original process was outlined by Lopez et al. [17] in Python, and is available in a public Github repository.
Conditioning raw static force spectroscopy datasets
Traditionally, AFMSFS experiments generate a variety of data streams for postprocessing, including ZRaw (expected Zposition based on the voltage applied to the Zpiezo), ZSensor (measured Zposition), deflection, amplitude, and phase. The ZSensor measurement corresponds to the position of the piezo that controls the position of the base of the cantilever, and the timederivative of the ZSensor data should be equal to the experimental approach velocity. To perform the fit procedure, only the deflection and ZSensor datasets are necessary. Figure 2a illustrates the shape of a common result as read from a force–distance curve output file.
Because the Lee and Radok derivation used previously is restricted to a monotonically increasing surface area, the approach portion of the dataset must be isolated for analysis. The maximum ZSensor point marks the beginning of retraction, and therefore all points occurring after the time where that reading occurs (t_{max(}_{z}_{)}, marked as t_{2} in Figure 2a) are removed from the dataset. The resulting curve is illustrated in Figure 2b.
Due to a number of possible factors, e.g., calibration and noise, the deflection in the region t_{0} to t_{1} can be nonzero despite the probe not being in contact with the surface. To correct for this error, one needs to first locate the minimum deflection point. From there, the user needs to determine a sufficient step backwards in time such that the probe has yet to experience both longrange attractive forces and “snapin”. The resulting time is marked as t_{1}, and the average deflection offset is calculated using the deflection data from t_{0} to t_{1} (recall that the deflection should be set to zero at t_{0}). The offset is then subtracted to produce the corrected deflection (d_{c}). The dataset should now show the deflection as roughly zero for the period before snapin, as seen in Figure 3a.
It remains to adjust the dataset such that the point of minimum deflection during snapin coincides with the origin of the deflectionvsZSensor curve. This step simplifies calculating the slope of the curve where the force is increasing as a function of the time. While not explicitly necessary, smoothing the deflection data with a simple firstorder Butterworth filter can remove some of the measurement noise and make finding the deflection minimum easier. While searching for the smallest value without filtering, noise can cause neighboring points to appear lower than the actual minimum and introduce error. It is convenient at this stage to store both the ZSensor and the deflection values at the time of minimum deflection as z_{0} and d_{0}. The corrected ZSensor values (z_{c}) are then calculated according to:
The applied tip force (F(t)) and indentation depth (h(t)) vectors can then be calculated using the cantilever stiffness (k_{c}), the minimum deflection offsets (z_{0}, d_{0}), the deflection (d(t)), and the following relations:
where k_{c} in Equation 14 is the AFM cantilever stiffness. This data is not used directly in the fit, but can be useful for showing the indentation and force during the approach of the tip into the sample. The final conditioning step requires extracting the repulsive portion from both of the corrected datasets (z_{c}(t) and d_{c}(t)). While there are a variety of methods for identifying where the repulsive data starts, an easily implemented method involves creating a linear approximation of the repulsive ZSensor data z_{rep}. The slope is taken from the corrected ZSensor data within the region noted in Figure 3b, which begins at the time where the deflection minimum occurs (t_{min(}_{d}_{)}), and the number of points n_{p} in that region is counted. Generating the z_{rep} vector is then performed by creating an array from 0 to n_{p} in steps of 1 and multiplying by . Subtracting the repulsive potion of the deflection (d_{rep}) from z_{rep} provides the tip position relative to its neutral point:
The data utilized for fitting is the region where z_{tip} is greater than zero, indicating the tip has begun penetrating the surface and coinciding with the start of force application on the sample. Trimming data before this point provides the final form of the time array t, the repulsive deflection d_{rep}, and the repulsive ZSensor values z_{rep} utilized in the following sections. In the following, the subscript “rep” will refer to these final forms where the force application has begun. To conclude the data conditioning, Equation 14 and Equation 15 are used again to calculate the experimental force F and the indentation depth h with repulsive data. At this point it is recommended to save resampled, logscale versions of both F and h for use during the fit process. Because Equation 9 will be parameterized using a nonlinear leastsquares approach, the number of data points per decade will drive the fit quality on that timescale. Ensuring that there are an equal number of data points in each log order will lead to consistent weighting across the entire experimental time scale [17].
The conditioning thus far has focused on preparing a single AFMSFS dataset. However, most experiments involve multiple runs at a single site, which is then repeated at multiple locations on the sample to acquire a broader description of the surface properties. The mechanical responses can be averaged to both reduce computational overhead during fitting and give a prediction that is less sensitive to inhomogeneities. However, care must be taken to avoid issues that can arise when each dataset is taken at slightly different times, for varying overall lengths of time, and especially for different sampling rates. To create a valid “averaged” estimate from a large number of datasets, each curve must be used independently to predict the response for a consistent set of input conditions. Implementing this approach first requires the user to find the range of “common” ZSensor values that all of the datasets contain, and generate a list of numbers between those values. This array will act as the “average ZSensor” in the following. Next, the ZSensor and the deflection vectors from each dataset are provided to a 1D interpolation function (such as “interp1()” in MATLAB), which accepts a sample set of xvalues (the ZSensor dataset), a sample set of yvalues (the deflection dataset), and a set of “query points” to evaluate the interpolation with. The average ZSensor array created previously will act as the “input” (i.e., query points) to interp1(), and the output will be the interpolated deflection values as predicted by a given dataset. This process continues for every AFMSFS file used in the average, with each supplying a new set of xvalues, a new set of yvalues, and each being evaluated with the same input ZSensor array. A time array is then generated for each run by using the sampling frequency and number of points in that dataset as an estimate. Then, a spline is performed to resample the estimated time array such that it contains the exact same number of points as contained in the average ZSensor array. At the end of this step, the user can directly average the predicted deflection and time arrays from all datasets to obtain the “average deflection” and “average time”.
Clearly, when creating a list of numbers for the average ZSensor array at the beginning, it is critical to use an appropriate number of steps. Depending upon the number of points used, it is possible to artificially change the sampling rate to appear faster or slower, thus erroneously changing the temporal data contained within the averaged result. For example, if all of the individual experimental datasets contain roughly 1000 data points in the range of valid ZSensor values, and the average ZSensor array is made using an arbitrary 10,000 points, the average time vector will appear to have a timestep that is ten times faster than expected. This simple mistake can lead to the averaged dataset containing information that is not actually present in the force–distance files being processed. The number of points used to generate the average ZSensor array should be calibrated to match the number of points expected for a dataset of that length in the given experiment, such that the effective sampling frequency does not appear to change drastically. The most convenient means of verifying that the averaged dataset is valid involves looking at the average time array generated from the spline process, and ensuring that its timestep is close to the timestep of the individual data files. For example, if using a sampling frequency of 2 kHz for the AFMSFS experiments, the average time array should appear to have a timestep of roughly 5 × 10^{−4} s.
Estimating material model parameters using a nonlinear leastsquares approach
Modeling the force–indentation relationship with Equation 9 can be performed using a variety of techniques. In their paper, Lopez et al. utilize the costfunction minimization approach. The method involves comparing the lefthand side of Equation 9 to the predicted righthand integral, which is calculated using a specific set of model parameters (J_{g}, τ_{n}, J_{n}). These parameters are then varied by the builtin algorithm (specifically Levenberg–Marquardt in the case of Python’s minimize() function) to solve the nonlinear leastsquares minimization problem. The function searches for the “optimal” parameter set that minimizes the relative error between the dataset and integral prediction.
Due to differences in the fit functions available, it is inconvenient to implement the Levenberg–Marquardt algorithm in MATLAB with the necessary upper and lower bounds on the model parameters. Instead, several nested “anonymous” functions (i.e., functions defined on a single line instead of in a separate file) are used to simulate a single experimental dataset. The toplevel function is passed to lsqcurvefit(), which will call all of the subfunctions while fitting the parameters. The MATLAB function lsqcurvefit() is also capable of taking an input matrix as opposed to simply an array of data, meaning it could fit multiple datasets simultaneously if properly implemented. While this capability is not explored here, the ability to fit multiple averaged data sets from different experiments could prove useful in the future.
The series of nested anonymous functions are configured according to the following workflow:
 F_{conv}: Performs the convolution of using the “full” setting, resulting in an array that is N + M − 1 in length, where N is the number of points in U(t), and M is the number of points in F(t);
 subref: Accepts the convolved result of F_{conv} and trims all data but a subset equal to the size of t_{rep};
 F_{conv,wrap}: Scales the output of subref by the time step dt, adds the glassy compliance term, and resamples the result to logarithmic scaling. This is the function provided to lsqcurvefit().
Because lsqcurvefit() calculates the cost internally, it is not necessary to define a function for this purpose. However, other approaches may require an additional “comparison” function. The default algorithm used by lsqcurvefit() is “trustregionreflective”, which is a nonlinear leastsquares method, similar to Levenberg–Marquardt. A comparison of each approach and their efficacy is beyond the scope of this manuscript, but the reader is directed to the literature for explanations of both algorithms [27,28]. It is important to note that the parameter set obtained from these algorithms is not necessarily unique, and repeatedly performing the fit for a statistically significant number of times is recommended before evaluating its performance. Furthermore, the solutions are sensitive to the initial guess and parameter bounds. To counteract these effects the starting point for each parameter can be generated randomly within the upper and lower bounds for a fitting run, which is in turn performed many times to develop an understanding of the parameter space.
Beyond the function definitions and solver algorithms, it remains to appropriately determine and enforce the upper and lower bounds for each parameter, and also to choose starting “guesses” for each leastsquares run. According to the specification of Lopez et al. [17], each successive term in the generalizedVoigt compliance function is intended to represent the material response during one order of magnitude in time. The time constants (τ_{n}) are thus limited to a factor of ten lower and higher, with each additional term accounting for the next largest scale in the series. For example, a “oneterm” series would contain the glassy compliance term and the compliance of a single Voigt element with a characteristic time centered of the order of the data acquisition frequency. Such a series would allow the characteristic time (τ_{1}) to vary from 10^{−5} to 10^{−3} in the case of sampling at 10 kHz, and the compliance (J_{1}) to be restricted to all reasonable values in the range (0,1). These terms should theoretically account for experiments with a maximum effective time of 0.01 to 1 ms, but for cases where the contact duration reaches 1 s or more this term alone would be insufficient to reproduce the material response. To adequately model these datasets, five or more terms would be necessary reaching characteristic times centered on 10^{0} or greater.
As mentioned previously, selecting an appropriate range of characteristic times is a critical step in extracting the most applicable parameter set. Without including a sufficient number of terms, the material response at large scales would dominate the lower orders, thus forcing each characteristic time parameter to tend toward their respective upper bounds. If such action is observed during the fit process, it is recommended to include an additional term and repeat the process until at most one order larger than the experiment length has been added, or the parameters are no longer converging on their upper or lower bounds.
It is evident from Equations 10–12 that each parameter set will also generate unique harmonic quantities. By examining the loss angle (δ(ω)), parameter groups could be evaluated and eliminated based upon their feasibility. In the case where a wellstudied material is being investigated, this task can be relatively straightforward. However, for new materials without significant background literature it remains difficult to discard parameter sets based upon their harmonic quantities. The following section examines the fit procedure applied to nylon, and the effects of experimental settings on the quality of fit provided by the presented methodology.
Extracting viscoelastic parameters from experimental data for nylon
A hollow nylon tube (source: McMasterCarr, P/N 8628K48, Grade 6/6) was cut, ground smooth with a handheld rotary tool, and cleaned in isopropyl alcohol (IPA) followed by a bath in deionized (DI) water. The sample was then allowed to dry sufficiently and mounted on an AFM metal specimen disk using doublesided carbon tape. The force mapping (six scan lines, six points per scan line) technique was implemented using an Asylum MFP3D AFM to accommodate a variety of sample configurations, which yielded 36 indentation datasets at 36 different locations. Furthermore three approach velocities were selected in a logarithmic distribution: 10, 100, and 1000 nm/s. The probe utilized was an OLYMPUS AC 240TSR3, featuring a tip radius of roughly 10 nm. Before measurement, the tip was calibrated using the thermal noise method [29] in which a hard silicon sample was used after sonicating using Mucasol, followed by IPA, and then DI water for a duration of 10 min per step. The resonant frequency and spring constant were found to be 68.953 kHz and 1.70 nN/nm, respectively.
The curves are aligned as a result of the previously described conditioning steps for each velocity, and are visualized in Figure 4 and Figure 5. For clarity, the repulsive portions of each curve have also been isolated to show the exact data utilized for fitting Equation 9.
As expected, the fit qualities vary based upon the approach speed and number of fit terms used. This is due to varying amounts of longtimescale characteristics being present. For shorter experiments, using fewer terms resulted in a closer approximation while keeping characteristic time parameters away from their upper and lower bounds. The compliances generated in those cases are more reasonable, indicated by viscoelastic moduli of the order of several gigapascals, which agree with stiffnesses reported in the literature [30]. The parameter sets that were rated optimal for each approach velocity are given in Table 2.
Table 2: Generalized Kelvin–Voigt viscoelastic model parameter sets resulting from nonlinear leastsquares fits of AFMSFS force–indentation data. Characteristic times are given in units of [s], and compliances are provided in units of [Pa^{−1}]. Note that the parameter sets utilized for generating the final fit are marked with an asterisk (*).
velocity [nm/s]  parameter  number of elements  

1  2  3  4*  5  
10  J_{g}  3.76e11  3.87e12  3.79e11  1.51e11  2.65e14 
J_{1}  3.4e09  2.22e14  7.39e11  3.5e12  4.98e12  
τ_{1}  9.56e04  2.22e05  2.98e04  5.1e05  1.1e04  
J_{2}  3.94e09  2.75e10  1.09e09  4.62e10  
τ_{2}  9.99e03  1.01e03  8.89e03  4.53e03  
J_{3}  5.4e09  2.11e09  1.98e09  
τ_{3}  9.94e02  7.55e02  3.45e02  
J_{4}  8.24e09  7.51e09  
τ_{4}  7.52e01  5.01e01  
J_{5}  2.49e09  
τ_{5}  8.05  
1  2  3*  4  5  
100  J_{g}  4.86e14  2.88e13  2.15e11  6.14e11  1.87e10 
J_{1}  2.8e09  2.22e14  1.06e10  3.2e14  3.16e10  
τ_{1}  9.71e04  2.32e04  7.58e04  5.15e04  9.89e04  
J_{2}  4.63e09  3.27e10  5.1e10  7.05e11  
τ_{2}  9.99e03  6.33e04  1.06e03  8.13e03  
J_{3}  4.94e09  4.61e09  4.68e09  
τ_{3}  1.72e02  1.7e02  1.66e02  
J_{4}  4.66e09  2.23e09  
τ_{4}  7.45e01  9.02e01  
J_{5}  3.02e09  
τ_{5}  3.76  
1  2  3*  4  
1000  J_{g}  5.1e14  3.73e14  3.65e11  2.48e14  
J_{1}  3.18e09  1.24e09  1.64e09  6.92e10  
τ_{1}  9.96e04  5.82e04  7.35e04  9.7e04  
J_{2}  7.64e09  2.22e14  9.77e10  
τ_{2}  9.46e03  6.26e03  6.11e04  
J_{3}  3.51e08  1.71e08  
τ_{3}  6.08e02  3.85e02  
J_{4}  1.25e07  
τ_{4}  9.3e01 
After evaluating the fit quality for each velocity, it was determined that four terms provided the best data approximation for the 10 nm/s dataset, while three terms were sufficient for 100 and 1000 nm/s. Figure 6 shows the fitted viscoelastic model against each corresponding dataset. Note that the parameterization was unsuccessful for the 5Voigt series fit on the 1000 nm/s averaged data because the conditioned repulsive dataset contained less than eleven data points. This is significant because lsqcurvefit() will not operate on models where the number of parameters within (eleven terms for the 5Voigt series) is greater than the number of data points used for comparison.
Harmonic quantities have also been calculated using the fit parameters from Figure 6, providing the predicted relationship between the elastic and viscous material response. The results are provided in Figure 7. It is critical to reiterate that each set of optimization parameters is not unique. There is a theoretically unlimited number of distinct parameter groupings that could recreate the material response for a given experiment. With each parameter set, the calculation of Equations 10–12 will indicate different distributions for the observed harmonic quantities. To illustrate the possible variations, Figure 8 visualizes the predictions generated by the topten ranked parameter sets and their apparent fit to the averaged 10 nm/s approach velocity dataset.
Clearly, estimating the harmonic properties of a material sample can be difficult using a single AFMSFS experiment. As the number of test sites on the sample surface increases, the larger data population should theoretically limit the number of suitable parameter sets that recreate all of the material responses. Similarly, it could also be worthwhile to attempt to fit an even larger number of datasets using a variety of approach velocities simultaneously. Considering multiple approach velocities would reduce the parameter space, eliminate the need for a comparison of the results between approach velocities, and could provide a more generally applicable set of parameters after one script iteration. The implementation of a multiplevelocity fit process would be significantly more complex, and would require dynamically eliminating model terms for different velocities, such that each dataset would be compared to a corresponding model with an appropriate number of characteristic times. To avoid adding this additional degree of freedom to the problem, it would be necessary to have a priori knowledge of the optimal number of terms to use for each velocity. This would remove the need to vary the number of terms dynamically and drastically reduce the associated computational overhead. Regardless of the approach, determining an effective method for extracting the most broadly applicable parameters, and by extension their corresponding harmonic quantities, remains an attractive subject for further research.
Conclusion
A methodology for extracting viscoelastic parameters from a sample material response in AFMSFS experiments has been discussed and applied to the specific case of nylon 6,6 using multiple approach velocities. Beginning with the correspondence principle for the spherical indentation of a viscoelastic halfspace, the presented inversion approach has been shown to provide mechanical and temporal insight into the viscoelastic response of nanoscale surfaces. As stated previously, the assumptions utilized during the objective function derivation places inherent limits on the appropriate experimental settings. Most critically, the indentation strain measures must remain small to prevent violation of the underlying theory. Common pitfalls and limitations have been discussed, with particular emphasis on both data conditioning and mechanicalequivalent model design. One important issue discussed was that the number of characteristic terms included in compliance series models must be intentionally prescribed based upon the length and data acquisition rate of a given experiment. If the effects of experiment length and the desired application time range are not considered, an otherwise appropriate model could be incorrectly parameterized and perform poorly, or be applied outside its valid temporal range. There remain many attractive questions worthy of investigation, including further restriction of the bestfit parameters using their implied harmonic quantities, extending the theoretical implementation to include multiaxial datasets, and the application of this inversion approach to nonlinear compliance models – the last of which is the subject of a future paper by the authors.
Supporting Information
The Supporting Information includes a directory containing all necessary MATLAB scripts used for AFMSFS file inversion, and an instructional document to explain the operation and settings available in the script. Additionally, the full derivations mentioned in section “Theoretical Background” are included in Supporting Information File 3.
Supporting Information File 1: Readme file for MATLAB script operation.  
Format: PDF  Size: 800.9 KB  Download 
Supporting Information File 2: Compressed directory containing the MATLAB script and supporting functions.  
Format: ZIP  Size: 1.8 MB  Download 
Supporting Information File 3: More detailed derivations of the relationships presented in section ‘‘Theoretical Background’’.  
Format: PDF  Size: 319.5 KB  Download 
References

Dittmer, J. J.; Lazzaroni, R.; Leclère, P.; Moretti, P.; Granström, M.; Petritsch, K.; Marseglia, E. A.; Friend, R. H.; Brédas, J. L.; Rost, H.; Holmes, A. B. Sol. Energy Mater. Sol. Cells 2000, 61, 53–61. doi:10.1016/s09270248(99)000963
Return to citation in text: [1] 
Mayol, L.; Quaglia, F.; Borzacchiello, A.; Ambrosio, L.; La Rotonda, M. I. Eur. J. Pharm. Biopharm. 2008, 70, 199–206. doi:10.1016/j.ejpb.2008.04.025
Return to citation in text: [1] 
Lekka, M.; Laidler, P. Nat. Nanotechnol. 2009, 4, 72. doi:10.1038/nnano.2009.004
Return to citation in text: [1] 
Eslami, B.; LópezGuerra, E. A.; Diaz, A. J.; Solares, S. D. Nanotechnology 2015, 26, 165703. doi:10.1088/09574484/26/16/165703
Return to citation in text: [1] 
Lekka, M. Bionanosci. 2016, 6, 65–80. doi:10.1007/s1266801601913
Return to citation in text: [1] 
Guerrero, C. R.; Garcia, P. D.; Garcia, R. ACS Nano 2019, 13, 9629–9637. doi:10.1021/acsnano.9b04808
Return to citation in text: [1] 
Attard, P. J. Phys.: Condens. Matter 2007, 19, 473201. doi:10.1088/09538984/19/47/473201
Return to citation in text: [1] 
López Guerra, E. A. Analytical Developments for the Measurement of Viscoelastic Properties with the Atomic Force Microscope. Ph.D. Thesis, The George Washington University, Washington, DC, USA, 2018.
Return to citation in text: [1] 
Proksch, R.; Kocun, M.; Hurley, D.; Viani, M.; Labuda, A.; Meinhold, W.; Bemis, J. J. Appl. Phys. 2016, 119, 134901. doi:10.1063/1.4944879
Return to citation in text: [1] 
Efremov, Y. M.; Wang, W.H.; Hardy, S. D.; Geahlen, R. L.; Raman, A. Sci. Rep. 2017, 7, 1541. doi:10.1038/s41598017017843
Return to citation in text: [1] 
Garcia, P. D.; Garcia, R. Nanoscale 2018, 10, 19799–19809. doi:10.1039/c8nr05899g
Return to citation in text: [1] 
Benaglia, S.; Amo, C. A.; Garcia, R. Nanoscale 2019, 11, 15289–15297. doi:10.1039/c9nr04396a
Return to citation in text: [1] 
Tschoegl, N. W. The phenomenological theory of linear viscoelastic behavior: an introduction; Springer Science & Business Media, 2012.
Return to citation in text: [1] 
Brinson, H. F.; Brinson, L. C. Polymer Engineering Science and Viscoelasticity; Springer US: Boston, MA, U.S.A., 2008. doi:10.1007/9780387738611
Return to citation in text: [1] [2] [3] 
LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2014, 5, 2149–2163. doi:10.3762/bjnano.5.224
Return to citation in text: [1] [2] [3] 
Chen, W.q. J. Zhejiang Univ., Sci., A 2014, 15, 231–240. doi:10.1631/jzus.a1400079
Return to citation in text: [1] 
LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327
Return to citation in text: [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] 
Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020
Return to citation in text: [1] [2] [3] [4] [5] [6] 
Graham, G. A. C. Q. Appl. Math. 1968, 26, 167–174. doi:10.1090/qam/99860
Return to citation in text: [1] 
Ting, T. C. T. J. Appl. Mech. 1966, 33, 845–854. doi:10.1115/1.3625192
Return to citation in text: [1] 
Hunter, S. C. J. Appl. Mech. 1961, 28, 611–617. doi:10.1115/1.3641792
Return to citation in text: [1] 
Kelly, P. Solid mechanics lecture notes. http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/Part_I/BookSM_Part_I/10_Viscoelasticity/10_Viscoelasticity_Complete.pdf (accessed May 19, 2020).
Return to citation in text: [1] [2] 
Schwarzl, F. R. Rheol. Acta 1971, 10, 165–173. doi:10.1007/bf02040437
Return to citation in text: [1] 
LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2017, 8, 2230–2244. doi:10.3762/bjnano.8.223
Return to citation in text: [1] [2] 
Liu, C.; Matseevich, T.; Askadskii, A.; Popova, M.; Serenko, O.; Zhang, X. AIP Conf. Proc. 2017, 1839, 020065. doi:10.1063/1.4982430
Return to citation in text: [1] 
Iqbal, S. Characterization of Viscoelastic Properties of a Material Used for an Additive Manufacturing Method. Ph.D. Thesis, University of North Texas, TX, USA, 2013.
Return to citation in text: [1] [2] [3] [4] 
Levenberg, K. Q. Appl. Math. 1944, 2, 164–168. doi:10.1090/qam/10666
Return to citation in text: [1] 
Yuan, Y.x. A review of trust region algorithms for optimization. In ICM99: Proceedings of the Fourth International Congress on Industrial and Applied Mathematics, 2000; pp 271–282.
Return to citation in text: [1] 
Hutter, J. L.; Bechhoefer, J. Rev. Sci. Instrum. 1993, 64, 1868–1873. doi:10.1063/1.1143970
Return to citation in text: [1] 
McKeen, L. W. 8  Polyamides (Nylons). In Film Properties of Plastics and Elastomers, 3rd ed.; McKeen, L. W., Ed.; Plastics Design Library; William Andrew Publishing: Boston, MA, USA, 2012; pp 157–188. doi:10.1016/b9781455725519.000086
Return to citation in text: [1]
26.  Iqbal, S. Characterization of Viscoelastic Properties of a Material Used for an Additive Manufacturing Method. Ph.D. Thesis, University of North Texas, TX, USA, 2013. 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
1.  Dittmer, J. J.; Lazzaroni, R.; Leclère, P.; Moretti, P.; Granström, M.; Petritsch, K.; Marseglia, E. A.; Friend, R. H.; Brédas, J. L.; Rost, H.; Holmes, A. B. Sol. Energy Mater. Sol. Cells 2000, 61, 53–61. doi:10.1016/s09270248(99)000963 
2.  Mayol, L.; Quaglia, F.; Borzacchiello, A.; Ambrosio, L.; La Rotonda, M. I. Eur. J. Pharm. Biopharm. 2008, 70, 199–206. doi:10.1016/j.ejpb.2008.04.025 
3.  Lekka, M.; Laidler, P. Nat. Nanotechnol. 2009, 4, 72. doi:10.1038/nnano.2009.004 
4.  Eslami, B.; LópezGuerra, E. A.; Diaz, A. J.; Solares, S. D. Nanotechnology 2015, 26, 165703. doi:10.1088/09574484/26/16/165703 
5.  Lekka, M. Bionanosci. 2016, 6, 65–80. doi:10.1007/s1266801601913 
6.  Guerrero, C. R.; Garcia, P. D.; Garcia, R. ACS Nano 2019, 13, 9629–9637. doi:10.1021/acsnano.9b04808 
16.  Chen, W.q. J. Zhejiang Univ., Sci., A 2014, 15, 231–240. doi:10.1631/jzus.a1400079 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
13.  Tschoegl, N. W. The phenomenological theory of linear viscoelastic behavior: an introduction; Springer Science & Business Media, 2012. 
14.  Brinson, H. F.; Brinson, L. C. Polymer Engineering Science and Viscoelasticity; Springer US: Boston, MA, U.S.A., 2008. doi:10.1007/9780387738611 
15.  LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2014, 5, 2149–2163. doi:10.3762/bjnano.5.224 
14.  Brinson, H. F.; Brinson, L. C. Polymer Engineering Science and Viscoelasticity; Springer US: Boston, MA, U.S.A., 2008. doi:10.1007/9780387738611 
15.  LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2014, 5, 2149–2163. doi:10.3762/bjnano.5.224 
22.  Kelly, P. Solid mechanics lecture notes. http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/Part_I/BookSM_Part_I/10_Viscoelasticity/10_Viscoelasticity_Complete.pdf (accessed May 19, 2020). 
9.  Proksch, R.; Kocun, M.; Hurley, D.; Viani, M.; Labuda, A.; Meinhold, W.; Bemis, J. J. Appl. Phys. 2016, 119, 134901. doi:10.1063/1.4944879 
10.  Efremov, Y. M.; Wang, W.H.; Hardy, S. D.; Geahlen, R. L.; Raman, A. Sci. Rep. 2017, 7, 1541. doi:10.1038/s41598017017843 
11.  Garcia, P. D.; Garcia, R. Nanoscale 2018, 10, 19799–19809. doi:10.1039/c8nr05899g 
12.  Benaglia, S.; Amo, C. A.; Garcia, R. Nanoscale 2019, 11, 15289–15297. doi:10.1039/c9nr04396a 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
7.  Attard, P. J. Phys.: Condens. Matter 2007, 19, 473201. doi:10.1088/09538984/19/47/473201 
8.  López Guerra, E. A. Analytical Developments for the Measurement of Viscoelastic Properties with the Atomic Force Microscope. Ph.D. Thesis, The George Washington University, Washington, DC, USA, 2018. 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
29.  Hutter, J. L.; Bechhoefer, J. Rev. Sci. Instrum. 1993, 64, 1868–1873. doi:10.1063/1.1143970 
19.  Graham, G. A. C. Q. Appl. Math. 1968, 26, 167–174. doi:10.1090/qam/99860 
20.  Ting, T. C. T. J. Appl. Mech. 1966, 33, 845–854. doi:10.1115/1.3625192 
21.  Hunter, S. C. J. Appl. Mech. 1961, 28, 611–617. doi:10.1115/1.3641792 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
30.  McKeen, L. W. 8  Polyamides (Nylons). In Film Properties of Plastics and Elastomers, 3rd ed.; McKeen, L. W., Ed.; Plastics Design Library; William Andrew Publishing: Boston, MA, USA, 2012; pp 157–188. doi:10.1016/b9781455725519.000086 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
27.  Levenberg, K. Q. Appl. Math. 1944, 2, 164–168. doi:10.1090/qam/10666 
28.  Yuan, Y.x. A review of trust region algorithms for optimization. In ICM99: Proceedings of the Fourth International Congress on Industrial and Applied Mathematics, 2000; pp 271–282. 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
14.  Brinson, H. F.; Brinson, L. C. Polymer Engineering Science and Viscoelasticity; Springer US: Boston, MA, U.S.A., 2008. doi:10.1007/9780387738611 
15.  LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2014, 5, 2149–2163. doi:10.3762/bjnano.5.224 
22.  Kelly, P. Solid mechanics lecture notes. http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/Part_I/BookSM_Part_I/10_Viscoelasticity/10_Viscoelasticity_Complete.pdf (accessed May 19, 2020). 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
26.  Iqbal, S. Characterization of Viscoelastic Properties of a Material Used for an Additive Manufacturing Method. Ph.D. Thesis, University of North Texas, TX, USA, 2013. 
26.  Iqbal, S. Characterization of Viscoelastic Properties of a Material Used for an Additive Manufacturing Method. Ph.D. Thesis, University of North Texas, TX, USA, 2013. 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
17.  LópezGuerra, E. A.; Eslami, B.; Solares, S. D. J. Polym. Sci., Part B: Polym. Phys. 2017, 55, 804–813. doi:10.1002/polb.24327 
24.  LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2017, 8, 2230–2244. doi:10.3762/bjnano.8.223 
18.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
23.  Schwarzl, F. R. Rheol. Acta 1971, 10, 165–173. doi:10.1007/bf02040437 
24.  LópezGuerra, E. A.; Solares, S. D. Beilstein J. Nanotechnol. 2017, 8, 2230–2244. doi:10.3762/bjnano.8.223 
25.  Liu, C.; Matseevich, T.; Askadskii, A.; Popova, M.; Serenko, O.; Zhang, X. AIP Conf. Proc. 2017, 1839, 020065. doi:10.1063/1.4982430 
26.  Iqbal, S. Characterization of Viscoelastic Properties of a Material Used for an Additive Manufacturing Method. Ph.D. Thesis, University of North Texas, TX, USA, 2013. 
© 2020 Parvini et al.; licensee BeilsteinInstitut.
This is an Open Access article under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0). Please note that the reuse, redistribution and reproduction in particular requires that the authors and source are credited.
The license is subject to the Beilstein Journal of Nanotechnology terms and conditions: (https://www.beilsteinjournals.org/bjnano)