Functional dependence of resonant harmonics on nanomechanical parameters in dynamic mode atomic force microscopy

Introduction

When an AFM cantilever oscillating freely and harmonically at a given frequency f and amplitude A₁ approaches a solid surface, the oscillation becomes anharmonic due to the non-linear interaction, represented by the force field F_ts, between the cantilever tip and the surface [1]. Thus, the time dependent trajectory a(t) of the cantilever tip, which can be expressed in the harmonic limit by a(t) = A₁cos(2πft), is transformed into a Fourier series with harmonic oscillations of amplitudes A_n and frequencies f_n = nf [2,3]:

where a is a constant amplitude value. The A_n coefficients can be analytically calculated, in the limit of weak tip–sample interaction (A₁ >> A_n for n > 1), by integrating the force field F_ts that is modulated by weighted Chebyshev polynomials of the first kind, T_n(u) [4,5]. The simplified expression is:

where k_c stands for the cantilever stiffness, z is the distance between the cantilever base and the sample surface and T_n(u) = cos(ncos⁻¹(u)). Note that A_n decreases for increasing order number (n) and k_c values.

The tip–surface interaction F_ts can be expressed as a function of experimental parameters, such as the tip radius (R) and the Young’s modulus of the sample (E). A list of well-accepted models can be found in the literature, including the most widely used Hertz, Derjaguin–Muller–Toporov (DMT) and Johnson–Kendall–Roberts (JKR) models, describing the analytical dependence on such parameters [6-10]. The DMT model, which will be used in this work, has the following expression in the repulsive regime:

where H is the Hamaker constant, a₀ the intermolecular distance, d the tip–sample gap (related to z) and E* the reduced Young modulus, which includes the contribution from the tip. Thus, by combining Equation 2 and Equation 3 it is possible to determine the dependence of A_n as a function of all relevant parameters. However, A_n can be hardly solved analytically, so that numerical simulations are required. On the other hand, the dynamics of the oscillating cantilever cannot be oversimplified and flexural eigenmodes must be considered, in particular when they are located close to one of the harmonics [11].

Here, we report a combined theoretical and experimental study on the functional dependence of the amplitudes of higher harmonics on relevant parameters such as the tip radius, free oscillation amplitude, cantilever stiffness and sample Young’s modulus. Because of the low amplitudes of the involved harmonics (well below 1 nm), we concentrate on the repulsive regime of the tip–sample interaction and on those harmonics close to flexural eigenmodes of rectangular cantilevers, hence the term resonance, so that their intensities can be reliably determined [11,12]. Previous results have shown that the intensity of the 6th harmonic is noticeably larger that the intensities of the neighbouring 5th and 7th harmonics, respectively, using cantilevers with 350 kHz resonant frequency and polystyrene samples [13]. Such study has provided a practical qualitative method to continuously monitor changes in the shape and radius of a cantilever tip in amplitude modulation AFM mode. Here, we will also focus on the 6th harmonic since it provides the larger amplitude, but the methodology can now be extended to other harmonics.

The work presented here provides a tool to help to prepare experiments. Such tool is a mathematical expression describing the evolution of the 6th harmonic on experimental parameters such as tip radius, free amplitude, cantilever and sample’s stiffness. For given experiments directed to, e.g., the determination of the evolution of the tip radius evolution and of sample stiffness, the mathematical tool should help in the selection of the favourable range of parameter values (including amplitude setpoint).

Results and Discussion

Simulations

Let us analyse first the dependence of the amplitude of the 6th harmonic, A₆, on the independent parameters R, E and A₁, respectively, for a fixed k_c value. We will analyse the dependence of A₆ on k_c at a later stage. Note that A₁ can be externally and continuously varied by selecting the excitation amplitude and frequency of the cantilever base and that R and E depend both on the materials used. Figure 1a and Figure 1b show the calculated approach curves (as a function of z) corresponding to the amplitude and phase lag φ of the fundamental mode, respectively, for A₁ = 30 nm, R =10 nm and E = 3 GPa using a 25 N/m silicon cantilever with a resonant frequency f₀ = 300 kHz and a Q factor of 400. The threshold of the repulsive region (φ > 0 degrees) is represented by a vertical discontinuous grey line. Figure 1c displays the A₆ approach curves corresponding to the parameters used in Figure 1a and Figure 1b (continuous black line), as well as those obtained by increasing R (continuous red line), E (continuous blue line) and A₁ (continuous magenta line), respectively. The black line shows the shape of the A₆(z) curve, with A₆ > 0 above the repulsive regime threshold, exhibiting a maximum value at z ≈ A₁/2 and decreasing back to zero for sufficiently small z values, which correspond to small oscillation amplitudes, as depicted from Figure 1a.

