Sensitivity of viscoelastic characterization in multi-harmonic atomic force microscopy

Quantifying the nanomechanical properties of soft-matter using multi-frequency atomic force microscopy (AFM) is crucial for studying the performance of polymers, ultra-thin coatings, and biological systems. Such characterization processes often make use of cantilever's spectral components to discern nanomechanical properties within a multi-parameter optimization problem. This could inadvertently lead to an over-determined parameter estimation with no clear relation between the identified parameters and their influence on the experimental data. In this work, we explore the sensitivity of viscoelastic characterization in polymeric samples to the experimental observables of multi-frequency intermodulation AFM. By performing simulations and experiments we show that surface viscoelasticity has negligible effect on the experimental data and can lead to inconsistent and often non-physical identified parameters. Our analysis reveals that this lack of influence of the surface parameters relates to a vanishing gradient and non-convexity while minimizing the objective function. By removing the surface dependency from the model, we show that the characterization of bulk properties can be achieved with ease and without any ambiguity. Our work sheds light on the sensitivity issues that can be faced when optimizing for a large number of parameters and observables in AFM operation, and calls for the development of new viscoelastic models at the nanoscale and improved computational methodologies for nanoscale mapping of viscoelasticity using AFM.

The lift motion denotes the motion of the cantilever at a position close to the surface. It provides a measure to compensate the contribution of long-range linear forces due to squeeze-film damping or electrostatic interactions. These effects are embedded in the linear transfer function of the so-called background forces χ BG [1]. From the measurements of ( d free , d eng , d lift ) we estimate the tip-sample nonlinear force at intermodulation frequencies by applying [2]: in which the last term corresponds to the background force compensation, with its associated linear transfer function defined by and approximated on the narrow frequency band with the polynomial [1]: The coefficients a and b come from the fit of Eq. (2) at the two drive frequencies (ω 1 , ω 2 ). In addition, we apply the following phase rotation to compensate the phase shift potentially caused by a time delay inherent to the processing equipment [3]: where ω c = 1 2 (ω 1 +ω 2 ) ≈ ω 0 , and F (c) ts denotes the tip-sample intermodulation components with the rotation coefficients (R 0 , R 1 ) adjusted such that arg(d eng (ω 1 )) = arg(d eng (ω 2 )) = 0: R 0 = arg(d eng (ω 1 )) − arg(d eng (ω 2 )) − arg(d eng (ω 1 )) ω 2 − ω 1 ω 1 (5) R 1 = arg(d eng (ω 2 )) − arg(d eng (ω 1 )) The phase equalization procedure defined by Eqs. (4)- (6) is also applied on the simulated components for comparison purposes. The driving force signal F d (t) used in the simulations is defined specifically for the experimental data considered in the study. In particular, it is estimated for the set of frequency, stiffness and quality factor of the first resonance of the cantilever f 0 , k and Q obtained from the thermal calibration. The excitation signal is obtained from the free motion frequency components as : with where the ω IM denotes the pulsation of intermodulation [4]. The time signals are simulated using the following dimensionless values: in which the displacement of reference is the amplitude of the engaged motion at the second drive frequency A = |d c | ω=ω 2 . The following dimensionless design parameters are considered in the numerical procedure: Thus, the equation of motion (Eq.(1) of the main manuscript) is in which the overbars are dropped for the sake of brevity. The time signals are computed by simulating Eq. (11) using a Runge-Kutta scheme. At low sample relaxation times τ s = η s /k s < 10 −3 , a scheme designed for stiff systems is employed (the ode23s function of Matlab is used, instead of the classical ode45 time integration solver). The signals for d c , d s and F ts are simulated on 8 ms, which corresponds to four intermodulation beatings since ∆ f = ω 2 −ω 1 2π = 500 Hz is applied in experiments. A zero initial condition for displacements and velocities is applied.
In order to convert the simulated tip-sample force signal F ts,sim from time to frequency domain at intermodulation frequencies, we extract two beat periods in steady state oscillations ( Fig. S2.1 (a-b)). The amplitude and phase components (| F ts,sim |, φ F ts ) of 32 intermodulation frequencies are estimated using a synchronous detection scheme [5]. Next, we use a sliding window with a length equal to one beat period as shown in Fig. S2.1 (c) and take 10 estimations of the phase and amplitude components. The estimations are then averaged to reduce numerical noise as shown in Fig. S2.1(d). Finally, the spectral components of the interaction force are stored in the same way as in experiments, in complex form like F ts,sim = | F ts,sim |e jφ F ts .
The objective function used for estimating the viscoelastic parameters is defined by [4,6]: In order to minimize Eq. (12) we use Levenberg-Marquardt algorithm [7] and combine it with nonlinear least squares (lsqnonlin function) in Matlab. This least-square minimization is performed with an iterative procedure which involves the computation of the partial derivatives (gradient) at each iteration. A parallel implementation on a small cluster was used to perform multiple minimization routines: approximately 10 nodes and 36 hours in total were needed to obtain the results shown in Fig. 5 of the main manuscript. We show in table S2.1 the lower and upper limit of parameter values defined for the optimization. These parameter ranges are deliberately wide because we assume that we have no prior knowledge of the material properties, except in the case of the probe height for which a first approximation can be extracted from the force quadrature curves (see section S2.C).
from windows averaged  Parameter Maximal value 100 10k 10k/ω 0 20k 20k/ω 0 45 In this section we discuss the use of a global optimization procedure for parameter estimation and further elaborate on the limitations of the procedure. In general, global optimization techniques such as Particle swarm optimization do not rely on gradient descent method used by local optimization techniques like the Levengerg-Marquardt method, and hence don't require a differentiable objective function. Such a characteristic helps to determine if the lack of sensitivity of surface motion can be attributed to the chosen optimization algorithm or it is linked to model parameters. Additionally, a global optimization method has the advantage that a large parameter space can be searched from different initial starting points without having prior knowledge on the optimum solution. However, in order to obtain a physically interpretable solution and to reduce the computational time, it is necessary to restrict the search range. We achieve this by assigning values for each of the model parameter from previous experimental characterizations and then extending their ranges by an order of magnitude [8][9][10].
In particular, we choose the sample parameters suitable for PS-LDPE material and generate synthetic data sets based on the interaction with a Silicon cantilever. The sample properties used for the simulations is provided in table S2.2 together with the following cantilever properties: f 0 = 163 kHz, Q = 491, k = 23.95 N m −1 , the effective driving force F d = 1.39 nN and the unperturbed height h = 22.6 nm. Next, we use random sampling to select different starting parameter sets. A total of 15 different parameter sets are created and simulated with the moving surface model to generate the amplitude and phase frequency components which are then used as references (i.e. parameter sets for which the objective function is zero) for the Particle swarm based global optimization. Table S2.2 shows the optimization results for 4 randomly chosen parameter sets out of 15 simulated data sets. The results show that tip-sample dynamics is well approximated with low error values E, but the identified parameter values are far from their true values. This deviation is far more significant for surface parameters in comparison with bulk parameters. Once again, we attribute this issue to non-convexity and lack of sensitivity of surface parameters as discussed in the main manuscript. Additionally, Figs. S2.2 and S2.3 show the time trace of the cantilever and the associated surface motion together with the force quadratures for both the original dynamics coming from the model simulations and the identified dynamics resulting from optimization. In both the figures, while we observe a good agreement for the force quadratures and cantilever motion, the identified motion of the sample surface does not match with the simulated motion (See Figs S2.2(g)-(h) and S2.3(g)-(h) ). This further confirms the trivial contribution of the surface motion on amplitude and phase of intermodulation components.   The probe height h is included in the set of unknown parameters (see main manuscript modeling section III). In general h varies with the working height of the cantilever which in turn depends on how much the feedback control moves the z-piezo during the scanning operation. By taking advantage of the conservative quadrature, in phase with the cantilever motion, it is possible to estimate an approximate value for h based on the onset of repulsive forces.
We suggest two criteria for extracting h from force quadratures as illustrated in Fig. S2.4. We assume the maximum of the in-phase force component (related to adhesion) is achieved closely after the tip starts to penetrate the sample. Thus, the first criterion (denoted by red crosses in Fig. S2.4 (b)) is taken at the middle of the increasing part of F I , whereas the second one corresponds to the amplitude where the in-phase component starts to increase. We browse and apply these two criterion on all pixels of the black line displayed in Fig. S2.4(a). The comparison shown in Fig. S2.4(b) highlights a better match between the heights corresponding to LDPE pixels using the first criteria, when the second criteria seems more suited for the pixels related to PS material. That can be explained by the different material properties, for instance the larger stiffness for PS causes a faster increase of F I , whereas in case of the softer material the short-range adhesive force is more significantly involved before the tip starts to indent the sample. The analysis of these force quadrature curves could be further developed using a more accurate tip-sample force model such as Attard's model [11][12][13][14][15], in order to describe first the transition between the non-contact and adhesive regime, and secondly the transition between the adhesive and repulsive regime.

