Abstract
Viscoelastic characterization of materials at the micro and the nanoscale is commonly performed with the aid of force–distance relationships acquired using atomic force microscopy (AFM). The general strategy for existing methods is to fit the observed material behavior to specific viscoelastic models, such as generalized viscoelastic models or powerlaw rheology models, among others. Here we propose a new method to invert and obtain the viscoelastic properties of a material through the use of the Ztransform, without using a model. We present the rheological viscoelastic relations in their classical derivation and their zdomain correspondence. We illustrate the proposed technique on a model experiment involving a traditional rampshaped force–distance AFM curve, demonstrating good agreement between the viscoelastic characteristics extracted from the simulated experiment and the theoretical expectations. We also provide a path for calculating standard viscoelastic responses from the extracted material characteristics. The new technique based on the Ztransform is complementary to previous modelbased viscoelastic analyses and can be advantageous with respect to Fourier techniques due to its generality. Additionally, it can handle the unbounded inputs traditionally used to acquire force–distance relationships in AFM, such as ramp functions, in which the cantilever position is displaced linearly with time for a finite period of time.
Introduction
Atomic force microscopy (AFM) is a prominent technique for investigating material properties at the micro and the nanoscale [13], within which a wide variety of instruments, probes, and analysis techniques have been developed to attempt meaningful material property extraction [411]. With regards to viscoelasticity, efforts that incorporate classical viscoelastic theory [1216] rely on force–distance curves [1726], which describe the dependence of the probe–sample interaction force with respect to the probe–surface distance for a particular location on the sample [27]. Force–distance analysis provides direct information on the force and indentation history with respect to time, which makes it appropriate for viscoelastic material property inversion. Existing inversion methods are based on viscoelastic models, which have been developed to describe specific kinds of relaxation behavior in the material. For example, in some cases, the analysis involves an assumption regarding discrete characteristic relaxation timescales in the material, which are represented using spring–dashpot models [20,24], while in other cases a continuous distribution of characteristic times is assumed via powerlaw rheology models [17,22]. Regardless of the model chosen, the strategy encompasses fitting the properties implied by the model to the experimental force–indentation data.
In order to enable a new route to viscoelastic material property inversion, which is complementary to previous strategies, here we propose a paradigmatically different approach that is not based on the choice of a model. Instead of fitting force–distance data to specific functions dictated by approximate models, we transform the experimental information using the Ztransform mathematical technique into the socalled zdomain, which is analogous to the Laplace domain, but applicable to discrete finite signals [28]. This enables the extraction of the viscoelastic transfer functions of the material, bypassing the need for any viscoelastic model assumptions (see Figure 1). These transfer functions are the viscoelastic relaxance and retardance, where the relaxance describes the time response of a viscoelastic material to a unit impulse excitation (Dirac delta function) of strain and the retardance describes the response of the material to an impulsive excitation of stress.
Our proposed modelfree approach follows the spirit of related Fourierbased methods, where viscoelastic material extraction has been achieved under some limited circumstances [2931]. However, in this work, we exploit the advantages of the Ztransform, which we believe is more appropriate for the analysis of viscoelasticity. Guided by the same motivation that has led the classical rheology field to rely on the more general Laplace transform technique for analytical treatments, we rely on the Ztransform to accomplish viscoelastic extraction for finite discrete experimental data. The selection of the Ztransform over the more widely exploited discrete Fourier transform (DFT) will become evident in later sections of the manuscript, where it is shown that our approach has enough generality to obtain meaningful material property information in a modelfree fashion, without the periodicity and convergence constraints associated with the DFT, which is better suited for harmonic and/or steadystate excitations. Although the proposed technique can be utilized to analyze different types of experiments, in this work we demonstrate it using classical force–distance experiments (see Figure 2), for which we carry out the force analysis in the complex domain (force spectroscopy).
Theoretical Background
It is well known that the behavior of viscoelastic materials is historydependent, as a result of which, the stress–strain relationships governing their deformation are functionals (not functions). More specifically, the stress at a given instant depends on the total previous history of strain and vice versa [13]. This history dependence is often expressed in the form of convolution integrals:
where Q(t) and U(t) are known as relaxance and retardance, respectively. As already stated, relaxance and retardance describe the time response of a viscoelastic material to a unit impulse excitation (Dirac delta function) of either strain or stress, respectively [13]. Theoretically, knowledge of Q(t) or U(t) fully characterizes the viscoelastic behavior of the material, so we refer to them as “source” functions.
Linear viscoelasticity exploits the use of the Laplace transform, whose advantage becomes evident upon transforming the previous equations into the Laplace domain and observing the simplicity afforded by this treatment. It is a wellknown property of the Laplace transform that convolutions in the time domain, such as the righthand side of Equation 1 and Equation 2, transform into simple multiplications in the Laplace domain (sdomain) [32,33]. Thus, Equations 1 and 2 transform as follows:
From Equation 3 and Equation 4, we can write the relaxance (Q(s)) and retardance (U(s)) as operators or transfer functions:
From Equation 5 and Equation 6, it is clear that the relaxance is a transfer function with which the stress can be calculated using a given strain input. Likewise, by using the retardance, strain can be calculated for a specific stress input. Clearly also, the relaxance and retardance are inverses of one another in the Laplace domain.
Unit impulse excitations (Dirac delta functions) are mathematically convenient for defining material behavior but are experimentally impractical. Therefore, rheologists normally use a different type of “standard” excitations. The material responses to these standard excitations are known as “standard viscoelastic responses” and are used widely to characterize viscoelastic materials. A very common standard excitation is the harmonic excitation. For example, a harmonic stress input in the time domain that is governed by the angular frequency ω would be represented as:
which can also be expressed using sine and cosine functions. The corresponding standard viscoelastic response to this harmonic stress excitation in the steady state is calculated by multiplying the appropriate transfer function (the retardance in this case, since the input is stress and the output is strain) by the above stress input. We perform this operation for the case where the Laplace variable s equals iω:
In the above expression, J*(ω) is the complex compliance, the sinusoidal steadystate strain response to a unit sinusoidal steadystate stress input [12]. From the above equation, it can also be seen that the complex compliance J*(ω) is the Fourier transform of the retardance U(t). Also, it is a special case of the Laplacetransformed retardance when the complex variable s is regarded as purely imaginary:
Furthermore, this steadystate complex compliance can be separated into its real and imaginary components:
where J′(ω) and J″(ω) are the storage and loss compliance, respectively. The wellknown storage and loss moduli can be defined in a similar fashion for the case where a sinusoidal input strain is applied to the material. It is very important to recall the fact that Equations 7–10 and the frequencydependent quantities included in them refer only to the steadystate case of sinusoidal excitation. Most types of AFM characterization do not apply sinusoidal stresses or strains to the material. This is clear with regards to the acquisition of a quasistatic force–distance curve, where the cantilever position above the sample follows a ramp function. In the case of intermittentcontact methods (e.g., tappingmode AFM), the cantilever tip oscillates nearly sinusoidally, but since tip–sample contact is intermittent, the sample does not experience purely sinusoidal stresses and strains.
To analyze the case of an AFM tip penetrating a viscoelastic surface we need an equation relating force with sample penetration (indentation). Equation 1 and Equation 2 relate stress and strain but do not consider the geometrical aspects of our boundary value problem. Invoking the correspondence principle for the case of a parabolic tip indenting a viscoelastic halfspace yields the following force–distance relationship for the time domain [14,15]:
where h is the sample indentation, F is the probe–sample force, and R is the radius of the indenter. Transforming this expression into the Laplace domain we obtain:
In a force–distancecurve experiment, we have access to both force and indentation history, such that we can, in principle, solve for U(t). We have previously used this strategy to obtain the timedomain representation of U(t) within a specific viscoelastic model [19,20]. In the interest of maintaining generality, we have used the generalized Voigt model, for which the retardance has the following form (the generalized Voigt and Maxwell–Wiechert models are discussed in Supporting Information File 1):
where δ(t) is the Dirac delta function, J_{g} is the glassy compliance, and J_{n} and τ_{n} are the compliance and retardation time, respectively, of the nth Voigt unit in the generalized Voigt model. Physically, these retardation times are the characteristic times at which the molecular rearrangements occur within the viscoelastic material. Replacing the above expression in the time convolution of Equation 11 and simplifying, we obtain the equation we have previously used as the basis for extracting the viscoelastic model parameters and corresponding properties using nonlinear leastsquares optimization in our previous work:
An alternate strategy is to turn to the Laplace domain equation, working with the special case when s = iω, such that the equations are handled more conveniently in Fourier space instead of Laplace space. This seems convenient because, having force and indentation time history, one could envision transforming both of them into Fourier space [29] to obtain the complex compliance. Specifically, letting s = iω in Equation 12 and rearranging we obtain:
where J*(ω) is the Fourier transform of the retardance, also sometimes denoted as (ω). The convenience is also apparent as many computer packages can efficiently compute the Fourier transform of discrete (experimental) data through the efficient FFT algorithm [34]. However, as we have already pointed out, this only defines the steadystate harmonic response of the material and is therefore not applicable to transient or nonsteadystate applications, such as a quasistatic force–distancecurve experiment. This approach is well suited for harmonic excitations but is inappropriate for various types of common experimental inputs such as step excitations or ramp excitations, which is why the more general Laplace transform is widely used in rheology. Methodologies have been developed to circumvent the above restriction of the Fourier transform by applying time derivatives to the experimental quantities and making use of its wellknown mathematical properties [29], but although the approaches are theoretically feasible, the existence of noise in AFM signals generally precludes their successful experimental implementation.
ZTransform Approach
In this paper, we exploit the wider generality of the Laplace transform to deliver a modelfree method to extract the viscoelastic retardance and relaxance of the material from AFM force curves. We focus on the Ztransform, which can be regarded as the discretetime counterpart to the Laplace transform, and which is, therefore, more general than the discretetime Fourier transform (a more detailed description of these integral transforms can be found in Supporting Information File 1). In the Laplace transform, a function is mapped into the complex domain via the variable s = α + iω, whereas in the Ztransform the function is mapped into the complex domain via the variable z = e^{α}e^{i}^{ω} (note that the real part of the complex variable s is customarily denoted by σ, but here we represent it with α to prevent any confusion with the stress, which is also denoted by σ). Figure 3 provides a comparison between the splane of the Laplace transform and the zplane of the Ztransform. Based on the relationship between the variables s and z and disregarding (for now) the fact that the Ztransform is a discrete transform, we see that the zplane is a polar version of the splane, whereby circles on the zplane correspond to vertical lines on the splane. For example, the imaginary axis of the splane, which corresponds to the Fourier transform (s = iω, with α = 0), maps to the unit circle on the zplane. It is especially easy to relate the unit circle of the zplane to the discrete Fourier transform. On the unit circle, the angular position of a point can be directly assigned to a frequency value in the discrete Fourier transform. The arc between 0° and −180° (from z = 1 to z = −1, clockwise) corresponds to frequencies between 0 Hz and the Nyquist frequency (half of the sampling frequency). The complementary arc between 0° and +180° (from z = 1 to z = −1, counterclockwise) corresponds to frequencies between 0 Hz and the negative of the Nyquist frequency. Other vertical lines on the splane correspond to nonunit circles on the zplane.
The use of the Ztransform reflects the fact that the experimental data obtained from an AFM experiment always consist of discrete signals. Furthermore, in quasistatic force–distance experiments, the deformation imposed on the material consists of a nonperiodic, nonsteadystate excitation (e.g., a ramp function). Therefore, the most appropriate transform is one that is discrete and reflects the capabilities of the Laplace transform [3537].
Reflecting the fact that the stress and strain now consist of the discrete signals σ[n] and ε[n], respectively, we can write discrete relationships that are analogous to the convolutions presented in Equation 1 and Equation 2:
where N corresponds to the number of points in the signal. As previously done for writing Equation 3 and Equation 4, we make use of the convolution properties of the Ztransform, which are similar to those of the Laplace transform, to obtain:
where all of the above variables are transformed variables in the zplane. We now write Q(z) and U(z) as operators or transfer functions as:
We see that the relaxance is a transfer function in the zplane with which stress can be calculated for a given strain input. Likewise, by using the retardance, strain can be calculated for a given stress input. The Ztransforms of the observables are calculated using the definition of the Ztransform:
Following Equation 12, which reflects the geometry of a spherical indentation experiment, we can write:
Similar to our previous methods, these expressions allow for the calculation of the desired transfer functions from experimental data.
Results and Discussion
Calculation of the source functions in the zdomain
First, to illustrate the proposed technique, we have calculated the retardance of a simulated material. The material is modeled using a generalized Voigt model with a single Voigt unit in series with a spring. Details about the simulation and the material parameters can be found in the Methods section. Linear ramp stress is applied to the material, from which the corresponding strain is calculated. This calculation gives the stress vs strain as a timedependent array, analogous to an AFM force–distance experiment, which is discussed below. The stress and the strain time series are then transformed into the zdomain using Equation 22. Finally, in the zdomain, we divide the strain by the stress to obtain the retardance, as prescribed by Equation 21. For comparison, we have also evaluated the retardance from the theoretical equation (Equation 29 in the Methods section) to compare it with our calculated retardance and both are plotted in Figure 4a (calculated) and Figure 4b (theoretical) in the zdomain. Apart from the unit circle of the zdomain, which corresponds to the Fourier transform of the retardance, the calculated and theoretical retardance are in good agreement (the discrepancy at the unit circle will be discussed below). For zdomain points outside of the unit circle, the error between the simulated retardance and the theoretical prediction is less than 0.25%. This is easier to visualize in Figure 4d–f, which do not include the unit circle.
One of the most significant advantages of working in the complex domain is the ease of inversion of the operators. Unlike the time domain where the inversion of retardance into relaxance and vice versa is challenging and requires either fitting routines for specific models [3841] or convoluted methods [42], in the zdomain, once the retardance is calculated, the relaxance is available through simple inversion of numerator and denominator, as stated in Equation 20 and Equation 21. This inversion is illustrated in Figure 5 for the data provided in Figure 4. Furthermore, although we have used a viscoelastic model to generate the material behavior for our simulation, the calculated retardance and relaxance (Figure 4 and Figure 5) are directly derived from the response of the material and are modelfree. Therefore, they do not inherently carry model assumptions or limitations. One may use these material operators in a modelfree fashion to predict the material response to an arbitrary stress or strain input, as will be discussed later. Alternatively, it is also possible to fit the calculated retardance (or relaxance) to a viscoelastic model and calculate relevant material constants, characteristic times, and/or loss and storage moduli, as suggested in the bottomright corner of Figure 1.
Relationship between the real axis of the zdomain and the timescale of the response
An interesting feature of the zplane is its real axis. Consider the point z = 1, which implies e^{α}e^{i}^{ω} = 1, or α = 0 and ω = 0. This point thus corresponds to the zero frequency, or DC component of the Fourier transform. Likewise, other points on this axis correspond to the zerofrequency components of a modified Fourier transform. For example, the point z = b, where b is a real number different than unity, requires ω = 0 (similar to z = 1) but corresponds to a nonzero value of α, such that b = e^{α}, where α is positive for b > 1 and negative for b < 1. If, say, b < 1, the point z = b corresponds to the zero frequency of the modified Fourier transform e^{−α}e^{i}^{ω}, which represents a decaying (damped) oscillation. Thus, depending on their position on the real axis, the points z = b correspond to either increasing or decreasing exponentials with different time constants, α = ln(z). For our purposes, these time constants can be thought of as different timescales for material behavior. For example, at z = 1 the time constant is zero, hence our exponential yields a constant steadystate DC component. Increasingly larger z values can be thought of as material responses for shorter and shorter timescales. Figure 6 shows the retardance of our example material plotted along the real axis. For this material, the theoretical retardance for z = 1 (this is the verticalaxis intercept, since the horizontal axis is ln(z) and ln(1) = 0) is the same as the equilibrium compliance of the material, also known as the creep compliance, for which the simulation provides a fair approximation according to Figure 6. As the value of z increases along the real axis, the value of the retardance approaches the theoretical glassy compliance, which refers to the material’s infinitely shorttimescale response (i.e., instantaneous response), and which is approximated very well by the simulation. We can thus obtain two very important properties of the material quite easily. It is of course not unexpected to obtain the glassy and equilibrium compliances from the real axis, as inspection of Equation 29 (Methods section) shows that those are the expected limits of the retardance:
Limitations of Fourier techniques
The Fourier representation of the retardance and the relaxance of the material is commonly used in rheology, where the real components of the operators are referred to as storage modulus or storage compliance, respectively, and the imaginary components as loss modulus or loss compliance, respectively. Although directly obtaining the Fourier representation of the source operators is possible with the proposed method when harmonic or bounded inputs are used, this is not the case for the classical force–distance curve approach, which is based on a ramp function for which the Fourier integral does not converge. Although mathematical tricks can be used for a welldefined function [29,43], additional experimental complications emerge since the signals in an experiment are discrete. When dealing with a bounded function, it can be reasonable to treat the sampling time window as periodic, such that the use of the discrete Fourier transform is warranted. However, this is not appropriate for an unbounded function, which does not possess a finite bandwidth, and for which assuming periodicity leads to signal aliasing. From this consideration we see that the discrete Fourier transform does not necessarily represent the continuous analytical Fourier transform [30], and therefore, a viscoelastic source function calculated using the discrete Fourier transforms of the observables for experiments with unbounded input functions does not correspond to the theoretical harmonic response of the material (Equations 7–10). This is, in fact, the underlying reason for the discrepancy between the theoretical retardance and the calculated retardance in Figure 4 for the unit circle (the Ztransform on the unit circle corresponds to a discrete Fourier transform). In Supporting Information File 1, Figure S6, it is demonstrated that this discrepancy also occurs for the FFT algorithm.
Since nonunit circles of the zplane correspond to the modified Fourier transforms of the system, which contain a timedependent exponential component in addition to their harmonic parts, they can properly represent the characteristics of an experiment with unbounded input. The modified Fourier transforms can properly handle nonsteadystate, nonbounded, nonperiodic systems, accommodating transient and nonharmonic responses. Figure 7 compares the simulated calculated retardance with the corresponding theoretical values demonstrating good agreement.
As already stated, the storage and loss compliance are not directly accessible from the force–distance experiments, for which the input is unbounded (see also Supporting Information File 1, Figure S6). However, it may be possible to approximate them in some cases by inspecting the real and imaginary components of the neighboring modified Fourier transform regions (i.e., nonunit circles on the zplane which are close to the unit circle). An example is shown in Figure 8 for a circle on the zplane with a radius r = 1.002. The error between the estimated storage and loss compliances and the analytical result is small for high frequencies and larger for low frequencies, although in all cases it remains within the order of magnitude of the desired quantities. A much more detailed mathematical analysis of the differences between the Fourier and modified Fourier loss and storage compliances is provided in the Supporting Information File 1. For the generalized Voigt model, the modified Fourier transform version (nonzero time constant) of the magnitudes of the storage and the loss compliance plotted against frequency are reduced by a factor of with respect to the Fourier result (zero time constant). Similarly, the width of the peaks is increased by a factor of with respect to the Fourier result and we can observe a compliance behavior that is more spread out over different frequencies.
As illustrated in Figure 8, the loss and storage compliance and their estimation from the modified Fourier transform converge for high frequencies. Therefore, it is possible to accurately estimate the glassy compliance of a material and its immediate response from a force–distance AFM experiment. Furthermore, the position of the local maxima of the loss compliance with respect to frequency changes slightly between the Fourier transform and the modified Fourier transform for the generalized Voigt and Maxwell–Wiechert models, and is given by Hz (a mathematical derivation and a discussion are provided in chapter 6 of Supporting Information File 1). Therefore, the characteristic times can be accurately pinpointed from the modified Fourier transform estimation. Once the characteristic times of the material are estimated, the loss and storage compliances can be calculated from the modified Fourier transforms in light of the peak geometry correction factor discussed above, . For this calculation Δt is a known experimental parameter, and r is chosen by the researcher.
Demonstration with AFM contact mechanics
So far, we have demonstrated our method for stress–strain inputs using the generalized Voigt model. However, in AFM experiments one observes the deflection of the AFM cantilever as a function of the cantilever base position instead of directly observing stress and strain. From the available observables, one can calculate the tip–sample force and the indentation. In order to account for the AFM probe and sample geometry and the nature of the corresponding observables, it is necessary to invoke the correspondence principle, through which Equation 11 and Equation 24 are derived. We have simulated an AFM experiment (see Methods section for details), calculating the retardance and comparing it with the theoretical behavior, as shown in Figure 9. As before, we observe a good agreement between simulation and theory, which suggests that the method is also suitable for AFM analysis. Retardance profiles for different circles of the zplane and along the real axis, similar to the data of Figure 6 and Figure 7, can be found in Figures S7 and S8 in Supporting Information File 1.
Calculation of responses to standard or arbitrary inputs through the inverse Ztransform
So far we have only discussed the behavior of the retardance and relaxance in the zdomain. However, we recall that these operators define the reaction of the material to stress or strain, noting that with the proposed method the representation is modelfree. The next logical question is what to do with these operators. Naturally, one option is to fit an established viscoelastic model to them, such as the generalized Voigt or Maxwell–Wiechert models, or any other suitable model. As long as there is a complex domain correspondence between the zplane and the domain of the model, any model can be fit, which yields the parameters of the model. These parameters can then be used to express, for example, compliances, moduli, and characteristic times in analytical form. By extension, one can also obtain the storage and loss moduli. Alternatively, one may continue with a modelfree approach, using the operators as they are to predict material responses to specific inputs. For this, all one needs to do is transform the userdefined input into the zdomain and multiply with the appropriate operator (calculated relaxance or retardance). The result of this multiplication is the material response in the zdomain, from which the timedomain response can be obtained via the inverse Ztransform. Figure 10 shows the result of calculating the creep (strain) response of our model material from the calculated retardance for a unit stress input. The inverse Ztransform is calculated by taking a counterclockwise closed integral along a contour in the zplane (see Equation 32 in the Methods section) [35]. In this work, we have calculated the integral numerically around a nonunit circle, counterclockwise, after evaluating the retardance in the zplane, for convenience.
Discussion and future work
So far the proposed method seems to be quite promising, but there is one very important consideration that needs to be highlighted, namely computer precision in the calculations. Specifically, modern computational coding languages use variable precision. For example, for a number that is of order unity, the computational precision is ≈10^{−16}, but for a number that is of order 10^{100} the precision is ≈10^{84}. Therefore, in calculations with large numbers, one can accumulate a massive amount of error, and there are often situations where positive and negative numbers that should cancel each other do not do so due to lack of precision. Unfortunately, this issue plagues the Ztransform due to its numerical nature. While calculating the Ztransform values on the zplane, one calculates the zvalue’s power −n, where n scales with the number of sample points in the signal. Likewise, while calculating the inverse Ztransform, one raises z to the nth power. During the calculation of the Ztransform, numbers raised to very large negative powers converge to zero, which although not catastrophic, does result in loss of accuracy. On the other hand, while calculating the inverse Ztransform, calculating very large powers can inaccurately diverge to infinity. This can place practical limits on the sample size (number of points in the signal) and may also require very mindful coding.
There are also areas of future work concerning the methods presented here. We have focused on the nonunit circles of the zplane, corresponding to modified Fourier transforms of the operators, where each circle represents a combination of a harmonic response with a timedependent exponential coefficient having a specific time constant. As previously stated, these exponential coefficients represent nonsteadystate behaviors and are crucial to our method. One area of further investigation concerns the relationship between the time constants (or zdomain radii) and the time scale of the experiment. Since these timedependent coefficients lay along the real axis of the zdomain (without any harmonic component), the real axis of the zdomain contains very important information concerning nonharmonic and nonsteadystate behavior in such cases. Nonunit circles with different radii (and hence different time constants) may be more or less important for different operations with different time scales. Furthermore, the nonunit circles of the zdomain also contain information about the harmonic behavior of the material (although under damped conditions). Thus, the storage and loss of energy by the material and its relationship to the modified harmonic components should also be investigated further.
One final area of future work concerns the effect of electronic or thermal noise on the quality of the viscoelastic analysis. Generally speaking, analyses in the complex domain can handle noise more robustly than analyses in the time domain, especially when specific features are expected in the plots of the viscoelastic functions. For example, the fact that in the Kelvin–Voigt model, nonzerocentered peaks are not expected in the storage compliance (see Figure 11) could be used to discriminate noise from legitimate material behavior. One can envision the development of smart algorithms that “clean up” the viscoelastic functions based on their expected behavior. However, this is clearly more challenging in our modelfree framework, where the features of the material behavior may not be known a priori. Further research is encouraged on this topic.
Conclusion
A novel method for obtaining the viscoelastic properties of a material using atomic force spectroscopy has been proposed and demonstrated computationally. The method utilizes Ztransform techniques and yields modelfree viscoelastic information, such as the material retardance and relaxance, as well as standard viscoelastic responses and information on material behavior at different timescales. The method has advantages over discrete Fourier transform methods, in that it can handle unbounded and nonperiodic signals, such as the inputs that are traditionally used in quasistatic AFM force–distance experiments. Furthermore, the acquisition of a transfer function that defines the linear viscoelastic behavior of the material enables the generation of any standard and nonstandard viscoelastic material response from it once the input is defined by the user. Although the method provides a way to perform modelfree viscoelastic analysis, it is also possible to fit the transfer functions to specific models in order to obtain parametrized analytical expressions for them.
Methods
The viscoelasticity of the material is simulated by using a single Voigt unit plus a residual spring within the generalized Voigt model (see Supporting Information File 1, Figure S2, Equation S42), which contains one characteristic retardation time τ. The material properties for this model are provided in Table 1.
To calculate the material retardance and relaxance in the zdomain, a ramp input stress is first defined (Equation 27) and the corresponding strain response is then calculated using Equation 28 below, which gives the theoretical response for the chosen model [13]. The duration of the simulated experiment is 0.1 s and the data is discretized using a timestep of 0.00001 s:
The stress and strain signals are subsequently transformed into the zdomain (see Equation S80, Supporting Information File 1). One can visualize the process by considering the signals as a succession of timeshifted delta functions, which can be easily transformed into the zdomain. The theoretical retardance is evaluated using Equation 29 (see Supporting Information File 1 for its derivation):
A spherical contact, appropriate for an AFM experiment, was simulated using the same material parameters, with a tip radius of 10 nm, which is common in AFM (See also Figure 2). Analogous to the previous simulation, a ramp indentation with a slope of 1 nm/s was applied to the material for 0.1 s. However, this time we have used Simulink to evaluate Equation 12. The equation was solved using a fourthorder Runge–Kutta solver with an integration timestep of 10^{−9} s, which was downsampled to a timestep of 10^{−5} s to simulate the experimental AFM signal. Note that in the case of an AFM experiment the user does not generally prescribe a variation of indentation over time. Instead, one prescribes the timedependent displacement of the cantilever base towards and away from the sample while measuring the force (via the deflection). However, the simpler simulation described here is appropriate to evaluate the Ztransform methodology because the data acquired through the AFM experiment needs to be transformed into force vs indentation information before it can be used as input for a viscoelastic inversion analysis [20].
To calculate the creep response, stress and strain were evaluated for a ramp input as in Equation 27 and Equation 28, for a shorter time period with a larger time step. The strain was calculated for 0.045 s with a timestep of 0.001 s. The subsequently calculated retardance was multiplied with a unit step input in the zdomain:
The above corresponds to the creep response in the zdomain. The unit step input in the zdomain can be written as:
To calculate the creep response in the time domain, the zdomain creep response was obtained by numerically integrating the contour along the circle with R = 1.14 counterclockwise. The theoretical creep response was evaluated using Simulink with a fourthorder Runge–Kutta solver using a timestep of 10^{−6} s, taking Equation 6 as the transfer function, in the form of Equation S42, with the parameters from Table 1. As in the previous case, the results were downsampled for demonstration purposes.
The inverse Ztransform of an arbitrary function in the zdomain, X(z), can be calculated using the following contour integral [35]:
Supporting Information
The Supporting Information features derivations of generalized Voigt and Maxwell–Wiechert models in the Laplace and the zdomains; a brief overview of Fourier, modified Fourier, Laplace, and Ztransforms; additional data illustrating the misrepresentation of the system in the Fourier domain and additional data from AFM simulations; and, finally, a detailed analysis of loss and storage as a function of frequency and their estimation from modified Fourier transforms.
Supporting Information File 1: Additional methodical data.  
Format: PDF  Size: 1.6 MB  Download 
References

