Robustness of attractors in tapping mode atomic force microscopy

In this work, we perform a comprehensive analysis of the robustness of attractors in tapping mode atomic force microscopy. The numerical model is based on cantilever dynamics driven in the Lennard–Jones potential. Pseudo-arc-length continuation and basins of attraction are utilized to obtain the frequency response and dynamical integrity of the attractors. The global bifurcation and response scenario maps for the system are developed by incorporating several local bifurcation loci in the excitation parameter space. Moreover, the map delineates various escape thresholds for different attractors present in the system. Our work unveils the properties of the cantilever oscillation in proximity to the sample surface, which is governed by the so-called in-contact attractor. The robustness of this attractor against operating parameters is quantified by means of integrity profiles. Our work provides a unique view into global dynamics in tapping mode atomic force microscopy and helps establishing an extended topological view of the system.


Introduction
Atomic force microscopy (AFM) has become one of the most prominent characterization tools in modern science. It is ubiquitously used to characterize and manipulate surface properties of materials down to atomic resolution in both air and liquid environments [1]. Among the various AFM operational modes, tapping mode AFM (TM-AFM) [2] also known as amplitude modulation AFM (AM-AFM) [3] is one of the most extensively used techniques to obtain high-resolution images of wide variety of samples. Owing to reduced lateral forces and high-phase sensitivity, TM-AFM is widely popular in soft matter [4], biological samples [5] and polymer [6] applications.
TM-AFM is based on the near-resonant excitation of a microcantilever with a sharp tip at its free end that is vibrating in the vicinity of the sample. The vibrations in the cantilever are influenced by the tip-sample interaction forces, which modifies the beam dynamics. In general, the tip-sample interaction is nonlinear and comprises long-range attractive Van der Waals forces, shortrange quantum mechanical repulsive forces, adhesive and contact forces. The nonlinear response due to tipsample interaction is even more involved in the presence of electrostatic and capillary forces [7]. In the existence of such complex nonlinearities, a comprehensive understanding of the multi-stable response is crucial, since these nonlinearities can be efficiently utilized for extracting several nanomechanical properties of the sample [3,[8][9][10]. Furthermore, the accuracy of imaging and nanomechanical characterization of sample surfaces from measured data depends crucially on the deconvolution of data with appropriate models.
In recent years, the underlying dynamics of AFM cantilever and its exploitation have been investigated by many authors [11][12][13][14][15][16][17][18]. The vast majority of these studies have been dedicated to study the microcantilever dynamics by utilizing a simplified single-degree-offreedom point mass model [16][17][18]. In this approach, first a static analysis is performed to determine the cantilever stiffness and then the equivalent mass is calculated based on the experimentally evaluated fundamental resonance frequency. Furthermore, the excitation is modelled as an external force acting on the point mass [16][17][18]. Such a lumped parameter model does not represent the conditions encountered in reality, where the microcantilever is subjected to base excitation at the clamped end by means of a piezoelectric actuator. The base excitation induces linear and nonlinear parametric excitations that are typically not captured by the lumped parameter model, thus failing to accurately describe the dynamics of the microcantilever [13,14].
Continuous beam models on the contrary have proven to predict the nonlinear aspects of AFM cantilever dynamics accurately [12][13][14][15]. These models are able to provide precise and deeper insights into the physics behind the nonlinear phenomena such as amplitude jumps, period-doubling, and grazing bifurcations [14,19]. Therefore, in spite of the complexities involved in modelling of the AFM cantilever as a continuous beam, it allows for capturing the overall nonlinear aspects of the AFM dynamics [20].
In practical operation, the dynamics of AFM cantilever is influenced by several operating parameters such as tip radius, excitation amplitude, excitation frequency, and feedback values. The real-time variation of these parameters during an AFM operation can lead to unwanted dynamical phenomena such as bifurcations and unstable and aperiodic motions, which can decrease the reliability of results and strictly limit the operating ranges of the AFM. The dynamic models as well as the feedback strategies implemented in standard AFM systems tend to focus on governing the local dynamics of the system, leaving the impact on the global dynamics largely unknown. Therefore, in order to predict and control these dynamical events as well as preserving the stability of operation in TM-AFM, it is important to study and understand the nonlinear responses from a global perspective.
Currently, there are no detailed works on global dynamics of TM-AFM, which (i) evaluate the escape boundaries, (ii) estimate dynamical integrity, and (iii) perform detailed analysis of bifurcations. Existing literature has focused on the dynamical integrity and bifurcation scenarios of non-contact AFM [19,21]. But basins of attraction and erosion process of basin portraits in TM-AFM are lacking in the literature. The erosion of uncorrupted basins of attraction surrounding each main solution as a function of AFM operation parameters is of paramount importance from both theoretical and experimental perspective.
In this article, we elucidate the global dynamics and robustness of attractors in TM-AFM. Differently from existing works, which are based on limited analysis on the local dynamical behaviour, this work makes systematic use of bifurcation diagrams to highlight the appearance and disappearance of steady-state solutions. The latter offers an overall interpretation of the dynamic response with respect to operational parameters, namely excitation frequency and forcing amplitude. Additionally, we show that, by changing the operation parameters, a microcantilever initially in the primary resonant branch can escape from its local potential well and get captured by a second coexisting local potential well close to the sample surface. This regime of oscillation close to the sample surface, referred in this article as 'in-contact' attractor, is unexplored in the literature. In this work, we examine its evolution and robustness properties with frequency response curves and basins of attraction. Furthermore, the envelopes of local bifurcation boundaries are built to understand the escape scenarios of the solutions. Finally, integrity analysis is performed to quantify the steady-state solutions associated with all the attractors present in TM-AFM.
Based on these motivations, the article is organized as follows: the modelling of the system is discussed briefly in Sect. 2. This is followed by a detailed analysis of frequency response of the system, bifurcation charts and response scenarios including escape of solutions in an excitation parameter space in Sect. 3. Particular attention is paid to the frequency response analysis of the in-contact attractor around its primary and parametric resonance. The results of this section are then utilized to build basins of attraction and dynamical integrity curves for the main attractors present in the system in Sect. 4.

Numerical model
The classical beam theory, based on the Euler-Bernoulli assumptions, is used to obtain the continuous model for the AFM microcantilever shown in Fig. 1. The nomenclature used to describe the equations in this article is identical to the one described by Ruetzel et al. [12]. The deflection of the cantilever towards the sample is treated as positive, and the rest position of the cantilever is taken as reference. The considered microcantilever has a length L, mass density ρ, Young's modulus E, area moment of inertia I, and cross-sectional area A. The beam is clamped at x = 0 and free at x = L. The tip-sample separation distance in the reference configuration is denoted by Z , and the total deflection of the microcantilever w(x, t) can be expressed as w( is the deflection of the microcantilever relative to a noninertial reference frame attached to the base and w * (x) is the static deflection towards the sample due to tipsample interaction. The base excitation is generated by a dither piezo and is assumed to be harmonic, i.e. y(t) = Y sin(Ωt), where Y and Ω are the amplitude and frequency of excitation, respectively.

Tip-sample interaction
In TM-AFM, the microcantilever oscillates in close proximity to the sample surface. In our work, we use the Lennard-Jones (LJ) potential to describe the tipsample interactions [22]. Although the model doesn't take into account the real contact mechanics encountered in TM-AFM, it represents a generic tip-surface interaction potential which mimics qualitatively, the more detailed and computationally expensive models [13]. The LJ potential models the non-retarded dispersive Van der Waals forces as well as the short-range repulsive exchange interactions between two molecules. Assuming a spherical tip apex with radius R and a flat sample surface, the interaction potential and the force are: where z is the instantaneous tip-sample separation gap. A 1 and A 2 are the Hamaker constants for the repulsive and attractive potentials, respectively. A positive interaction force implies repulsion. The Hamaker constants are A 1 = π 2 ρ 1 ρ 2 c 1 and A 2 = π 2 ρ 1 ρ 2 c 2 , where ρ 1 and ρ 2 are the number densities of molecules in the interacting media and c 1 and c 2 are the interaction coefficients of intermolecular pair potential [22].

Equation of motion
The nonlinear static deflection of the microcantilever in the absence of base excitation is computed by solving for the equilibrium gap between the tip and the sample shown in Fig. 1. The static equilibrium gap η * at the free end is calculated as a function of the approach distance Z through static balancing of the cantilever restoring force and the tip-sample interaction forces. The dynamic equation of motion of the tip deflection u(x, t) about its nonlinear equilibrium subjected to base harmonic motion is then derived through a single mode discretization of the Euler-Bernoulli beam equation. The interaction forces given by Eq. (1b) are assumed to be acting on the free end of the cantilever and are mathematically achieved through the use of Kronecker delta (δ) function. Writing the equation of motion of the vibrating cantilever with respect to a (a) (b) Fig. 1 Schematic of the AFM cantilever in a static deflection configuration and b configuration at which AFM cantilever exhibits oscillations around the nonlinear static equilibrium. The sample is atη 1 = 1 non-inertial frame of reference leads to the following governing equation, Equation (2) is a non-autonomous and nonlinear equation. The equation is discretized by projecting the dynamics onto the system's linear modes of vibration. The natural modes and frequencies are obtained using the Galerkin approach [12,23]. The frequency range of the analyses in this paper spans around the neighbourhood of the fundamental resonance, where the contribution of higher modes is substantially negligible. Based on this assumption, a single-degree-of-freedom model is used. Assuming u(x, t) = φ 1 (x)q 1 (t) (where φ 1 is the first approximate eigenfunction around the static deflected configuration) and using the Galerkin approach, the following nonlinear equation can be derived [12]: The dimensionless variables and the corresponding coefficients are described in "Appendix". The microcantilever tip deflection towards the sample is denoted byη 1 . In addition, the equation is made dimensionless with respect to equilibrium gap width (η * ) and the fundamental frequency of the free microcantilever (ω 1 ) in the absence of tip-sample interaction forces. The amplitude of the dither piezoelectric actuator is denoted byȳ. The dotted quantities represent derivatives with respect to rescaled time τ (τ = ω 1 t). Finally, the modal damping d 1 is explicitly introduced in Eq. (4) and is related to the quality factor Q of the cantilever by the relation Q = 1/d 1 .

Numerical analysis
In order to investigate the dynamical behaviour of the TM-AFM, the simulations of the model given by Eq. (4) are performed in this section. The entire analysis is carried out for the interaction of a soft monocrystalline silicon microcantilever with the (111) face of flat silicon sample. The cantilever and interaction properties are listed in Table 1. Furthermore, the analysis is performed for tip-sample gap (η * ) in the bi-stable region of nonlinear elastostatic equilibrium curve. Thus, a value of η * = 6.542 nm 1 is chosen for the rest of the analysis in this article. However, a similar analysis can be carried out for any other tip-sample gap values. Numerical simulations are performed by using a pseudo-arc-length continuation technique [24]. We also make use of basins of attraction (phase space) in order to illustrate the presence of various attractors (steady-state solutions). A basin of attraction is a set of possible initial conditions about an equilibrium point in phase space that assures a specific response from the cantilever. In other words, any chosen initial condition within the phase space will be 'attracted' to a particular steady-state motion of the cantilever. Unless specified, the basins are evaluated in a phase-space grid ofη 1 = [−0.9, 0.9] andη 1 = [−0.9, 0.9] as it contains all the main attractors involved in the system potential well.  Table 4. (Color figure online)

Frequency response and bifurcation scenarios
The interaction between tip and sample gives rise to different nonlinear frequency responses depending on the tip-sample separation distance. The response displayed in Fig. 2a shows an initial softening behaviour when the tip is far away from the sample. This region is dominated by attractive Van der Waals forces. However, when the tip-sample separation reaches the order of the interatomic distance, the response exhibits hardening behaviour and this region is dominated by repulsive forces (see Fig. 2a). Figure 2b shows the corresponding basin portrait associated with the frequency response atΩ = 0.8. The blue and crimson basins together form the attractive region, while the purple basin belongs to the repulsive region.
The nonlinear frequency response shown in Fig. 2a is well known and studied extensively by many authors [10,12,25]. However, there exists another overlooked steady-state response when the cantilever escapes the local potential well and gets trapped by a subsequent attractor very close to the sample. Figure 3a shows the frequency response for the model given by Eq. (4), and Fig. 3b shows the position of the in-contact attractor in the basin portrait. Note that Fig. 3a is made up of two different solutions belonging to the attractive (lower frequency response curve) and the in-contact attractor (upper frequency response curve). The two solutions can be obtained by using different initial conditions in the numerical integration. In the next section, the detailed analysis of the in-contact attractor's frequency response and corresponding bifurcation scenarios is presented parametrically.

Frequency response curves
Local dynamic analysis is performed using frequency response curves together with the bifurcation charts. The analysis offers a complete overview into the bifurcations and escape scenarios of the system. Furthermore, the unstable solution branches shown in the frequency response curves are not discussed in detail but are reported for the sake of completeness. tively. In both figures, red and blue circles indicate period-doubling and saddle-node bifurcation points, respectively. Although the dimensionless natural frequency of the system oscillating in attractive regime is found to beΩ 0 = 0.83 (Ω 0 is the natural frequency affected by the system potential well), it can be observed in Fig. 4 that the first natural frequency of the system oscillating in the in-contact regime is Ω 0 = 3.72 with corresponding parametric resonance at 2Ω 0 = 7.44. This shift in resonance frequency is due to the presence of strong repulsive forces which act as a hard spring connecting the cantilever to the sample. This can be visualized as a change in the boundary conditions of the cantilever similar to that of a clamped-clamped beam. Interestingly, the system also exhibits softening nonlinearity in spite of the presence of repulsive forces. This is due to the fact that the cantilever is oscillating in the potential well with a duration of oscillation longer in the attractive regime. In addition, multi-stability can be observed with different solutions overlapping in several discrete ranges of frequencies (see Fig. 4a). At lower values of forcing frequencyΩ ≤ 1.85, the system has only one non-resonant low-amplitude solution (LP 1 ). The corresponding basin associated with LP 1 solution at Ω = 0.5 is shown in light brown colour in Fig. 4b.
The figure illustrates the low-amplitude attractor being the dominant solution in the in-contact regime at lowfrequency values. The LP 1 solution eventually gives rise to a superharmonic branch (Su 1 HP 1 ) atΩ = 1.85 via saddle-node bifurcation, and later a resonant highamplitude solution (HP 1 ) atΩ = 3.3. Figure 4c reports the orange basin belonging to HP 1 solution arising from the boundaries of Su 1 HP 1 solution (green basin) in the in-contact regime. Furthermore, the HP 1 branch in Fig. 4a destabilizes with the inception of a pair of period-2 branches via flip/period-doubling bifurcations. One of the perioddoubling bifurcation occurs close to the low-frequency saddle-node bifurcation atΩ = 0.3, while the other period-doubling bifurcations occur atΩ = 1.56 (zoomed part of Fig. 4a). This behaviour is similar to the nonlinear cantilever response seen in attractive regime as illustrated earlier in Fig. 3a. Moreover, in Fig. 4a the period-2 branch continuation shows stable motion over a short frequency range before undergoing further period-doubling bifurcation cascade. The subharmonic response associated with the period-doubling bifurcation can be observed around the principal parametric resonance frequency of 2Ω 0 = 7.44. The stable large amplitude period-2 solution (referred to as Pa-P 2 in Fig. 4a), arising from one of the period-doubling points, is found to be stable over a wide range of excitation frequency (Ω ∈ [4.8, 7.7]). In addition to the above analysis, referring to Fig. 5, at larger excitation amplitudes (ȳ ≥ 0.0118), the softening behaviour in the nonlinear response of the system increases along with a larger field of existence of the superharmonic response (Su 1 HP 1 ). Furthermore, a second superharmonic branch (Su 2 HP 1 ) bifurcates through a saddlenode atΩ = 1.79. Here, four forced period-1 solutions coexist out of which only two are stable (Su 2 HP 1 and Su 1 HP 1 ) and two are unstable. The superharmonic branch for larger excitations eventually joins the main branch of the in-contact response via the saddle-node at Ω = 2.83. In addition, the Su 1 HP 1 branch is destabilized over narrow frequency ranges,Ω ∈ [0.92, 0.95] andΩ ∈ [2.41, 2.42], by the occurrence of perioddoubling bifurcations.
It is observed in Fig. 5 that the low-frequency perioddoubling point (Ω = 0.94) has a stable period-2 solution over a short frequency range (see zoomed part of the figure) and the high-frequency period-doubling (Ω = 2.42) presents mostly an unstable bifurcated period-2 solution except for few initial continuation points that are stable. These few stable points are not shown in the frequency response curves but are described in the bifurcation maps. In an analogous way, a period-doubling bifurcation atΩ = 2.90 is observed on the HP 1 branch having a stable period-2 solution limited in frequency rangeΩ ∈ [2.7, 2.9]. Moreover, period-doubling cascades are present in both stable period-2 solutions arising from Su 1 HP 1 and HP 1 branches. This period-doubling cascade can lead to chaos in a similar fashion as encountered in standard TM-AFM systems [26].   Tables 2  and 3, respectively the bifurcation point (saddle-node or period-doubling) with respect to operating parameters, namelyΩ and y. Furthermore, by assembling all the local bifurcation envelopes together, the global response and behaviour map of the entire system in the excitation amplitude and frequency control space are obtained. The bifurcation map has been obtained numerically for the incontact attractor over a wide range of frequencies which includes the fundamental (Ω 0 ) and principal parametric resonances (2Ω 0 ). Thus, it is convenient to analyse the global dynamics by dividing the bifurcation map into two separate regions: the first region focuses around the fundamental resonance frequencyΩ 0 = 3.72 illustrated in Fig. 6, whereas the second region analyses the principal parametric resonance frequencyΩ 0 = 7.44 as shown in Fig. 7. In addition, Table 2 outlines the data concerning the various bifurcation envelopes of Figs. 6 and 7. Furthermore, Table 3 summarizes the various dynamic regions formed by these envelopes and the corresponding solutions involved. From an experimental perspective, the response scenario map provides qualitative information on the form of cantilever response expected for the chosen set of excitation amplitude and excitation frequency.

Analysis of the bifurcation map around the fundamental resonance
Around the fundamental resonance frequency of the in-contact attractor, the dynamics is more involved than in the case of the attractive region [19]. The period-doubling/flip bifurcations appear not only on the main branch of the solution (HP 1 branch in Fig. 5), but also on the first superharmonic branch (Su 1 HP 1 branch in Fig. 5). This drastically increases the possibility of global escape through crisis and also chaotic behaviour through period-doubling cascade. In an analogous way, the saddle-node bifurcations that arise from the superharmonic branch increase the complexity of the response. Figure 6 reports all the possible regions of motions for the cantilever in the range of frequencies surrounding the fundamental resonance. The different regions are named with pink labels and accordingly numbered. For both low-amplitude excitations and frequencies up toȳ = 0.004 andΩ < 1.85, there exists only the LP 1 motion indicated by region R1 in Fig. 6. Hereafter, with the increase inΩ ≥ 1.85, the response consists of both the HP 1 and superharmonic high-amplitude (Su 1 HP 1 ) solutions bound by loci HP 1 -SPF as shown in the region R2 of Fig. 6. Moving to even larger values of excitation frequency, the solution HP 1 governs the behaviour of the system (R9). In particular, special attention should be given to the regions bounded by SPF (red) and SPF 2 (green) loci. These regions exist over a narrow frequency range and invade into other period-1 regions resulting in new multi-stable regions which consist of both period-1 and period-2 solutions. Similarly, in addition to the discussion of aforementioned regions, the presence of other response scenarios and dynamic  A-Pa-P 2 and Pa-P 2 regions data is tabulated in Table 3. Moreover, in contrast to the various stable responses seen in Fig. 6, the bifurcation analysis of in-contact attractor shows the presence of multiple strange attractors leading to crisis scenarios and global escape. The operating parameters leading to escape are depicted by grey regions (R11) in Fig. 6. The crisis scenario also highlights appearance of several rare attractors (period-5 and above) along with period-doubling cascades, which lead the solution to escape from the local potential well.
The detailed instructions on reading the bifurcation chart together with an example are provided in "Appendix B".

Analysis of the bifurcation map around the
principal parametric resonance Figure 7 shows that, below the parametric resonance frequency,Ω < 7.44 , only HP 1 solution exists (R9) until the Pa-P 2 solution from the parametric resonance (2Ω 0 = 7.44) overlaps with the region R9 resulting in two coexisting period-1 (HP 1 ) and period-2 (Pa-P 2 ) responses as shown in region R9 + R12. In between the classic parametric instability tongue (V-shaped region nearΩ = 7.44) formed by the parametric subcritical (Pa-SBF) and supercritical (Pa-SPF) bifurcations, the HP 1 solution becomes unstable and there exists only period-2 (Pa-P2) solution (R12). As expected in any of the archetypal parametric oscillators, the system requires a minimum excitation threshold to achieve parametric resonance. This is indicated by the lifting of the instability tongue (V-shaped region) along theȳ axis and in our system the critical threshold is atȳ = 0.00014. Finally, for excitation frequencȳ Ω > 7.44, only HP 1 solution is present as shown in region R9.

Dynamical integrity and robustness of attractors
In the previous section, we provided insight into diverse solutions and bifurcation scenarios including escape thresholds. However, the analysis did not furnish details on the various instability paths and eventual escape of the steady-state solutions. The information on the instability path (escape from local potential well, cross-well chaos) taken by the cantilever response is of utmost importance in practical applications of AFM. This helps to disentangle image artefacts from factual data.
From an experimental perspective, if the system perturbations can be quantized, then basins of attraction provide insight into the evolution of various steadystate responses (system attractors) and instabilities (erosion profiles) occurring in the system. Furthermore, measures of the basin portraits, the dynamical integrity of the system, are able to quantify the robustness of different attractors. This can be realized through various scalar integrity measures [27]. These integrity measures provide information on the strength of such quantized perturbations required to destabilize the corresponding system response. Therefore, basin portraits together with integrity measures provide a means to track the basin erosion process with respect to changes in operating parameters. Hence, in order to advance the dynamical analysis of the AFM cantilever in the in-contact regime of oscillation, this section focuses on the global topology analysis by means of basins of attraction [28,29].
There have been multiple integrity measures introduced in the literature [30], and this section makes use of two integrity indicators to measure the evolution of phase-space topology, namely local integrity measures (LIM) [27] and integrity factor [31] (IF). The LIM is defined as the normalized radius of the largest hypersphere (circle in 2D), centred on the safe attractor and entirely belonging to the safe basin. It is used to analyse the robustness of the attractor of interest against perturbations. On the other hand, the IF is defined as the normalized radius of the largest circle entirely belonging to the compact part of safe basin. The IF is suitable to study the dynamical integrity of the attractors subject to perturbation around its initial equilibrium condition. The reason to choose these measures with respect to others such as global integrity measures (GIM) relies on the fact that IF and LIM can disentangle the fractality of basin since they focus only on its compact part [27]. In our case this is a serious advantage since the homoclinic tangling of the saddle results in fractalization of the low-amplitude attractive regime basin.

Basin portraits and evolution as a function of tip radius
One of the common causes for image artefacts during AFM scanning operation is the degradation of tip radius R due to repetitive impacts with the sample surface [32]. Furthermore, the correct and reliable operation of the AFM is dependent on the status of the probe tip, since it is responsible for resolving the topography of the sample [10]. The change in radius value during the aforementioned scanning operation thus causes a sensible variation in the system response and corresponding global integrity values. Therefore, in this section, we make an effort to elucidate the variation of system integrity via basin portraits by considering the AFM tip radius (R) as the corresponding varying parameter. Figure 8 outlines the variation of dynamical integrity of system attractors as a function of AFM cantilever tip radius R. The numerical simulations are performed for tip radius ranging from 105 to 225 nm in steps of 3 nm, with excitation frequencyΩ = 0.85 being close to resonance, and excitation amplitudeȳ = 0.005.  Table 4. g Frequency response of microcantilever with radius R = 105 nm and R = 182 nm oscillating with initial condition in attractive region. h Frequency response of microcantilever with radius R = 136 nm and R = 182 nm oscillating with initial condition in the in-contact region. (Color figure online) Fig. 9 Phase-space topology evolution with respect to tip radius (R) from 105 to 225 nm, at fixed excitation amplitude and frequency ofȳ = 0.020,Ω = 0.8, respectively uous and dotted lines in Fig. 8a-c belong to LIM and IF measures, respectively. It is observed in Fig. 8a that in the attractive regime around the radius value R = 114 nm, there is a sharp increase in the high-amplitude resonant solution (shown in Fig. 8g as A-HP 1 ) and the integrity of lowamplitude non-resonant solution (shown in Fig. 8g as A-LP 1 ) rapidly decreases to zero, indicating the complete erosion of its basin from the system. This phenomenon is observed in the basin portrait of Fig. 8d, where the blue basin representing A-HP 1 solution completely dominates the attractive region. Physically, the disappearance of A-LP 1 solution marks the transition from a bi-stable to a monostable cantilever response in the attractive regime. Interestingly, in Fig. 8a, by considering an even blunter tip, with radius in the range of R ∈ [114, 117] nm, the robustness of A-HP 1 suddenly drops, due to the appearance of a novel competing attractor in the system. The new attractor is indicated by the light brown basin in Fig. 8e. This novel attractor, which is the in-contact attractor, has a smaller growth rate at lower radius values and does not affect the erosion of the A-HP 1 solution rapidly. This is observed in Fig. 8a by a steady increase in the integrity measure of A-HP 1 solution between radius values R ∈ [117, 185] nm. With the further increase in the blunting of the tip, the A-HP 1 basin is eroded along its boundaries smoothly by the in-contact attractor as seen in Fig. 8f. The robustness characteristics shown by LIM and IF measures for attractive regime follow similar trend. However, the IF safe basin measure is larger in magnitude compared to LIM which is due to the smooth erosion of the basin without in-well fractality at low excitation amplitudes (ȳ < 0.015).
The repulsive attractor, unlike attractive, shows a steady growth in basin size for increasing radius values as shown in Fig. 8b. This trend of increasing basin size is due to the fact that an increase in radius value will increase the area over which repulsive forces are perceived by the system. However, it is worth to observe in Fig. 8d-f that the size of the purple basin (strength of integrity measure) associated with repulsive attractor is very small compared to the attractive and the in-contact regimes. Thus, at lowȳ, the repulsive attractor although resilient to changes in radius has a smaller influence on the basin erosion process as compared to other two attractors. Moreover, at low radius values, a perturbation inside the attractive regime (blue basin) can lead the system towards the repulsive attractor leading to hardening behaviour as illustrated in Fig. 8g. In this case, for a cantilever oscillating initially in the attractive region, the system frequency response shows both softening due to attractive forces and hardening due to repulsive forces as shown in Fig. 8g (black frequency response curve). However, with a further increase in radius values, the attractive and repulsive basins are separated and the system oscillates purely in attractive regime showing only softening nonlinearity (green frequency response curve in Fig. 8g).
The robustness characteristics displayed by the incontact attractor in Fig. 8c are unique, and present features are not observed in the other two attractors. For low radius values, e.g. R < 117 nm, the in-contact attractor does not exist as seen in Fig. 8c and further illustrated by the absence of light brown basin in Fig. 8d. This highlights the existence of a critical radius value for the manifestation of the in-contact attractor. This critical radius value for our system is at R = 117 nm. Around this radius value, there is a sudden appearance of the in-contact attractor (shown in Fig. 8e as light brown basin) and it grows steadily with the increase in the value of R. Interestingly, the critical radius value remains the same for higher excitation amplitudes and excitation frequencies. Furthermore, from Fig. 8c it is observed that, the growth rate of low-amplitude solution (LP 1 ) is not as steep as the attractive and repulsive attractors, but at large radii, the in-contact attractor eventually becomes dominant with no in-well fractality or boundary erosion. Figure 9 illustrates the evolution of the basin portrait with respect to the tip radius. The LIM integrity measure is utilized to characterize the robustness of the attractors and track the changes in the basin portraits as the tip deteriorates. The basin portraits are analysed for constant parametersȳ = 0.020 andΩ = 0.8. The figure reinforces the previous discussion pictorially. It depicts the disappearance of A-LP 1 basin, fractalization of the attractive regime triggering the erosion process and finally the requirement of critical radius for manifestation of the in-contact attractor.

Basin erosion as a function of excitation frequency and excitation amplitude
In Sect. 4.1, the evolution of phase-space topology as a function of radius was showcased. Accordingly, in this section the dynamical integrity analysis aims at quantifying the extent and evolution of the basins along with their erosion process as a function of the excitation amplitude (ȳ) and excitation frequency (Ω). This is established in Fig. 10 for excitation amplitudes y = 0.005,ȳ = 0.010 andȳ = 0.020, respectively, whereas Figs. 11 and 12 report the snapshots of basin portraits at crucial excitation frequencies (Ω) near the neighbourhood of fundamental and principal parametric resonance frequencies, respectively. All simulations are performed at a constant radius value of R = 150 nm and finally, LIM and IF measures are calculated for the aforementioned data. However, we observed that the strength of LIM and IF measure as a function ofΩ remain approximately the same. Thus, we can argue that the excitation frequency in the selected interval does not modify the global shape of the portrait and we do not experience a significant subdivision of the basin. Therefore, in order to simplify the analysis, only LIM is used to quantify the robustness of dominant attractors present in the system.

Analysis of basin erosion profiles around the fundamental resonance frequency
The fundamental ideology of TM-AFM is based on the near-resonant excitation of the microcantilever. Therefore, it is of significant interest to study the system topology with regard to basin erosion profiles since, these profiles are indicative of solution instabilities. In this respect, the following section utilizes the integrity measures shown in Fig. 10 to discuss the various system attractors in the neighbourhood of their respective resonance frequencies.
In addition, to further delineate the behaviours observed in the integrity profiles, the basin portraits are reported in Fig. 11 at specific frequencies.
The integrity profiles of attractive regime attractors are illustrated in Fig. 10a and, similar to Sect. 4.1, the crimson and blue colours belong to the low-amplitude non-resonant solution (A-LP 1 ) and the high-amplitude resonant solution (A-HP 1 ), respectively. Figure 10a shows that the attractive regime is dominated by A-LP 1 non-resonant solution at low excitation frequen-ciesΩ < 0.77. This is illustrated through the basin portrait of Fig. 11a in which the crimson basin corresponding to A-LP 1 is the dominant solution. By further increasing the excitation frequency value abovē Ω > 0.77, the A-LP 1 solution exhibits a sharp decrease in its integrity value (illustrated in Fig. 10a). The sharp decline is attributed to the sudden appearance of the A-HP 1 resonant attractor from inside the local potential well via saddle-node bifurcation. The appearance of saddle-node triggers the erosion process of the A-LP 1 basin from outside its boundaries. The appearance of A-HP 1 (blue basin) inside the A-LP 1 local potential well (crimson basin) is shown in Fig. 11b. Any further increment in excitation frequency causes the complete erosion of the low-amplitude attractive (A-LP 1 ) solution as seen in Fig. 11c, leaving large amplitude oscillations of A-HP 1 as the dominant solution in the attractive potential well. Hereafter, the A-HP 1 solution remains largely stable as illustrated in Fig. 10a, and its robustness is mainly affected when the excitation frequency reaches the in-contact fundamental resonance or when the system is driven with an excitation amplitude above the parametric threshold of the system. This trend of solution instabilities observed around the fundamental resonance frequency remains unaffected irrespective of excitation amplitudes.
Contrary to the above discussion on attractive regime, where steady-state solutions exist over a large range of frequencies, the repulsive basin (purple) exists in a narrow frequency range around the fundamental resonance frequency and shows sharp increase/decrease as we move closer/farther away from resonance. This is illustrated in Fig. 10b forȳ = 0.005, y = 0.010 andȳ = 0.020, respectively. Moreover, the repulsive basin size is small compared to other two attractors in case ofȳ = 0.005 as shown in Fig. 11c. But it displays a sharp increase in size for higher values ofȳ as seen in Fig. 11d forȳ = 0.020. This is due to the fact that the harder the cantilever is driven, the deeper the oscillations penetrate into the repulsive regime, and the time period of oscillations spent in repulsive regime increases. Therefore, contrary to the observation in Sect. 4.1, the repulsive basin at higher excitation amplitudes significantly constricts the attractive regime basin in the neighbourhood of the fundamental resonance frequency. On the other hand, the in-contact regime solution displays rich nonlinear behaviour absent in the case of the attractive and repulsive regimes. Figure 10c-e illustrates the integrity profiles of various solution branches namely LP 1 , Su 1 -HP 1 , and HP 1 that are observed in the in-contact regime. Similar to Sect. 4.1 and summarized in Table 4, the light brown colour corresponds to nonresonant LP 1 solution branch, the green colour belongs to Su 1 HP 1 solution branch and orange colour belongs to HP 1 resonant solution branch. The appearance and  Table 4 disappearance of several rare attractors, together with period-2 responses appearing not only on HP 1 solution branch but also on Su 1 HP 1 superharmonic solution causes the robustness of in-contact attractor to vary rapidly. This behaviour can be observed in the integrity profiles of Fig. 10d, e in the form of sharp peaks and valleys as the excitation amplitude is increased. At low excitation frequenciesΩ < 1.5, the incontact basin is dominated by LP 1 solution as shown in Fig. 10c. This is further demonstrated in the basin portrait of Fig. 11a where only the light brown basin dominates the in-contact region. Further increasinḡ Ω > 1.5, we observe the drop in integrity of LP 1 solution due to the appearance of Su 1 HP 1 attractor. This is visualized by comparing Fig. 10c, d between frequency rangesΩ ∈ [1,2]. The drop in integrity is sharp for small values ofȳ and slowly smoothens with increasinḡ y amplitudes. This effect is due to the increased softening effect of the in-contact response which causes the Su 1 HP 1 branch to overlap with LP 1 branch over larger Ω values. Along with the increased overlapping effect, at higher amplitudes of excitationȳ > 0.00118, the Su 1 HP 1 solution grows more robust. The appearance of Su 1 HP 1 attractor (green basin) from within the LP 1 basin forȳ = 0.020 is shown in Fig. 11e, f. Furthermore, withΩ growing closer to resonance (i.e.Ω = 3.72), the high-amplitude resonant solution (HP 1 ) appears via saddle-node bifurcation and grows in dominance as depicted in Fig. 10e. In a similar fashion to the attractive regime, the saddle-node triggers the erosion of the previously dominant Su 1 HP 1 solution from outside the boundary. This is observed in basin portraits from Fig. 11g-i where the orange basin belonging to HP 1 solution is growing along the periphery of green basin. Interestingly, at low amplitudes of excitation,ȳ < 0.010 the erosion of Su 1 HP 1 basin is smooth with no influence of rare attractors on its robustness. After the complete erosion of Su 1 HP 1 basin (see Fig. 11i), the in-contact response is completely dominated by HP 1 solution untilΩ reaches close to parametric resonance frequency. This is illustrated by the gradual drop in robustness measure in Fig. 10e for values ofΩ > 3.7 . The increase inȳ has a peculiar effect on the HP 1 solution since the robustness decreases in contrast to the expected increasing trend. This is due to the large number of rare attractors which tend to appear at higher excitation amplitudes. This peculiarity makes the experimental investigation of the in-contact attractor a challenging task.

Analysis of basin erosion profiles around the principal parametric resonance frequency
The aforementioned discussion was focused on the robustness and erosion profiles of system attractors for excitation amplitudes below the parametric threshold and excitation frequencies around the fundamental resonance. However, the system excited parametrically exhibits different dynamics that are not observed through direct excitation. In addition, the theoretical analysis of parametrically driven AFM has shown the added benefits such as high quality factor, lower imaging forces and reduced cantilever transients [33]. These advantages, if harnessed, can be of significant interest in areas such as soft polymers and biological specimens. Therefore, the current section focuses on the basin erosion profiles of system attractors in the neighbourhood of respective principal parametric resonance frequencies. Similar to the previous section, the analysis utilizes integrity measures illustrated in Fig. 10 to discuss the evolution of system responses, whereas Fig. 12 is used to understand the metamorphoses of basin erosion graphically. Considering the attractive force regime, an excitation amplitude ofȳ = 0.020 is utilized to study the robustness of attractors oscillating near the parametric resonance frequency. The chosen excitation amplitude is above the required threshold amplitude ofȳ = 0.01515 needed to excite the system parametrically. The solution belonging toȳ = 0.020 case is drawn by dotted lines in Fig. 10a. From Fig. 10a there is a sudden decline in the A-HP 1 attractor's integrity value around the parametric resonance frequency of Ω ∈ [1.4, 1.75]. This sudden decrease in robustness is promoted by the appearance of period-2 attractor (A-Pa-P 2 ) within the compact part of A-HP 1 basin. This is illustrated in the basin portraits of Fig. 12a-c, where the A-Pa-P 2 attractors shown by dark red basin are surrounding the A-HP 1 blue basin. The rate of erosion of A-HP 1 basin is directly proportional to the nearness ofΩ to the principal parametric resonance frequency. For values ofΩ away from the parametric resonance the A-HP 1 integrity shows a steady increase back to its original value. Hereafter, the A-HP 1 solution remains almost independent ofΩ until the in-contact resonance frequency is reached.
In case of the in-contact regime, the parametric threshold amplitude is found to be as low asȳ = 0.00014. The low threshold amplitude also suggests  Table 4 the potential application of stochastic resonance of cantilever oscillating in the in-contact regime. However, in order to facilitate the study, an excitation amplitude of y = 0.005 is utilized. The corresponding integrity profile is illustrated in Fig. 10e by the continuous line, and the relative basin portraits are shown in Fig. 12d, e. By considering Fig. 10e, the integrity measure of HP 1 solution displays a steady decline asΩ is brought close to principal parametric resonance frequency (Ω = 7.44). Finally, aroundΩ = 7.1 the HP 1 basin completely disappears, indicating that the system is in the parametric instability region. The appearance of paramet-ric period-2 attractors (Pa-P 2 ) in the basin portrait is indicated by the dark red basin in Fig. 12d, e. The figures illustrates the appearance of period-2 attractors (dark red) from the compact part of HP 1 basin (orange) along with period-6 rare attractors (dark blue). Finally, the erosion process of HP 1 solution is non-smooth with the interacting basins featuring severe fractality.

Conclusions
The global dynamics of TM-AFM has been investigated with the aim of evaluating the robustness and dynamical integrity of coexisting attractors of the system. Extensive numerical analyses have been carried out to show the existence and the properties of the incontact attractor. The frequency response curves and basin portraits are obtained, showing that the in-contact attractor is highly sensitive to the main driven parameters with large amplitude superharmonic branches appearing for small excitation amplitudes. In order to unveil the entire bifurcation scenario of the so-called in-contact attractor around the primary and parametric resonance, several local bifurcation envelopes are combined together to build global bifurcation maps. Utilizing the bifurcation maps, the escape thresholds along with various response scenarios in the excitation parameter space are analysed in detail around the direct and parametric resonance frequency. The outcome of the analysis shows new routes to crisis, escape scenarios via appearance of strange attractors and multiple period-doubling cascades. In addition, the robustness of attractors has been analysed by making use of basins of attraction and integrity measures such as local integrity measure (LIM) and integrity factor (IF). The analysis has focused on the basin erosion with respect to variation in excitation frequency, excitation amplitude and the AFM probe tip radius. The results highlight the appearance of in-contact attractor for a critical radius value. In addition, the parametric resonance and its effect on basin erosion via fractalization are discussed for both attractive and in-contact attractors. It is seen that the period-2 attractors arising from parametric resonance decrease the robustness of attractive regime in a smooth fashion, whereas in case of in-contact attractor the period-2 solution together with higher-order strange attractors, erodes the basin through fractalization. In conclusion, our analysis of basins of attraction, global bifurcation charts and integrity profiles provides a method to study the complex dynamics involved in TM-AFM. , Appendix B: Instructions on reading the bifurcation chart of Fig. 6 The instructions on reading the complex bifurcation and response scenarios showcased in Fig. 6 are explained below with the help of a frequency response curve described in Fig. 5. In order to assist with easier understanding, the aforementioned figures are reintroduced in this section and marked with coloured circles in Fig. 13 and coloured lines in Fig. 14 at specificΩ values. These coloured lines depict regions of interest. Furthermore, Table 5 delineates the various dynamic response regions seen in Fig. 14.
Referring to Fig. 13, the blue lines which represent the saddle-node bifurcation (SN) loci divide the bifurcation chart into several regions of period-1 solutions such as R1, R2, R4, R5, R6, R7, and R9. Furthermore, the aforementioned saddle-node loci meet at specific values of excitation parameters in pairs of two and give rise to new multi-stable regions atΩ = 1.24 (intersection of R5 and R10),Ω = 1.851 (intersection of R2 and R4), andΩ = 3.721 (intersection of R1 and R9). The dynamic behaviour of system changes depending on whether the excitation amplitude is above or below the point of intersection of saddle-node loci. For instance, in Fig. 13 atΩ = 1.851 region R2 exists above and region R4 below the intersection point.
In particular, special attention should be given to regions where period-2 solution exists. These regions are bounded by red and green loci corresponding  to period-doubling/flip bifurcations. In general the period-2 solutions are bound to narrow regions over a small range of excitation parameters (except for parametric resonance case). These narrow regions invade into larger regions of period-1 solutions resulting in new intersection regions where combination of period-1 and period-2 solutions exists. An example of such an intersecting region is indicated in Fig. 13 as R6 + R3 where HP 1 -P 2 and Su 1 HP 1 solutions coexist. Finally, by utilizing the aforementioned instructions the bifurcation behaviour observed in Fig. 5 is described below. An excitation parameter equal toȳ = 0.020 is utilized in simulating Fig. 14. This value is further highlighted in Fig. 13 with a light brown horizontal line spanning across the entire range of the excitation frequency. Referring to Fig. 14, the first region of interest is marked by a red vertical line and spans region up toΩ ≤ 0.93. From Fig. 13, this region starts with LP 1 solution marked by R1 and ends with LP 1 + Su 1 -P 2 (R1 + R8) solution atΩ = 0.93. The second region of interest is between the red and green vertical lines in Fig. 14 and ranges fromΩ ∈ (0.93, 1.236]. From Fig. 13, this region starts with LP 1 + Su 1 HP 1 solution indicated by R5 fromΩ > 0.93 and ends with multi-stable region R7 consisting of LP 1 + Su 1 HP 1 + Su 2 HP 1 solutions atΩ = 1.236. Similarly, the third region can be found between the green and blue vertical lines in Fig. 14 and ranges fromΩ   ∈ (1.236, 1.795]. This region shows the presence of Su 1 HP 1 + Su 2 HP 1 solutions marked by R10 from Ω > 1.236 up toΩ = 1.795. The next region of interest is bounded by blue and purple vertical lines in Fig. 14 and ranges fromΩ ∈ (1.795, 2.413]. Referring to Fig. 13, forΩ > 1.795, the region exhibits monostability consisting of only Su 1 HP 1 solution marked by region R6 and ends with region R8 consisting of Su 1 -P 2 solution atΩ = 2.413. Continuing along the same trend, the other sections highlighted in Fig. 14 can be traced and are shown in Table 5.