S3. HIGH VOLUME GRADIENT BASED OPTIMIZATION AND INITIAL POINT SELECTION PROCEDURE
In this section, we discuss the results obtained using the Levenberg-Marquardt algorithm from multiple initial points for both models with and without sample's surface motion. This is done to analyze the sensitivity of the model on initial starting points for the optimization. We begin by creating a numerical range for each parameter based on previous literature studies. Then, a grid of initial starting points is chosen and for each initial point we perform the optimization routine. The distribution of the identified parameters is analysed with histograms and by fitting Gaussian function to extract statistics. The distribution are discussed for each model separately in the following sections.

A. Piecewise linear model with surface motion
Using the moving surface model, we run multiple gradient-based optimizations for pixel (i) and pixel (iii) of Fig. 2 in the main manuscript with the grid of initial parameters defined in table S3.3. The grid includes 3 different values per parameters, chosen in such a way that the parameter exploration recovers a large parameter space (including notably at least one order of magnitude in the case of the viscoelastic properties), and that all routines are performed within a reasonable computational time. In total, 3 6 = 729 optimizations were performed, starting from all the combinations of the grid. In this section we present the histograms used to extract the values reported in table I of the main manuscript.   Here, we report the results and histograms obtained from the large set of optimizations carried out using the piecewise linear model without surface motion. We begin with a set of 3 4 initial parameters defined by the grid presented in table S3.4, and analyze the parameter distributions in the same way as outlined in the previous section.
With the 4 parameters model, statistic for the identified parameters depicts well defined Gaussian distributions that are specific for each type of material. Additionally, the mean of the Gaussian distributions correspond to the lowest values of the objective function. This is shown in Figs. S3.5 and S3.6 for PS and LDPE material sampled at pixel locations (i) and (iii) of Fig. 2 in the main manuscript. The parameter values from the optimization procedure are reported in table S3.5.  Fig. S3.7). The uncertainties are estimated with a 95% confidence interval.
From this statistical analysis, we extract a reduced set of starting parameters. The parameters summarised in table S3.6 have been used to obtain the results showcased in Fig. 5 of the main manuscript. The two first initial points in Table S3.6 were selected by identifying the mean values (also corresponding with the lowest error) among the final results displayed in Figs. S3.7 and S3.8. In addition, we add a third initial point leading to identified parameters within the confidence intervals for all parameters and both the pixels. We detail the final parameters and errors obtained on pixels (i) and (iii) with these three initial points in