From Figure 1c we readily observe that A₆ increases for increasing R, E and A₁ values at a given z value, respectively. When R increases (continuous black and red curves) we observe, in addition, that the maximum of the curve shifts towards lower z values, a behaviour that is also observed for increasing E values (compare continuous black and blue lines). In fact the weight of the area under the curves is shifted towards lower z values in both cases.

The comparison between the continuous black and red curves indicates that the variation of the radius has an impact in the value of A₆ only for certain values of the amplitude setpoint (A_sp): for example at A_sp = 25 nm (z = 24.2 nm) the difference between the black and red curves is clearly smaller than for A_sp = 15 nm (z = 14.2 nm). The same argument applies for the case of increasing E values. It means that the previous modelling of the system is necessary in order to select the appropriate experimental conditions.

Let us now explore in detail the dependence of A₆ as a function of R, E and A₁. The corresponding curves are shown in Figure 2, where the filled black points correspond to calculated values and the continuous lines to least-square fits. The simulations have been performed for silicon with the following parameters: k_c = 26 N/m, f₀ = 300 kHz and Q = 400, with A_sp = 0.5A₁. All curves confirm the increase of A₆ for increasing values of R, E and A₁ but with different evolutions. In the case of the phase of the 6th harmonic the behaviour is opposite, i.e., it decreases for increasing R values (see Supporting Information File 1, Figure S1).

In all cases the points closely follow a gun-shaped function with the expression:

where g₁, g₂ and g₃ are real numbers. This function covers the limiting cases of G(x) = x/g₁ for g₁ ≠ 0 and g₂ = g₃ = 0, G(x) = 1/g₂ for g₂ ≠ 0 and g₁ = g₃ = 0 and G(x) = √x/g₃ for g₃ ≠ 0 and g₁ = g₂ = 0. Figure 2c exemplifies the case where the dependence is nearly linear, an approximation certainly valid for smaller intervals around a central amplitude value. In the case of Figure 2a and Figure 2b we clearly observe two well defined regions corresponding to lower and higher R and E values, respectively. At lower values A₆ exhibits a strong variation as a function of R and E (large slope) while at higher values A₆ tends towards an asymptotic limit (small slope). Thus, for the k_c values used in this simulation (26 N/m) A₆ is most sensitive for values of tip radius below ca.10 nm (ultrasharp tips) and sample’s Young’s modulus below ca. 20 GPa (typically polymers).

The dependence of A₆ on the cantilever properties is more complex because the cantilever shape, the elastic constant k_c, the resonance frequency f₀ and the quality factor Q are interrelated parameters that all contribute to the dynamics of the oscillation. For simplicity, we characterize the cantilever by k_c. Sader’s formula, also termed Sader’s method, gives a well-accepted expression of the dependence of k_c on the other parameters [14]:

where ρ_f is the density of the fluid (in our case air), S_c is the plan view area of the cantilever (width × length) and Γ_i the imaginary component of the hydrodynamic function [15]. This expression is valid for Q >> 1.

In order to obtain an approximate manageable expression of the functional dependence of A₆ on k_c we have used the nominal geometrical and resonance frequency values of different commercial rectangular cantilevers, as provided by the manufacturers and Q factor values in the range 200 ≤ Q ≤ 600. The resulting k_c values obtained from Equation 4 have been used as input to the VEDA code [16]. The result from this calculation is shown in Figure 3. The selected cantilevers with their values are shown in the caption to the figure.

The calculated points in the figure also follow quite closely Equation 4, with g₁ = −63.092 nN·nm⁻², g₂ = 42.371 nm⁻¹ and g₃ = −42.223 nN^1/2·nm^−3/2.

Thus, the simulations provide a functional dependence of A₆ on the relevant R, E, A₁ and k_c parameters, which can be stated as:

where G_r, G_e, G_a and G_k represent the decoupled gun-shaped functions with their corresponding coefficients g. The apparently different behaviours can be described by a single universal curve, where the magnitudes and signs of the g coefficients determine the final shapes. This expression can substantially simplify the analysis of partial contributions. In addition, it paves the way to define a methodology for finding the optimal experimental conditions to monitor or characterize a specific magnitude from the acquisition of the amplitude of a higher harmonic.

Experiments

In the following section we compare the proposed functional dependence with experimentally derived results to validate the modelling and simulations described in previous sections.

Dependence of the amplitude of the 6th harmonic on tip radius