Binnig, G.; Quate, C. F.; Gerber, C. Phys. Rev. Lett. 1986, 56, 930–933. doi:10.1103/physrevlett.56.930
Return to citation in text: [1] 
García, R.; Pérez, R. Surf. Sci. Rep. 2002, 47, 197–301. doi:10.1016/s01675729(02)000778
Return to citation in text: [1] 
Garcia, R.; Herruzo, E. T. Nat. Nanotechnol. 2012, 7, 217–226. doi:10.1038/nnano.2012.38
Return to citation in text: [1] 
Herruzo, E. T.; Perrino, A. P.; Garcia, R. Nat. Commun. 2014, 5, 3126. doi:10.1038/ncomms4126
Return to citation in text: [1] 
Amo, C. A.; Perrino, A. P.; Payam, A. F.; Garcia, R. ACS Nano 2017, 11, 8650–8659. doi:10.1021/acsnano.7b04381
Return to citation in text: [1] 
Rico, F.; RocaCusachs, P.; Gavara, N.; Farré, R.; Rotger, M.; Navajas, D. Phys. Rev. E 2005, 72, 021914. doi:10.1103/physreve.72.021914
Return to citation in text: [1] 
LópezGuerra, E. A.; Banfi, F.; Solares, S. D.; Ferrini, G. Sci. Rep. 2018, 8, 7534. doi:10.1038/s41598018258284
Return to citation in text: [1] 
LópezGuerra, E. A.; Somnath, S.; Solares, S. D.; Jesse, S.; Ferrini, G. Sci. Rep. 2019, 9, 12721. doi:10.1038/s41598019491041
Return to citation in text: [1] 
Radmacher, M.; Tillmann, R. W.; Gaub, H. E. Biophys. J. 1993, 64, 735–742. doi:10.1016/s00063495(93)814334
Return to citation in text: [1] 
Yuya, P. A.; Hurley, D. C.; Turner, J. A. J. Appl. Phys. 2008, 104, 074916. doi:10.1063/1.2996259
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] 
Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley & Sons, 1980.
Return to citation in text: [1] [2] 
Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior; Springer: Berlin, Heidelberg, 1989. doi:10.1007/9783642736025
Return to citation in text: [1] [2] [3] [4] 
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] 
Graham, G. A. C. Int. J. Eng. Sci. 1965, 3, 27–46. doi:10.1016/00207225(65)900182
Return to citation in text: [1] [2] 
Ting, T. C. T. J. Appl. Mech. 1966, 33, 845–854. doi:10.1115/1.3625192
Return to citation in text: [1] 
de Sousa, J. S.; Santos, J. A. C.; Barros, E. B.; Alencar, L. M. R.; Cruz, W. T.; Ramos, M. V.; Mendes Filho, J. J. Appl. Phys. 2017, 121, 034901. doi:10.1063/1.4974043
Return to citation in text: [1] [2] 
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] 
LópezGuerra, E. A.; Shen, H.; Solares, S. D.; Shuai, D. Nanoscale 2019, 11, 8918–8929. doi:10.1039/c8nr10287b
Return to citation in text: [1] [2] 
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] 
Garcia, P. D.; Guerrero, C. R.; Garcia, R. Nanoscale 2017, 9, 12051–12059. doi:10.1039/c7nr03419a
Return to citation in text: [1] 
Garcia, P. D.; Guerrero, C. R.; Garcia, R. Nanoscale 2020, 12, 9133–9143. doi:10.1039/c9nr10316c
Return to citation in text: [1] [2] 
Rajabifar, B.; Jadhav, J. M.; Kiracofe, D.; Meyers, G. F.; Raman, A. Macromolecules 2018, 51, 9649–9661. doi:10.1021/acs.macromol.8b01485
Return to citation in text: [1] 
Chyasnavichyus, M.; Young, S. L.; Tsukruk, V. V. Langmuir 2014, 30, 10566–10582. doi:10.1021/la404925h
Return to citation in text: [1] [2] 
Parvini, C. H.; Saadi, M. A. S. R.; Solares, S. D. Beilstein J. Nanotechnol. 2020, 11, 922–937. doi:10.3762/bjnano.11.77
Return to citation in text: [1] 
Forstenhäusler, M.; LópezGuerra, E. A.; Solares, S. D. Facta Univ., Ser.: Mech. Eng. 2021, 19, 133–153. doi:10.22190/fume200920009f
Return to citation in text: [1] 
Garcia, R. Chem. Soc. Rev. 2020, 49, 5850–5884. doi:10.1039/d0cs00318b
Return to citation in text: [1] 
Oppenheim, A. V.; Willsky, A. S.; Nawab, S. H. Signals & Systems, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, USA, 1997.
Return to citation in text: [1] 
Evans, R. M. L.; Tassieri, M.; Auhl, D.; Waigh, T. A. Phys. Rev. E 2009, 80, 012501. doi:10.1103/physreve.80.012501
Return to citation in text: [1] [2] [3] [4] 
Aspden, R. M. J. Phys. D: Appl. Phys. 1991, 24, 803–808. doi:10.1088/00223727/24/6/002
Return to citation in text: [1] [2] [3] 
Holly, E. E.; Venkataraman, S. K.; Chambon, F.; Henning Winter, H. J. NonNewtonian Fluid Mech. 1988, 27, 17–26. doi:10.1016/03770257(88)800028
Return to citation in text: [1] 
Gardner, M. F.; Barnes, J. L. Transients in Linear Systems Studied by the Laplace Transformation; John Wiley & Sons, 1956.
Return to citation in text: [1] 
Stroud, K. A.; Booth, D. J. Advanced Engineering Mathematics; Palgrave Macmillan, 2011.
Return to citation in text: [1] 
Cooley, J. W.; Tukey, J. W. Math. Comput. 1965, 19, 297. doi:10.1090/s00255718196501785861
Return to citation in text: [1] 
Oppenheim, A. V.; Schafer, R. W. Digital Signal Processing; PretinceHall: Englewood Cliffe, NJ, USA, 1974.
Return to citation in text: [1] [2] [3] 
McClellan, J. H.; Schafer, R. W.; Yoder, M. A. Signal Processing First; Pearson\Prentice Hall: Upper Saddle River, NJ, USA, 2003.
Return to citation in text: [1] 
Smith, S. W. The Scientist & Engineer’s Guide to Digital Signal Processing, 1st ed.; California Technical Publishing: San Diego, CA, USA, 1997.
Return to citation in text: [1] 
Schapery, R. A.; Park, S. W. Int. J. Solids Struct. 1999, 36, 1677–1699. doi:10.1016/s00207683(98)000602
Return to citation in text: [1] 
Park, S. W.; Schapery, R. A. Int. J. Solids Struct. 1999, 36, 1653–1675. doi:10.1016/s00207683(98)000559
Return to citation in text: [1] 
Mead, D. W. J. Rheol. 1994, 38, 1769–1795. doi:10.1122/1.550526
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, D.C., USA, 2018.
Return to citation in text: [1] 
Sorvari, J.; Malinen, M. Int. J. Solids Struct. 2007, 44, 1291–1303. doi:10.1016/j.ijsolstr.2006.06.029
Return to citation in text: [1] 
Kusse, B. R.; Westwig, E. A. Mathematical Physics; WileyVCH: Weinheim, Germany, 2006. doi:10.1002/9783527618132
Return to citation in text: [1]
35.  Oppenheim, A. V.; Schafer, R. W. Digital Signal Processing; PretinceHall: Englewood Cliffe, NJ, USA, 1974. 
1.  Binnig, G.; Quate, C. F.; Gerber, C. Phys. Rev. Lett. 1986, 56, 930–933. doi:10.1103/physrevlett.56.930 
2.  García, R.; Pérez, R. Surf. Sci. Rep. 2002, 47, 197–301. doi:10.1016/s01675729(02)000778 
3.  Garcia, R.; Herruzo, E. T. Nat. Nanotechnol. 2012, 7, 217–226. doi:10.1038/nnano.2012.38 
19.  LópezGuerra, E. A.; Shen, H.; Solares, S. D.; Shuai, D. Nanoscale 2019, 11, 8918–8929. doi:10.1039/c8nr10287b 
20.  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.  de Sousa, J. S.; Santos, J. A. C.; Barros, E. B.; Alencar, L. M. R.; Cruz, W. T.; Ramos, M. V.; Mendes Filho, J. J. Appl. Phys. 2017, 121, 034901. doi:10.1063/1.4974043 
18.  Efremov, Y. M.; Wang, W.H.; Hardy, S. D.; Geahlen, R. L.; Raman, A. Sci. Rep. 2017, 7, 1541. doi:10.1038/s41598017017843 
19.  LópezGuerra, E. A.; Shen, H.; Solares, S. D.; Shuai, D. Nanoscale 2019, 11, 8918–8929. doi:10.1039/c8nr10287b 
20.  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 
21.  Garcia, P. D.; Guerrero, C. R.; Garcia, R. Nanoscale 2017, 9, 12051–12059. doi:10.1039/c7nr03419a 
22.  Garcia, P. D.; Guerrero, C. R.; Garcia, R. Nanoscale 2020, 12, 9133–9143. doi:10.1039/c9nr10316c 
23.  Rajabifar, B.; Jadhav, J. M.; Kiracofe, D.; Meyers, G. F.; Raman, A. Macromolecules 2018, 51, 9649–9661. doi:10.1021/acs.macromol.8b01485 
24.  Chyasnavichyus, M.; Young, S. L.; Tsukruk, V. V. Langmuir 2014, 30, 10566–10582. doi:10.1021/la404925h 
25.  Parvini, C. H.; Saadi, M. A. S. R.; Solares, S. D. Beilstein J. Nanotechnol. 2020, 11, 922–937. doi:10.3762/bjnano.11.77 
26.  Forstenhäusler, M.; LópezGuerra, E. A.; Solares, S. D. Facta Univ., Ser.: Mech. Eng. 2021, 19, 133–153. doi:10.22190/fume200920009f 
29.  Evans, R. M. L.; Tassieri, M.; Auhl, D.; Waigh, T. A. Phys. Rev. E 2009, 80, 012501. doi:10.1103/physreve.80.012501 
12.  Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley & Sons, 1980. 
13.  Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior; Springer: Berlin, Heidelberg, 1989. doi:10.1007/9783642736025 
14.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
15.  Graham, G. A. C. Int. J. Eng. Sci. 1965, 3, 27–46. doi:10.1016/00207225(65)900182 
16.  Ting, T. C. T. J. Appl. Mech. 1966, 33, 845–854. doi:10.1115/1.3625192 
4.  Herruzo, E. T.; Perrino, A. P.; Garcia, R. Nat. Commun. 2014, 5, 3126. doi:10.1038/ncomms4126 
5.  Amo, C. A.; Perrino, A. P.; Payam, A. F.; Garcia, R. ACS Nano 2017, 11, 8650–8659. doi:10.1021/acsnano.7b04381 
6.  Rico, F.; RocaCusachs, P.; Gavara, N.; Farré, R.; Rotger, M.; Navajas, D. Phys. Rev. E 2005, 72, 021914. doi:10.1103/physreve.72.021914 
7.  LópezGuerra, E. A.; Banfi, F.; Solares, S. D.; Ferrini, G. Sci. Rep. 2018, 8, 7534. doi:10.1038/s41598018258284 
8.  LópezGuerra, E. A.; Somnath, S.; Solares, S. D.; Jesse, S.; Ferrini, G. Sci. Rep. 2019, 9, 12721. doi:10.1038/s41598019491041 
9.  Radmacher, M.; Tillmann, R. W.; Gaub, H. E. Biophys. J. 1993, 64, 735–742. doi:10.1016/s00063495(93)814334 
10.  Yuya, P. A.; Hurley, D. C.; Turner, J. A. J. Appl. Phys. 2008, 104, 074916. doi:10.1063/1.2996259 
11.  Attard, P. J. Phys.: Condens. Matter 2007, 19, 473201. doi:10.1088/09538984/19/47/473201 
14.  Lee, E. H.; Radok, J. R. M. J. Appl. Mech. 1960, 27, 438–444. doi:10.1115/1.3644020 
15.  Graham, G. A. C. Int. J. Eng. Sci. 1965, 3, 27–46. doi:10.1016/00207225(65)900182 
29.  Evans, R. M. L.; Tassieri, M.; Auhl, D.; Waigh, T. A. Phys. Rev. E 2009, 80, 012501. doi:10.1103/physreve.80.012501 
30.  Aspden, R. M. J. Phys. D: Appl. Phys. 1991, 24, 803–808. doi:10.1088/00223727/24/6/002 
31.  Holly, E. E.; Venkataraman, S. K.; Chambon, F.; Henning Winter, H. J. NonNewtonian Fluid Mech. 1988, 27, 17–26. doi:10.1016/03770257(88)800028 
13.  Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior; Springer: Berlin, Heidelberg, 1989. doi:10.1007/9783642736025 
28.  Oppenheim, A. V.; Willsky, A. S.; Nawab, S. H. Signals & Systems, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, USA, 1997. 
32.  Gardner, M. F.; Barnes, J. L. Transients in Linear Systems Studied by the Laplace Transformation; John Wiley & Sons, 1956. 
33.  Stroud, K. A.; Booth, D. J. Advanced Engineering Mathematics; Palgrave Macmillan, 2011. 
17.  de Sousa, J. S.; Santos, J. A. C.; Barros, E. B.; Alencar, L. M. R.; Cruz, W. T.; Ramos, M. V.; Mendes Filho, J. J. Appl. Phys. 2017, 121, 034901. doi:10.1063/1.4974043 
22.  Garcia, P. D.; Guerrero, C. R.; Garcia, R. Nanoscale 2020, 12, 9133–9143. doi:10.1039/c9nr10316c 
20.  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.  Chyasnavichyus, M.; Young, S. L.; Tsukruk, V. V. Langmuir 2014, 30, 10566–10582. doi:10.1021/la404925h 
13.  Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior; Springer: Berlin, Heidelberg, 1989. doi:10.1007/9783642736025 
35.  Oppenheim, A. V.; Schafer, R. W. Digital Signal Processing; PretinceHall: Englewood Cliffe, NJ, USA, 1974. 
36.  McClellan, J. H.; Schafer, R. W.; Yoder, M. A. Signal Processing First; Pearson\Prentice Hall: Upper Saddle River, NJ, USA, 2003. 
37.  Smith, S. W. The Scientist & Engineer’s Guide to Digital Signal Processing, 1st ed.; California Technical Publishing: San Diego, CA, USA, 1997. 
34.  Cooley, J. W.; Tukey, J. W. Math. Comput. 1965, 19, 297. doi:10.1090/s00255718196501785861 
29.  Evans, R. M. L.; Tassieri, M.; Auhl, D.; Waigh, T. A. Phys. Rev. E 2009, 80, 012501. doi:10.1103/physreve.80.012501 
13.  Tschoegl, N. W. The Phenomenological Theory of Linear Viscoelastic Behavior; Springer: Berlin, Heidelberg, 1989. doi:10.1007/9783642736025 
20.  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 
30.  Aspden, R. M. J. Phys. D: Appl. Phys. 1991, 24, 803–808. doi:10.1088/00223727/24/6/002 
35.  Oppenheim, A. V.; Schafer, R. W. Digital Signal Processing; PretinceHall: Englewood Cliffe, NJ, USA, 1974. 
42.  Sorvari, J.; Malinen, M. Int. J. Solids Struct. 2007, 44, 1291–1303. doi:10.1016/j.ijsolstr.2006.06.029 
29.  Evans, R. M. L.; Tassieri, M.; Auhl, D.; Waigh, T. A. Phys. Rev. E 2009, 80, 012501. doi:10.1103/physreve.80.012501 
43.  Kusse, B. R.; Westwig, E. A. Mathematical Physics; WileyVCH: Weinheim, Germany, 2006. doi:10.1002/9783527618132 
30.  Aspden, R. M. J. Phys. D: Appl. Phys. 1991, 24, 803–808. doi:10.1088/00223727/24/6/002 
38.  Schapery, R. A.; Park, S. W. Int. J. Solids Struct. 1999, 36, 1677–1699. doi:10.1016/s00207683(98)000602 
39.  Park, S. W.; Schapery, R. A. Int. J. Solids Struct. 1999, 36, 1653–1675. doi:10.1016/s00207683(98)000559 
40.  Mead, D. W. J. Rheol. 1994, 38, 1769–1795. doi:10.1122/1.550526 
41.  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, D.C., USA, 2018. 
© 2021 Uluutku et al.; licensee BeilsteinInstitut.
This is an Open Access article under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0). Please note that the reuse, redistribution and reproduction in particular requires that the author(s) and source are credited and that individual graphics may be subject to special legal provisions.
The license is subject to the Beilstein Journal of Nanotechnology terms and conditions: (https://www.beilsteinjournals.org/bjnano/terms)