Let us first explore the shape of the experimental approach curves and compare to the calculations shown in Figure 1. Figure 4 shows simultaneously acquired experimental approach curves, i.e., amplitude (a) and phase (b) of the fundamental mode and amplitude of the 6th harmonic (c), respectively, as a function of the piezo displacement in the direction perpendicular to the surface (z) for a nominally k_c ≈ 44 N/m cantilever with resonance frequency 293 kHz on silicon surfaces and A₁ = 34 nm. The black curve in Figure 4c was taken at the beginning of the experiment (fresh tip). The rest of the curves are acquired after the acquisition of several intermediate images (i.e., several images are taken in between each approach curve). The order in which the curves are taken was black, green, red and blue, respectively.

The first conclusion that we can extract is that the simulated curves (Figure 1a and Figure 1b) reproduce well the shape of the experimental ones (Figure 4a and Figure 4b). In addition, Figure 4a and Figure 4b provide evidence that both the amplitude and phase of the fundamental mode remain unchanged during the experiments, except for a small variation in the phase at the transition between attractive to repulsive regimes (φ ≈ 0). Conversely, the value of A₆ as a function of the piezoscanner displacement exhibits an increase over the noise level above the threshold corresponding to the onset of the repulsive regime. We clearly observe an increase of the maximum value of A₆ as well as an increase of the position of the maximum, a trend that is reproduced by the simulations, as depicted from Figure 1c. Assuming that both sample’s Young’s modulus and k_c remain constant, then the increase of A₆ can be attributed primarily to an increase in R.

The small variation of the phase (see Figure 4b) occurs at the transition between attractive to repulsive regimes, where the critical amplitude A_c is defined [17]. It turns out that A_c depends on R, closely following a power law function (A_c R^m), where the parameter m (m < 1) depends on the particular cantilever used [17-20]. This provides a further evidence that the observed increase of A₆ can be attributed to an increase in R. On the other hand, A_c is evaluated, in those works, at the sharp attractive–repulsive transition, which implies a rather involved experimental determination, while using approach curves one can select the setpoint and thus the A₆ intensity in a larger range (within the repulsive mode). However, larger repulsions may lead to wear, and thus to damage of the tip [21-25].

Additional information can be obtained from the acquired images. Figure 5 shows the evolution of the mean A₆ value acquired simultaneously with the topographic and phase images. As it is observed from the simulations and experiments, the value of A₆ is below 0.1 nm, which implies a very low signal to noise ratio. To overcome this difficulty, A₆ is acquired at each point of the AFM image (256 × 256 points) and then it is calculated by averaging all the values obtained at the image. In this way, the time evolution of the value of A₆ is expressed in terms of sequentially acquired images (each point corresponds to one image), where the experimental parameters such as A₁ and the amplitude setpoint do not change over time. We observe a rapid increase from image 1 to 10 followed by an increase with a smaller slope above image number 10. The figure resembles Figure 2a with a higher slope at the beginning and a lower slope afterwards. Because of the expected tip wear, the evolution observed in Figure 5 can be again ascribed to an increase in tip radius. This method has been proposed to monitor the stability of the tip in a continuous manner [13,26,27].

The quantitative determination of the actual tip radius at each image is a rather difficult task, since it depends critically in several experimental parameters. We have performed a parallel calibration of the tip radius with reference samples. Supporting Information File 1, Figure S2 shows the evolution of A₆ as a function of R obtained from commercial gold nanoparticles (5.5 ± 0.7 nm diameters) dispersed on a thin poly-lysine film grown on mica. From the figure we can observe the increase of A₆ for increasing R values.

Dependence of the amplitude of the 6th harmonic on bulk modulus

Figure 6 shows the evolution of A₆ for discrete values of Young’s modulus from different materials, namely PDMS (E = 0.0025 GPa), LDPE (E = 0.1 GPa), PS (E = 2.7 GPa), fused silica (E = 72.9 GPa), titanium (E = 110 GPa) and sapphire (E = 345 GPa), using a 10.9 N/m cantilever, as determined with the thermal tune method and Sader’s corrections [28].

In spite of the reduced number of experimental points, the curve can be compared to Figure 2b, with a sharp increase at low E values and a nearly constant value for larger E values. One has to take into account that during the experiments, where the tip has to be retracted and samples have to be changed, a variation of the tip radius cannot be excluded. The sample with the highest wear was titanium, because of its higher roughness as compared to the rest of the samples, and for this reason it was measured at the end of the cycles. From the figure it can be clearly observed that A₆ for the titanium sample shows the larger values.

Dependence of the amplitude of the 6th harmonic on free oscillation amplitude

Figure 7a shows the evolution of A₆ as a function of the free oscillation amplitude as determined experimentally with nominally 26 N/m cantilevers. The experimental points (full black circles) correspond to the mean A₆ values obtained from the average of a full image acquisition of the 6th harmonic, in analogy to the experiments described in Figure 5. The points have been fitted to the G function from Equation 4. The nearly linear behaviour agrees with the predictions shown in Figure 2c and strongly suggests that the tip radius has not changed during the measurements, since otherwise the behaviour of A₆ would have been nonlinear, with a slope increasing for increasing A₁ values.

In addition, Figure 7b shows the evolution of A₆ as a function of sequentially acquired images using softer cantilevers (k_c ≈ 3 N/m), where both A₁ and A_sp have been intentionally varied during the experiment. Note that here A₆ increases for increasing A₁ values, as observed in Figure 7a, and that for a particular couple of A₁ and A_sp values, A₆ increases due to tip wear, and thus to an increase in tip radius.

Note the higher A₆ values in Figure 7b, as compared to those in Figure 7a for similar A₁ values. This is essentially due to the softer cantilevers used in Figure 7b, an effect described in Figure 3.

Conclusion

Based on numerical simulations using the VEDA code we have proposed a functional dependence of the amplitude of the 6th harmonic of rectangular cantilevers on the tip radius, sample’s Young’s modulus, free oscillation amplitude and cantilever force constant. The 6th harmonic has been chosen because it is in close resonance with the 2nd flexural mode, leading to observable signals, and because its frequency is within the reach of the control electronics used in the experiments. The outcome of the simulations is that the 6th harmonic can be analytically expressed by the product of four identical decoupled gun-shaped functions, each one associated with a specific parameter and with its own coefficients. Thus, the partial evolution for a particular parameter can be traced using this universal behaviour.

The simulations have been validated with AFM experiments using rectangular cantilevers in the 3–45 N/m range and different samples in the 2–345 GPa range. The predicted trends are well reproduced by the experiments. As we can notice from the different trends, the 6th harmonic is most sensitive to changes in tip radius for values of tip radius below ca. 10 nm and Young modulus below ca. 20 GPa. If we consider the method for implementing tip radius real time monitoring, this will be more effective when sharp (new) tips are imaging stiff samples. Instead, in order to measure the surface Young modulus the modelling shows that the best results will be obtained using larger tip radius (10–15 nm) in order to have an almost constant radius-dependent contribution.

So far the results are only qualitative. In order to obtain trustworthy predictions, a more precise and more accurate calibration of the cantilever is still necessary. The performed simulations might also be improved including parameters which could be varied and measured experimentally in order to refine the method and make the results of the simulations more comparable to the experimental ones. The mathematical tool that we present here based on gun-shaped functions allows to identify the optimal conditions to obtain information about the properties of materials from the harmonic response in non-contact dynamic modes.

Simulations and Experimental Details

Simulations

Simulations have been performed using the Virtual Environment for Dynamic AFM (VEDA) open code, which takes into account the dynamics of oscillating rectangular cantilevers with multiple eigenmodes [16]. The frequencies, stiffness and quality factors of the 2nd flexural eigenmode have been approximated by the known relationships 6.27f, 6.27²k_c and 6.27Q, respectively, corresponding to a massless tip [29]. As mentioned above, we have used the DMT model to describe the tip–surface interaction and the tip has been approximated by a hemisphere with a well-defined radius R. For simplicity neither viscoelastic nor capillary forces have been considered. The negligible contribution of viscoelastic forces in PS with 25 N/m cantilevers is discussed in the Supporting Information File 1 (see Figure S3).

Experimental

AFM experiments have been performed with a Bruker DIMENSION ICON instrument hosted in a homemade controlled humidity environment with Nanoscope V control electronics. Commercial rectangular microfabricated silicon cantilevers with ultrasharp silicon tips (R < 10 nm) have been used with the following nominal values: OTESPA (Bruker) with k_c ≈ 44 N/m and 300 ≤ f₀ ≤ 400 kHz, OTESPA-R3 (Bruker) with k_c ≈ 25 N/m and 200 ≤ f₀ ≤ 400 kHz, and PPP-FMR (Nanosensors) with k_c ≈ 3 N/m and 70 ≤ f₀ ≤ 80 kHz, where f₀ stands for the resonance frequency. The amplitudes of the higher harmonics were registered using an internal lock-in amplifier. In general, such amplitudes will depend on the selected working frequency. The effect of working slightly off-resonance is discussed in Supporting Information File 1, Figure S4. The estimation of the A_n magnitudes (in nm) has been obtained by calibration of the laser-detector sensitivity, which is about 100 nm/V, as determined from force curves. Humidity was kept below 20%. The z motion of the piezoscanner has been calibrated using virtual standards [30]. Due to the value of the fundamental resonance frequency of the employed cantilevers, we have focused in the resonance of the 2nd flexural mode and the 6th harmonic, which frequencies are below 2.5 MHz, since the control electronics is limited to 5 MHz.