Local Versus Global Dynamics and Control of an AFM Model in a Safety Perspective

The role of local and global dynamics to assess a system robustness and actual safety in operating conditions is investigated, by also studying the effect of different local and global control techniques on the nonlinear behavior of a noncontact AFM. First, the nonlinear dynamical behavior of a single-mode model of noncontact AFM is analyzed in terms of stability of the main periodic solutions, as well as attractors robustness and basins integrity. To the same AFM model, an external feedback control is inserted during its nonlinear continuum formulation, with the aim to keep the system response to an operationally suitable one. The dynamical analysis of the controlled system is developed to investigate and verify the effects of control into the system overall behavior, which could be unexpectedly influenced by the local nature of the control technique. A different control technique is finally applied to the AFM model, acting on global bifurcation events to obtain an enlargement of the systems safe region in parameters space. The analytical procedure, based on Melnikov method, is applied to the homoclinic bifurcation involving the system hilltop saddle, and its practical effects as regards possibly increasing the system overall robustness are numerically investigated by means of a dynamical integrity analysis. Then, a fully numerical procedure is implemented to possibly control global bifurcations involving generic saddles. The method proves to succeed in delaying the drop down of the erosion profile, thus increasing the overall robustness of the system during operating conditions.


Introduction
Analyzing the nonlinear dynamical behavior of a mechanical system usually entails the investigation of the stability of its response under variation of some (characterizing and/or controllable) parameter. Bifurcation diagrams and eventually behavior charts in parameters space allow to verify the system robustness to possible changes of the operational parameter setup, and often reveal a wide variety of nonlinear phenomena, such as bifurcations and in-well instability regions which can lead to possible unstable, aperiodic or chaotic oscillations. However, besides these analyses, possible changes in the system initial conditions due to imperfections have to be taken into account, since it is nowadays ascertained that the safe operation of a nonlinear system depends not only on the local stability of its solutions but also on the global dynamics associated with the uncorrupted basin surrounding each solution. It appears therefore evident that tools and concepts of global dynamics and dynamical integrity represent crucial instruments not only to achieve a complete and accurate description of the system dynamics, but, from a practical viewpoint, to assess the actual safety of a dynamical system (Rega and Lenci, 2015). More importantly, they can be proposed as a comprehensive approach to be applied for an in-depth investigation of a nonlinear model, in order not only to theoretically analyze the robustness of competing attractors and the erosion processes that bring to the escape from bounded regions, but also to critically evaluate the effect of different control techniques and, from an operational viewpoint, to furnish hints useful for engineering design. Following these guidelines, this Chapter aims to present a general framework for studying local and global dynamics, as well as control, of a sample mechanical system, taking as reference a reduced order model of noncontact Atomic Force Microscope (AFM), which can indeed represent a large number of nonlinear models with some ensuing dynamic phenomena.
AFMs are powerful devices used for surface analysis in nano-electronics, mechanics of materials and biotechnology, as they permit to topologically characterize surfaces up to micro and nano resolution levels (Sarid, 1991;Morita et al., 2009;Eaton and West, 2010). In a typical AFM, the topography is imaged by scanning a sharp tip, fixed to the free end of a microcantilever vertically bending over the sample surface, and by measuring the tip deflection through a laser technology. The tip-sample interaction modifies the beam dynamics and allows not only to image surfaces, but also to measure some physical properties of the sample. As far as the AFM dynamics is concerned, the most common operation modes are the tapping mode, in which the tip operates in both attractive and repulsive force regions and touches the surface only for short time intervals, and the noncontact mode, in which there is absence of contact between the tip and the sample and their interaction is governed by a solely attractive potential. Consequently, for the latter type of AFM the tip has to maintain a design gap from the sample such to ensure that the beam elastic restoring force is stronger than the atomic attraction. Otherwise, instability of the equilibrium configuration occurs, with the so called "jump to contact", or escape (in dynamical systems terms), phenomenon. As the dynamic excitation tends to strongly reduce the equilibrium gap, the study of system stability as a function of a varying excitation amplitude, among other physical parameters, is a very important issue for noncontact AFMs, in order to reliably determine the escape threshold separating the region of bounded (i.e., noncontact) solutions from that of unbounded solutions, the latter corresponding to unwanted contacts between tip and sample.
During the scan operation, moreover, the sample roughness can modify the distance between the microcantilever tip and the sample to be scanned, and thus, the nonlinear atomic force interaction which is used to obtain the topography can lead to unstable dangerous motions. To prevent these undesirable motions and improve the microscope performances, several control techniques have been proposed in the field of AFMs during the last two decades, primarily based on the feedback control methods that work by keeping the microcantilever vibration to a selected reference one and allowing to simultaneously and reliably measure the sample surface. However, since these procedures are focused on controlling the system local dynamics, their impact on the system overall dynamics is generally unknown. Nevertheless, analyzing the latter is of great importance, especially for micro/nano mechanical systems like AFMs for which slight changes of the initial position and/or velocity at the nanoscale level, or modifications of the operational parameter setup, can produce dramatic modifications of the overall dynamics.
From a global dynamics perspective, it is thus of interest to clearly detect in the state plane the basins of attraction of the periodic (acceptable) solutions and that of the (undesirable) unbounded response compromising the device operation. Furthermore, the relevant erosion process due to variations of some system parameter can be followed and quantified by means of the construction of the so-called erosion profiles. As a practical consequence of such analyses, which belong to the field of dynamical integrity, the sensitivity of a system to variations of both operational parameters and initial conditions can be discussed, and some hints useful to define thresholds able to ensure acceptable safety targets can be achieved. The concept of dynamical integrity can be also referred to for developing a different control technique based on the global properties of the dynamical system, which acts on the homo/heteroclinic bifurcations involving the stable and unstable manifolds of the system saddles responsible for the erosion of the basins. It operates by properly shifting them thanks to the addition of some controlling superharmonics to the reference harmonic excitation, thus increasing the system overall robustness.
The Chapter is organized as follows. Nonlinear local and global dynamics of a reduced order model of noncontact AFM are presented in Section 2, while Section 3 investigates the dynamical effects of an external feedback control on the overall response of the AFM model. Section 4 deals with the implementation, into the same model, of a global control technique, whose effects are analyzed and critically discussed also in terms of safer AFM operation. Finally, in Section 5 some general conclusions are drawn.
2 Nonlinear dynamics of a reduced order model of noncontact AFM 2.1 Equations of motion, reduced order model and unperturbed system The physical model under consideration is a fixed-free Atomic Force Microscope (AFM) microcantilever, which is assumed to be planar, inextensible and horizontal, with length L and a sharp tip of height h T close to its free end, and with a distance g between its fixed side and the sample ( Fig. 1(a)) (Hornstein and Gottlieb, 2008). The beam material is considered linearly elastic, ho-mogeneous and isotropic, with Young modulus E. Reference is made to the general formulation based on the classical inextensional beam model of Crespo da Silva and Glynn (1978). Accordingly, the set of two coupled partial differential equations (PDEs) for the beam horizontal and vertical transverse vibrations is whereū(r, t) andv(r, t) are the horizontal and vertical displacements and subscript letters denote partial differentiation with respect to arc length r and time t. Coefficients EI, J z and m are beam stiffness, principal moment of inertia and mass per unit length, respectively, and Λ is a Lagrangian multiplier accounting for the inextensibility condition (1 +ū r ) 2 +v 2 r = 1. Generalized forces in horizontal and vertical direction are represented byQ u ,Q v , the former corresponding to a feedback control force depending on the horizontal displacement and the latter also accounting for the localized (at r = a T , see Fig. 1(a)) transverse atomic force interaction derived from a Lennard Jones potential for a sphere-plane system (Israelachvili, 1992;Sarid et al., 1996), localized at the tip and representing the interaction between the tip and the sample. Since the AFM under analysis operates in noncontact regime, the most important atomic force interaction is the attractive one and the repulsion interaction can be neglected. The set of non-homogeneous boundary conditions v(0, t) =V (t),v r (0, t) = 0,ū(0, t) =Ū (t), v rr (L, t) = 0,v rrr (L, t) = 0,ū r (L, t) = 0 ( 2) completes the formulation of the problem, withV (t) andŪ (t) being the prescribed vertical transverse and horizontal scan displacement, respectively. The incorporation of the holonomic constraint, which allows to express the horizontal displacement in terms of the vertical one, the isolation and expansion of the Lagrange multiplier and the nondimensionalisation and use of a moving reference frame (v(s, τ ) = w(s, τ ) + V (τ )) lead to the formulation of an initial-boundary value problem (IBVP) with a homogeneous set of boundary conditions for the transverse vibrating motion. A single mode assumption (w(s, τ ) = q 1 (τ )φ 1 (s)) and a Galerkin approximation (the basis function being that of a clamped-spring beam) reduce the IBVP to the second order ordinary differential equation (ODE), which has the following nondimensional form: 1 + α 2 x 2 ẍ + α 1 + α 2ẋ 2 + α 3 x 2 x = − Γ 1 (1 + x + V g ) 2 − ρ 1 + ρ 2 x 2 ẋ − ν 2 V g + ν 1Vg + µ 1 x + µ 2 x 3 Ü g + η 1Ug + η 2 U g where x(τ ) = q 1 (τ )φ 1 (α)/γ and t N = ω 1 τ ; α and γ are the nondimensional tip distances from the microbeam fixed end and the sample, rispectively; ρ 1 is the damping coefficient; Γ 1 is the atomic attraction parameter; α 1 , α 2 and α 3 are coefficients of the linear and nonlinear terms; Figure 1. Microcantilever model. and V g = V /γ represent the horizontal and vertical excitations; η 1 , η 2 and ρ 2 are internal feedback control parameters, the first two related to the time dependent horizontal scan and the latter to the cubic damping term. For the nondimensionalisations and detailed expressions of all coefficients, see Hornstein and Gottlieb (2008).
It is worth underline that the choice to refer to a single-degree-of-freedom model is here justified by the operating conditions dealt with in the following. In fact, even if it has been shown that for tapping AFMs a multimode Galerkin approximation is needed to detect nonlinear phenomena (e.g., grazing bifurcations) that a single-mode analysis does not match, however, in the noncontact operation range, a multimode discretization does not enrich the system response, and a singlemode approximation is sufficient to detect the main nonlinear aspects (Bahrami and Nayfeh, 2012). Moreover, the investigated frequency range spans around primary and secondary resonance of the first mode, where the contribution of the higher modes is substantially negligible (Hornstein and Gottlieb, 2012).
Accounting for the orders of magnitude of various coefficients in commercial AFMs, feedback controls and the nonlinear terms related to α 2 can be neglected, to obtain Note that the AFM scan process is performed by means of both vertical and horizontal excitations, the former allowing to quantify the interaction forces and the latter being necessary to obtain the three dimensional map of the sample. Yet, in existing commercial AFMs the horizontal scan (parametric) frequency is much smaller than the vertical (external) one, so that it appears reasonable to study the two forced cases in the neighborhood of resonances separately. The results presented in this chapter will thus refer to the sole parametric horizontal scan, which is supposed to be harmonic (U g = U sin(ω u t)), even if few considerations will be drawn about the role played by the vertical external excitation, as well.
The operation domain of noncontact AFMs must be such to avoid jump to contact with the scanned sample. In dynamical system terms, this is ascertained by considering the undamped, unforced version of Eq. (4)ẍ with the relevant Hamiltonian, where y =ẋ, being with the associated single, asymmetric, potential well with left (i.e., towards the sample position x = −1) contact direction ( Fig. 2(a)). The unperturbed state space is depicted in Fig. 2(b), where the two fixed points of the time-independent problem, i.e. the stable equilibrium (E) of the cantilever tip under elastic (α 1 , α 3 ) and atomic interaction (Γ 1 ) forces and the corresponding hilltop saddle (S H ), are reported for a given set of values of the governing parameters. The homoclinic orbit y h (t) of the saddle is also plotted in Fig. 2(b): it separates the inner region of bounded periodic solutions from the outer region of unbounded solutions, the former representing the safe domain for noncontact AFM operation. The study of the fixed points as a function of the atomic interaction coefficient Γ 1 is reported in Fig. 2(c) for α 1 = 1 and α 3 = 0.1, and provides the upper boundary for the stable equilibrium E existence at the limit value Γ 1 = 4/27 where a saddlenode bifurcation leads to the birth of the unstable branch. The lower stable equilibrium E 1 settles under the limit value x = −1, which corresponds to the position of the sample with respect to the cantilever tip, thus having no physical meaning. This entails that for Γ 1 values higher than that of the saddle-node the microbeam "jumps to contact" with the sampled surface. This phenomenon, associated with the escape of the response from the single potential well of Fig. 2(a), is precluded in noncontact AFMs: so, there is a great interest in analyzing the conditions for its impending , and unperturbed phase space (b) for α 1 = 1, α 3 = 0. occurrence in different ranges of frequency around parametric (and/or external) resonances, along with the features through which they are realized when varying a control parameter, typically an excitation amplitude characteristic of the system.

Bifurcation scenarios, response charts and escape threshold
The nonlinear response of the parametrically excited single-mode model is analyzed via continuation techniques and numerical simulation (Doedel and Oldeman, 2012), for the following set of parameter values (Rega and Settimi, 2013): The calculated natural frequency for these values is ω 1 = 0.835. Several bifurcation diagrams as a function of the forcing amplitude U have been obtained in a large range of frequencies including the fundamental (ω u = ω 1 ) and principal (ω u = 2ω 1 ) parametric resonances, and the main periodic solutions and local bifurcations have been detected. These diagrams exhibit a variety of response scenarios. The local bifurcation loci in the excitation parameter control space (forcing frequency vs. forcing amplitude) are summarized in the semi-logarithmic chart of Fig. 3. The system overall escape threshold separates the bounded solutions (below the curve) from the unbounded solutions (above the curve) and is obtained as the envelope of local bifurcation escape thresholds in different parameter ranges. From a physical viewpoint, it represents the (unacceptable) amplitude value that would bring the beam tip oscillation beyond the location of the sample (at x = −1); from the dynamical viewpoint it corresponds to the total annihilation of all basins of attraction (see Sect. 2.3 forward). Changes in the escape threshold slope correspond to changes in the kind of bifurcation event leading to escape and, apart from a localized exchange of the governing one in the frequency range between 0.4 and 0.5, four main different regions can be identified in the parameter control space (Fig. 3), from the escape viewpoint. Region I includes low frequencies set on the left (ω u = 0.5 − 0.72) of the downward vertex A of the overall escape threshold corresponding to nonlinear fundamental resonance. Here, the system displays coexistence, for low values of the forcing amplitude, of two stable 1-period solutions, with the initial Low-amplitude (nonresonant) P1L solution being connected to the (resonant) P1H solution of High-amplitude through the classical unstable branch in between a couple of saddlenode bifurcations. In this range the nonresonant solution governs the escape from the bounded region, which thus is determined by the occurrence of the relevant saddle-node bifurcation, while the resonant response exists in a limited range of forcing amplitude and becomes unstable through a supercritical period doubling. As the frequency approaches the nonlinear resonance (ω u = 0.72, A-vertex), the latter bifurcation moves to growing amplitude values, until at the right of the nonlinear resonance (region II, ω u = 0.72 − 1.12) it becomes the local event triggering escape, as shown in the sample bifurcation diagram reported in Fig. 4. Escape occurs after a sequence of period doublings ending up to chaos and disappearing through a successive boundary crisis, all of this being not represented in Fig. 3 because of occurring in very narrow ranges of the control parameter. The III and IV regions are characterized by the presence of a 2-period solution, typical of the dynamical behavior around the subharmonic resonance, which arises through a subcritical period doubling of the 1-period response (SbPD1) and disappears by means of a supercritical period doubling. In this respect, thus, the role played by the saddle-node bifurcation of the lowamplitude 1-period solution at left of the A-vertex is herein replaced by this subcritical period doubling bifurcation. The 2-period solution becomes stable via saddle-node (SN2) bifurcation at low values of forcing amplitude; the response pattern is now similar to the one in region I, apart from replacing P1H with P2.
As in that case, the range of coexistence of stable P1 and P2 solutions grows up with increasing forcing frequency, up to entailing (at ω u ∼ = 1.56, B-vertex) the exchange of the local event triggering overall escape from SbPD1 to the supercritical period doubling bifurcation (SpPD2) of P2. This is substantially the same pattern as the one on the right of the nonlinear fundamental resonance and, indeed, the cascade of period doubling bifurcations  arising from SpPD2 keeps governing the system escape in the whole region IV, which includes the principal resonance range (ω u = 1.56 − 1.8, see the bifurcation diagram of Fig. 5).
The bifurcation diagrams of Figs. 4-5 show also the short-range coexistence with the main periodic solutions of variable solutions of higher periodicity, which arise from saddle-node bifurcations and end up with local chaotic responses via series of period doublings occurring in very narrow ranges of amplitude values. These solutions are indeed of minor interest.
Focusing around the fundamental resonance ( Fig. 6), the system shows the same qualitative behavior of a number of softening oscillators subjected to primary external excitation, such as the   Helmholtz oscillator (Szemplinska-Stupnicka, 1992) and single-mode models of MEMs (Lenci and Rega, 2006;Alsaleem et al., 2010), especially for what concerns the V-shaped region of escape, its limiting boundaries, and the underlying triangle region with the two coexisting 1-period solutions. Figure 6 shows that such a coexistence in the region in between the two SN loci occurs up to ω u = 0.835 = ω 1 , where they collapse with each other. Right of this value, and above the corresponding value of forcing amplitude, the sole 1-period solution previously associated with P1H occurs (now called P1 solution) and its annihilation through SpPD (and the following series of period doublings) characterizes the smooth transition to escape for increasing amplitude. Note again that no loci of supercritical period doubling of higher periodicity solutions are reported in Fig. 6, left of (and above) the SpPD threshold, nor the ensuing locus of boundary crisis leading to escape, since these events occur in a very narrow range of control parameter values. Note also that, in contrast, transition to escape from the left side occurs via the typical "sudden" SN bifurcation.
For the sake of completeness, it is worth noting that the inspection of the bifurcation/response scenarios at low values of the forcing frequency ( Fig. 7) points out some interesting features of the system dynamical behavior. The global escape threshold, in fact, displays local minima at frequencies values corresponding to superharmonic resonances ω u ∼ = 0.41 = ω 1 /2 (also visible in Fig. 3), ω u ∼ = 0.27 = ω 1 /3 and ω u ∼ = 0.21 = ω 1 /4, where, as already seen at fundamental resonance, the dynamical response is characterized by a recurrent triangle region of coexistence of nonresonant P1L and resonant P1H solutions, and, according to the resonance frequency value, such solutions are now periodic responses with two harmonics of frequency ω and 2ω (3ω,4ω). Soon after and to the right of these triangle regions, the system displays the presence of two thresholds of period doubling which delimit a narrow region around ultrasuperharmonic frequencies at ω u = 0.218−0.229 ∼ = 2ω 1 /7, ω u = 0.301−0.326 ∼ = 2ω 1 /5 and ω u = 0.47−0.55 ∼ = 2ω 1 /3. Inside these ranges, the sole P1 solution becomes unstable via a subcritical period doubling which leads to the birth of a P2 solution, whose frequency components include, besides the ω/2 frequency characterizing the period-doubled solutions, the superharmonic frequency and its higher harmonics. Thus, summarizing the general dynamical behavior of the model, it can be stated that at fundamental and superharmonic resonance frequencies (i.e. ω u = ω 1 /n, n = 1, 2, 3, 4...) the system response is characterized by coexistence of resonant and nonresonant 1-period solutions, while at principal resonance and ultrasuperharmonic resonance frequencies (i.e. ω u = 2ω 1 /n, n = 1, 3, 5, 7...) the main and partially coexisting periodic responses are 1-period and 2-period solutions.
Finally, it is of interest to analyze the influence of the nonlinear interaction parameter variation Γ 1 on the system dynamical response; such parameter, which depends on the kind of tip and sample materials and their distance at nanoscale level, is in fact the characterizing ingredient of an AFM model, and introduces a nonlinear term of order -2 into the system (see Eq. (4)). The behavior chart in the Γ 1 − U plane reported in Fig. 8 for a forcing frequency ω u = 0.7 near the fundamental resonance shows that as the nonlinear interaction increases, with respect to the reference value Γ 1 = 0.1 used for the previous numerical analyses, the amplitude escape value for the P1L solution rapidly decreases, up to the disappearance of such periodic solution and consequently of the region of coexistence of resonant P1H and nonresonant P1L solutions (at Γ 1 ∼ = 0.13); after that, for higher values of the interaction parameter the sole P1H solution remains as system stable periodic solution, and its escape boundary represented by the supercritical period doubling threshold PD1H moves to higher values of the forcing amplitude U , enlarging the stability region up to Γ 1 ∼ = 0.14, when it drastically falls down causing the annihilation of the P1H stability region. The outcomes of Fig. 8 can be exploited for design purposes by referring to a dual reading key: on the one hand, in fact, they point out the crucial role played by the tip-sample nonlinear interaction on the response of the AFM cantilever, furnishing practical information about the limiting values of forcing amplitude to be used depending on the sample constitutive properties (i.e., with a given value of Γ 1 ), while on the other hand they can be used to calibrate the tip-sample interaction (e.g., tip material choice, or material of the sample to be possibly scanned) depending on the AFM operation settings (e.g., excitation amplitude). Both aspects are of considerable interest in the design stage.

Global dynamics and integrity
The previous investigations highlight that the lowest escape values of forcing amplitude occur at the two main resonance regions, and provide the escape profile as the envelope of forcing amplitude values at which bounded solutions disappear. Yet, such a global stability boundary does not furnish any information about the erosion process of the basins of attraction of the various solutions, which is indeed a critical issue corresponding to system impending jump to contact; therefore, it has no practical utility from the viewpoint of AFM safe operation mode. In this respect, the fundamental concept is the dynamical integrity of the system (see Chapter 2), which depends on the extent of the erosion of its safe basin represented by the region of the phase space bounded by the homoclinic orbit in the underlying unperturbed system of Fig. 2(a).
Tools for investigating the complete basin evolution under a control parameter variation up to escape are the erosion profiles, which allow us to quantify the varying level of basin erosion. Their construction is carried out by means of specific computational tools, based on the safe basin definition and the integrity measure concept. In this work, the safe basin is considered as the union of all classical basins of attraction of the bounded solutions belonging to the system potential well, without taking care of possible transient dynamics out of the safe basin. This means that, in the case of coexistence of more in-well attractors, the safe basin comprises two, or more, competing basins. However, safe basins in the more classical meaning of basins of attraction of given solutions will also be considered when being interested in evaluating the robustness of competing attractors, too, along with the respective roles in the erosion process. The integrity indicators used to build the erosion profiles are the global integrity measure (GIM), which represents the normalized hypervolume (area in 2D) of the safe basin, and the integrity factor (IF), which is the normalized radius of the largest hyper-sphere (circle in 2D) entirely belonging to the safe basin, and represents a measure of the sole compact part of the safe basin. In this last respect, note that IF is used instead of the local integrity measure (LIM) also accounting for the sole compact part, because the safe basin of interest here is the whole potential well. The phase space window x ∈ [−0.3, 0.3], y ∈ [−0.65, 0.25] has been considered in the numerical simulations, since it contains the compact part of the basin of each of the main attractors involved in the erosion/escape process. Focusing around the fundamental resonance, which from what obtained in Sect. 2.2 results to be an important critical region as regards the system stability, Figure 9 displays the evolution for increasing forcing amplitude of the basins of attraction before (ω u = 0.7, Fig. 9(a)), in the close neighborhood (ω u = 0.8 ∼ = ω 1 , Fig. 9(b)) and after (ω u = 0.9, Fig. 9(c)) the resonance, representing different features of the erosion process leading to escape. Before the resonance ( Fig. 9(a)), a weak competing (purple) basin corresponding to the resonant P1H solution appears within the in-well safe basin initially coinciding with the (green) basin of attraction of the nonresonant P1L response, close to its boundary. However, as the amplitude slightly increases, the two basins are separated by erosion tongues of the unbounded solution (white) basin surrounding the well (third panel of Fig. 9(a)), and the small basin of the resonant solution is rapidly eroded up to its complete disappearance. After that, the erosion of the main nonresonant basin processes up to the escape with smooth, i.e. uncorrupted, basin boundary.
When moving closer to the resonance condition, in between nonlinear and linear resonances ( Fig. 9(b)), the competition between nonresonant and resonant basins, whose robustness is now comparable, becomes definitively stronger with the newly born P1H basin which grows up swiftly, and with smooth boundary, within the in-well safe basin (second panel of Fig. 9(b)) up to full replacement of the original P1L basin. For increasing amplitude, it is then raggedly eroded from the unbounded solution basin (fourth panel of Fig. 9(b)).
Finally, after the resonance ( Fig. 9(c)), the erosion scenario exhibits a sequence of competing basins corresponding to higher periodicity solutions, each one of them lasting for a limited forcing amplitude range. Their distributed small sub-basins are located close to the boundary of the sole 1-period basin, a circumstance that entails an overall ragged aspect of the latter along the erosion process. The subsequent onset of the new competing 1-period P1 * solution, seen also in the final part of the bifurcation diagram of Fig. 4, for high values of the forcing amplitude (above U = 0.604), brings to an inside-the-well, secondary, seemingly double-well potential structure of the competing basins (third panel of Fig. 9(c)), which then split from each other with strongly fractal edges in between (fourth panel of Fig. 9(c)).
The richness of erosion patterns underlying escape is also reflected on the erosion profiles obtained by calculating IF and GIM integrity measures from the basins evolution of Fig. 9. In Fig. 10, reference is made both to the individual basins of main periodic solutions (light gray and black curves) and to the total in-well safe basin (thick black lines). GIM and IF profiles are plotted in solid and dashed lines, respectively, and are nondimensionalized with respect to the safe basin of the unforced system (i.e. for U = 0), to obtain dimensionless numbers.
At ω u = 0.7 ( Fig. 10(a)), the safe basin profile exhibits a sharp, though relatively limited, fall down at low amplitude values when the resonant and nonresonant basins separate inside the well (at around U = 0.0067). After that, its evolution, which now coincides with that of the P1L basin, evolves with a long and smooth decrease down to zero. Differently, close to the resonance ( Fig. 10(b)) the strong competition between the basins is highlighted by the crossing between the decreasing P1L (light gray) profile and the increasing (black) one related to the new growing up P1H solution. The internal competition occurs up to nearly complete substitution of the original P1L basin with the new P1H one, which is immediately followed by an ever sharper erosion of the in-well safe basin (thick black, now coinciding with the P1H one), due to the surrounding escape tongues. For frequencies higher than the resonance one, the safe basin profile coincides with that of the main 1-period solution, since the size of the other high period basins is negligible compared to the well. Its evolution develops smoothly for a wide range of forcing amplitude, up to the separation of the two 1-period basins inside the well which causes a sharp fall down to zero coinciding with the loss of safety of the system.
The integrity curves highlight also a different behavior of the two integrity measures, ensuing from their respective definitions and related to the different features of the erosion process. In fact, when the erosion develops smoothly and from the outer edge of the safe basin or, more generally, in case of weak competition between the basins, the volume (GIM) of the basin is affected in a major way than its compact core (IF), thus entailing a lower value (i.e., in design terms, a major conservativeness) of the GIM measure. In contrast, when the competing basins undergo a strong rolled up evolution, or in presence of very marked fractal eroding tongues, the compact part of the basins is reduced more than its entire volume, and the IF measure becomes more conservative (or even much more conservative, if looking at the single attractors' profiles) than the GIM one, from the safety/robustness viewpoint.
It is worth looking at a summary diagram of GIM profiles of the in-well safe basin for different forcing frequencies: Fig. 11 shows an erosion surface in the range ω u = 0.5 − 1.8, with several isointegrity curves obtained by expressing the erosion profiles in terms of remaining safe basin percentage. The profiles have the classical qualitative behavior of the so-called "Dover cliff" erosion curve (Thompson and Stewart, 2002), which is characterized by a slow decrease of the   Figure 11. Isointegrity curves (a, b) and overall erosion surface (b). uneroded volume of the safe compact region, followed by a sudden fall down to zero. Near the two main resonance frequencies, the surface shows two evident depressions, with the lowest percentage values of residual integrity before the V peak and the sharpest decrease of the profiles just after the peak, thus confirming what already noticed about the differences between the profiles of Figs. 10(a) and 10(b).
As a practical comment, it is important to point out that slight discrepancies between the escape threshold of Fig. 11(a) and that of Fig. 3 are due to different numerical procedures used in their realization. Thresholds of Fig. 3 are the result of a continuation analysis, for very small variations of the control parameters, applied to several local bifurcation points, which is able to catch small changes in the curves trend thus furnishing quite accurate and reliable outcomes. Differently, thresholds of Fig. 11(a) are obtained as envelope of the (discretized) integrity values deduced from the construction of the erosion profiles. The latter, which are computationally burdensome to build, have been realized at frequencies steps which are much wider than the previous case, thus furnishing less refined, though reliable, results.
From the viewpoint of AFM safe operation mode, some matters are of considerable interest as regards the evaluation of system theoretical and practical stability. The first one is concerned with the comparison of the overall escape threshold considered up to now and carried out by mapping the bifurcation diagrams (bd), to be considered as a global stability boundary, with the escape threshold obtained by looking at the response under numerical integration (ni) with fixed initial conditions (i.c.), to be considered as a local stability boundary. The comparison is reported in Fig. 12(a), where the dashed threshold, related to numerical simulation, represents the forcing amplitude values U , for several forcing frequencies, at which divergence of a specific system response does occur. Of course, the ni escape values strongly depend on the particular selection of the i.c. pair; here, values corresponding to the equilibrium position of the unperturbed system have been chosen. With respect to all of the local ni boundaries to be possibly constructed by considering divergence from different pairs of i.c. inside the safe basin, the global bd boundary represents the upper bound, since it corresponds to safe basin annihilation. For the considered selection of i.c., the figure highlights significant differences between the two boundaries, the local one lying below the global, with a gap of up to 25 times the ni absolute value. Moreover, ni local minima occur at primary and secondary resonance frequencies, while bd minima are shifted towards left (i.e., to nonlinear resonances) due to the softening behavior of the system. The apparent underestimation of system stability with ni ensues from the particular selection of the i.c. pair against the basin erosion scenario. Focusing on fundamental resonance and looking, e.g., at ω u = 0.8, a slight increase of the excitation amplitude (around the local escape vertex value U = 0.03) is seen to cause a basin boundary erosion which swiftly expands up to including the initial position (black point in Fig. 12(b)), while leaving more than 50% of residual integrity of the safe basin. In contrast, the actual basin annihilation (bd threshold) occurs only at U = 0.74. This clearly highlights how, in terms of overall system safety with respect to escape, consideration of the outcome of a single trajectory may furnish misleading and too conservative information, unless being specifically interested in the response ensuing from that particular set of i.c.
Anyway, even correctly referring to the global stability boundary in terms of overall escape, the major problem in a safety assessment perspective ensues from the associated total lack of information about the features of the underlying basin erosion. Hence, in practical applications, it is particularly important to refer to integrity evaluations in order to determine acceptable frequencydependent thresholds associated with a priori safe design targets. Figure 13 shows four isointegrity curves corresponding to increasing target values in the fundamental resonance region. Selecting for instance the 30% residual safety target, the corresponding threshold allows us to critically discuss the results furnished by numerical simulations (dotted red threshold). Away from the nonlinear resonances (left of them, for the principal resonance see Fig. 12(a)), though the ni thresholds underestimate system safety with respect to the bd threshold, they are definitely unacceptable because of corresponding to very low values of residual integrity (0 − 30%). Yet, the even more questionable point is that ni thresholds correspond to a residual integrity which is strongly variable over the control parameter range (see, e.g., the frequency range in between 0.7 and 0.8 in Fig. 13) -and thus unreliable in overall terms -, even though the associated threshold may become more and more over conservative (thus corresponding to a higher residual integrity) just where this is more needed. The above comparison highlights the importance of a global analysis which Figure 13. Comparison between theoretical (ni (dotted red line) and bd (continuous black line)) and practical (residual integrity, gray lines) boundaries close to fundamental resonance, with detection of some increasingly residual isointegrity curves.
solely accounts for the features of the progressive decrease of the system practical safety due to the erosion up to its final ending to the unwanted escape.
As a final remark, it is worth reporting here that the same analyses have been developed for the system under the sole vertical, harmonic, excitation V g (see Rega and Settimi (2013)). The results concerning the basin erosion process and the local escape thresholds exhibit the same qualitative behavior as the one obtained under parametric excitation, even if the absolute minimum of the total escape threshold is shifted from the subharmonic resonance range to the primary one, consistent with the well-known higher response amplitudes occurring at the latter for an externally driven system.
Furthermore, the influence of combined horizontal scan excitation and vertical beam excitation, which coexist in the noncontact AFM operation mode, on the overall dynamics has been investigated, highlighting a negligible difference of results with those obtained from the sole parametric excitation case, therefore confirming that, for the chosen set of parameters values, the dual excitation analysis is not even necessary in the primary resonance region in which the vertical excitation produces the major effects.

AFM with external feedback control
The results presented in the previous Section have highlighted that slight changes of parameters and/or initial conditions of the AFM model can produce dramatic modifications of its dynamical response, leading to possible unstable oscillations which restrict the operating range of the device. AFMs working in noncontact regime, specifically, can undergo the unwanted "jump-to-contact" phenomenon, or escape in dynamical terms, due to the atomic attraction between the cantilever tip and the sample to be scanned, which can become stronger than the beam restoring elastic force making the equilibrium configuration unstable and producing contacts between tip and sample responsible for errors in the topography process. To avoid these undesirable effects and improve the microscope performances, several control techniques have been proposed and successfully implemented in the last decades, mostly based on feedback control methods (Yamasue and Hikihara, 2006;Arjmand et al., 2008;Payton et al., 2011). Among them, a simple external feedback control can be implemented in the continuum formulation of the noncontact AFM model presented in Sect. 2. Thanks to the control of the cantilever base position (Fig. 14(b)), the procedure works by keeping the system response (which can change due to the ensuing variations in the sample surface position) to a reference one ( Fig. 14(a)), which represents a stable periodic motion of the uncontrolled system under the same setting of operating parameters. Following the model formulation reported in Settimi et al. (2015), the continuum equations for the controlled system read: whereξ(t) is the new control variable representing the distance of the fixed side of the microcantilever from the horizontal reference axis (see Fig. 14(b), which refers to the nondimensional system),v ref (r, t) is the reference vertical displacement, obtained from the uncontrolled system  (ξ = 0, Eq. (1)), andk is a feedback constant. As done for the uncontrolled model, the initialboundary value problem of the distributed parameter system is reduced to a system of two ODEs with one and a half degrees-of-freedom through a Galerkin procedure. After neglecting some terms related to internal feedback control and nonlinear geometry, and by considering the presence of the sole horizontal parametric scan excitation (which is assumed to be harmonic), the nondimensional reduced order system results: where x(t) and z(t) are the tip transverse displacement and control variable, respectively, and k g is the external feedback control parameter. Obviously, the system (9) has to be completed by the addition of Eq. (4) (with V g = 0 and U g = U sin ω u t) relevant to the uncontrolled system, which furnishes the reference response x ref essential for the control application. The control works when the response settles onto the reference one, i.e. when z is equal to the expected value z s , as shown in the sample case of Fig. 15. Here, the change of 3% in the tip-sample distance, represented by the increase of the z s parameter with respect to the reference configuration z s = 0.0 causes a modification in the system dynamical response, which passes from the reference 1-period solution (dashed thick black line in figure) to a high period one (gray thin line in figure). Starting from the latter configuration, the addition of the external feedback control is able to change the system dynamics by avoiding the occurrence of the highly periodic behavior; the controlled model, in fact, manages to properly reproduce the reference 1-period solution (orange thick line in figure), and to set the control variable z to the desired z s value, z = z s = 0.03. The study of the system equilibria points out that they are not influenced by the feedback control, as they result to be the same of the reference uncontrolled system. However, due to the presence of control, the arise of a Hopf bifurcation locus reduces the range of stability of the only acceptable stable equilibrium (the upper one in Fig. 2(b) uncontrolled) also for very low values of the feedback parameter k g . For the parameters choice (7), the asymptotic stability of the equilibrium solution occurs only for 0 < k g < 0.00223 (Settimi et al., 2015).

Local bifurcations and response scenarios
In order to investigate and verify the overall effects of the control introduction on the system behavior, a comprehensive analysis of its dynamical response is carried out and compared with the results already obtained for the uncontrolled system in Sect. 2.2 (Settimi and Rega, 2016c). For the set of parameters values reported in (7), Figure 16 summarizes the system dynamics under variation of the forcing amplitude U , with identification of the different local bifurcation thresholds. The local bifurcation loci which represent the escape thresholds for the underlying uncontrolled system (i.e. saddle-node (gray) and period doubling (black) thresholds in Fig. 16) are present also in the controlled case, but they do not correspond to the overall stability boundaries anymore. The control introduction in the model, in fact, causes the birth of new thresholds of torus (or Neimark-Sacker) bifurcation (TR) and of transcritical bifurcation (T) -the latter always appearing at much higher amplitudes than the former -, which considerably modify the system stability region. In the lower frequencies range up to the fundamental resonance (ω u = 0.835 = ω 1 ), the system escape threshold is governed by the torus bifurcation of the nonresonant 1-period (P1L, i.e., low-amplitude) solution which occurs before the saddle-node governing the onset of escape in the uncontrolled case, while the transcritical bifurcation of the resonant 1-period (P1H, i.e., high-amplitude) solution occurs before the period doubling sequence which instabilizes the response of the uncontrolled case. For higher frequencies which include the principal resonance (ω u = 1.67 = 2ω 1 ), the system shows the coexistence of 1-period and 2-period solutions -the latter governing the nonlinear response in the neighborhood of the 1/2-subharmonic resonance, as for the reference system -, which become unstable via a couple of period doublings and a torus bifurcation substituting the period doubling, respectively. It is worth underlining that the stable region of Fig. 16 refers to system solutions for which the feedback control works properly, i.e., stable 1-period responses of the type (x ref , z s ), such as the nonresonant one (P1L) obtained for the forcing frequency ω u = 0.8 before the torus bifurcation (TR) and the P1H solution delimited by two transcritical bifurcations (T points) marked in the sample bifurcation diagram of Fig. 17. However, the latter becomes stable only after the system has experienced a region of instability characterized by unbounded responses (white tongues inside the stable region in Fig. 16), a fact that makes it unsafe with respect to possible variations of the operational parameters. According with the control efficiency, stable motions exhibited by the system that do not settle onto the reference one are considered as unwanted outcomes and thus out of the stability region. This is the case, for example, of the new 1-period solution P1 born from the transcritical bifurcations T (represented with green lines in Figs. 17 and 18), which is a high-amplitude response coexisting with the main ones and displaying a rich scenario of local bifurcations. Indeed, bifurcation diagrams of the control variable z show that for these additional solutions the system fails to reach the reference position (i.e., z = z s = 0.01, see the right panel of Fig. 18(a)), actually making the z response periodic, so that they appear as unwanted responses. The reported bifurcation diagram highlights also the presence of stable quasi-periodic solutions arisen from the torus bifurcations (Fig. 17), which can be numerically detected and for which the control variable z does not actually reach the expected position z s = 0.01 (see Poincaré map of the z variable of Fig. 18(b)). Note that the quasi-periodic solution exists up to the occurrence of the saddle-node bifurcation that makes the reference solution disappear; since it represents the input for the numerical solution of the controlled system, this local bifurcation marks also the death of the stable quasi-periodic response.
Figs. 16 and 17 are obtained for low values of the feedback control k g ; indeed, it is interesting to investigate the influence of the increasing feedback control parameter on the stability of the main periodic solutions, since it is known that changes in such parameter can strongly modify the system response. Obviously, it must be taken into account that its maximum acceptable value corresponds to the one instabilizing the system equilibrium, i.e. k g = 0.00223. Thus, focusing around fundamental resonance, several behavior charts for increasing k g are reported in Fig. 19, to be compared with the one obtained for the uncontrolled system (i.e. k g = 0, Fig. 19(a)). Figure  19(b) shows the appearance of an unstable tongue delimited by the new torus and transcritical bifurcation thresholds which occur for low values of the forcing amplitude U , and even for weakly controlled systems. With the control introduction, in fact, the stability of the nonresonant P1L solution is reduced by the torus threshold, while the resonant P1H response loses its stability for low values of the forcing amplitude U , returning stable only after the transcritical bifurcation (Fig. 17). As a consequence, the typical triangle region below the V vertex displayed by the uncontrolled system with the coexistence of stable P1L and P1H solutions is not only reduced by the torus curve, but displays the presence of the sole nonresonant solution, up to its complete disappearance for higher k g values. For higher values of the feedback control parameter, (Fig. 19(c)), an unstable region confined by a torus bifurcation (TRH) arises below the upper transcritical threshold, and expands as the control increases. Simultaneously, the torus bifurcation which unstabilizes the P1L solution occurs for decreasing values of U . After that, stability regions reduce to narrow strips of existence of stable P1H solution, associated with very limited ranges of forcing amplitude U (see behavior chart at k g = 0.01 of Fig. 19(d)). All over Fig. 19, note that the frequency value ω u = 0.8358 marks the disappearance of the P1L solution, hence for higher frequencies the P1H solution becomes the sole 1-period solution for the system and its "high-amplitude" connotation becomes unnecessary (as for the reference system, see Fig. 6 and relevant comments); nevertheless, for the sake of readability of the charts its label has been kept unchanged.
To summarize the results concerning the controlled system, and in view of easily comparing them with the outcomes of the uncontrolled one (Sect. 2.2), the behavior chart in the frequencyamplitude plane is shown in Fig. 20, where the overall escape thresholds of the two systems (controlled and uncontrolled) together with the detection of the corresponding stability regions are reported. The most evident effect concerns the reduction of the dynamical stability of the system entailed by the external feedback control, mostly around the main resonance frequencies, i.e. ω 1 , 2ω 1 and ω 1 /2.  Figure 19. Behavior charts in the ω u -U plane around fundamental resonance for k g = 0 (a), k g = 0.001 (b), k g = 0.002 (c) and k g = 0.01 (d).
In the close neighborhood of these values, in fact, the feedback control causes the onset of instability tongues which dramatically decrease the escape value of the forcing amplitude, with reductions of about 99,9% around the fundamental (primary) resonance, of 99,4% around the principal (subharmonic) resonance and of about 94% around the superharmonic resonance frequency. This effect can be explained with the fact that the close proximity to the resonance frequencies leads to a substantial increase of the response amplitude of the resonant periodic solutions, that the feedback control is unable to dominate. In these regions, therefore, the escape threshold of the controlled system is governed by the nonresonant responses, which become unstable for considerably lower values of the forcing amplitude. Furthermore, the escape threshold minima of the controlled system are shifted towards frequency values related to the system natural frequency ω 1 (i.e. ω 1 /2, ω 1 , 2ω 1 ), since the nonresonant periodic solutions are not affected by the softening effect of the nonlinear resonance (as, in contrast, happens for the resonant solutions which govern the escape profile of the uncontrolled system).
It is important to remember also that for the controlled system the stability region does not include all the stable periodic responses, but only those associated with the reference periodic solution of the uncontrolled system, on which the control works properly. This causes an additional reduction of the limiting values of forcing amplitude. In terms of extent, this is evident particularly Figure 20. Behavior chart in the ω u -U plane with detection of the overall escape thresholds for the controlled (k g = 0.001, orange line) and uncontrolled (k g = 0, black line) systems under parametric excitation. Dark gray area represents the stability region of the controlled and uncontrolled systems, light gray area represents the stability region of the sole uncontrolled system. around the 2ω 1 frequency, where the escape boundary of the uncontrolled system is governed by the 2-period response representing the main reference periodic solution for the whole range of high frequency values which, however, the control fails to correctly reproduce.
It is worth noting that the above-mentioned response features are equally present both in the parametrically forced system herein summarized and in the externally forced one which could be considered, in the uncontrolled (see Eq. (4)) or controlled case, under vertical excitation of the microcantilever support. Furthermore, the qualitative behavior of the escape threshold is not modified by the control activation, and therefore maintains the typical features of systems under parametric (external) excitations, with V-shaped profiles around resonance regions and the absolute minimum at principal (fundamental) resonance, thus confirming that the applied control acts on the system irrespective of the forcing type. From a practical point of view, the highlighted behaviors have detrimental implications for AFM operation, as they entail limiting the range of both the forcing amplitude and the feedback control parameter to low values to ensure the effectiveness of control. Moreover, they suggest to select the low amplitude solution as the reference one, since the control of the high amplitude solution is effective only in limited regions of the parameter space, so that it would imply a very careful choice of the operational settings.
Finally, the detrimental effect of control on the system stability can be highlighted also when investigating the dynamical behavior as a function of system intrinsic parameters, e.g., the atomic interaction parameter Γ 1 governing the order -2 term of the equation of motion and representing the most significant system nonlinearity. The results, reported in Fig. 21(a) to be compared with those of Fig. 8, show that the addition of control into the model (k g = 0.001, Fig. 21), also in the Γ 1 -U plane, as already seen in the ω u -U plane, causes the disappearance of the region of coexistence of the P1 solutions with the lowering of the escape boundary of the P1L response due to the birth of a new torus TRL threshold. However, the control presence affects mostly the P1H solution, which in the uncontrolled system was confined by the period doubling threshold PDH; the onset of a new transcritical locus (TH) and of a new torus one (TRH) significantly reduces its stability range which passes from Γ 1 ∈ [0.127, 0.149] to Γ 1 ∈ [0.1305, 0.1313] at a forcing amplitude of U = 0.0005. Notwithstanding the small value of the feedback control gain (k g = 0.001), this prevents from operating the AFM for atomic interaction values (e.g., tip material choice, material of the sample to be possibly scanned) slightly higher than the reference one (Γ 1 = 0.1, see the  Figure 21. Behavior chart in the Γ 1 -U plane for ω u = 0.7 and k g = 0.001 (a); time histories and trajectories in the phase plane of the uncontrolled (reference) system (black line) and of the controlled system (orange line) at Γ 1 = 0.12 and U = 0.12 (b). comparison between the uncontrolled and controlled responses of Fig. 21(b) for Γ 1 = 0.12), also with nearly vanishing excitation amplitudes, thus furnishing further practical information of design interest as regards proper AFM operation settings.

Global dynamics and integrity
As previously stated, the investigation of the dynamical integrity of a system when varying some control parameter is strictly related to the definition and choice of specific concepts and tools, namely the safe basin and the integrity measure. From what said in the former Sect. 3.1, it follows that the safe basin for the AFM system with external feedback control must include the sole solutions for which the control procedure works properly, namely the nonresonant (i.e. lowamplitude) 1-period solution around the fundamental frequency. If compared with the reference system (Sect. 2.3), in which the safe basin coincides with the potential well, the system safe basin results to be strongly reduced by the control presence into the model, as it will be described in detail hereinafter.
Notwithstanding the considered minimal (single-mode) approximation of the microcantilever dynamics, the controlled system is five-dimensional in the state space, as it is defined by Eqs. (9) and (4), so that understanding and representing the relevant attractor-basin portraits becomes rather difficult. Therefore, to obtain readable and interpretable results, analyses are accomplished by systematically assuming that the initial conditions of the displacement variables of the controlled system are the same of the uncontrolled one, i.e. x(0) = x ref (0), y(0) = y ref (0) -which also corresponds to considering situations for which the control is more likely to work -, thus reducing the basin to three dimensions. Based on this, the erosion evolution is investigated constructing 2D cross-sections of the five-dimensional basins of attraction with fixed initial conditions of the control variable z.
In a comparison perspective, the tools selected to build the erosion profiles are the same already used to analyze the dynamical integrity of the uncontrolled system, i.e. the global integrity measure (GIM) and the integrity factor (IF) as integrity indicators, and x ∈ [−0.3, 0.3], y ∈ [−0.65, 0.25] as phase space window. Also in this case, normalization has been performed with respect to the integrity measure calculated for the reference safe basin of the unforced system (i.e. for U = 0), so that GIM and IF are dimensionless numbers.
The evolution of the system safe basin is analyzed for increasing forcing amplitude U , with fixed control parameter z s (equal to 0.01), which is related to the distance between microcantilever tip and sample surface, and with fixed initial condition of the control variable z(0) = 0. Several cross-sections of basins of attraction in the (x = x ref , y = y ref ) plane are built for different forcing frequencies around the fundamental resonance (ω u = ω 1 ) region, and the results are compared with those already obtained for the uncontrolled system presented in Sect. 2.3 to evaluate the influence of control on the dynamical safety of the system (Settimi and Rega, 2016a). As in Fig. 9, in Fig. 22 the basins of attraction of the controlled system in the (x = x ref , y = y ref , z = 0) plane are reported for increasing values of the forcing amplitude and for three forcing frequencies exemplifying the system behavior around the fundamental resonance, i.e. ω u = 0.7, ω u ∼ = ω 1 = 0.8 and ω u = 0.9. The projection of the relevant attractors is marked and the z value at which each attractor settles is indicated. To better frame the following global results refer also to Fig. 19(b).
The most evident consequence of the control presence is the generalized increasing of the erosion due to the enlargement of the basin of the unbounded solution (white basin in Fig. 22), which develops at the main detriment of the non-dominant responses. In particular, with respect to the uncontrolled case ( Fig. 9) for frequencies lower than the resonant one ( Fig. 22(a)) the basin of the resonant solution disappears, as well as the basins of the high period solutions for ω u values higher than ω 1 (Fig. 22(c)). At this frequency, moreover, the feedback control is able to properly reproduce the 1-period solution only in a confined range of the forcing amplitude, after which the periodic responses depicted by the system cannot be considered acceptable anymore. This is highlighted by the z values reported in each basin section of Fig. 22(c), indicating the z planes which the attractors belong to; after the critical value U = 0.598 (corresponding to a transcritical bifurcation) the stable solutions do not coincide with the reference ones (i.e. z = z s = 0.01), confirming the bad control operation. Near the fundamental resonance ( Fig. 22(b)), the control causes the onset of tongues of the unbounded solution basin (white) inside the potential well, whose penetration causes the separation of the resonant/nonresonant basins for low values of the forcing amplitude, and whose development occurs to the detriment of the nonresonant basin, which is strongly and swiftly reduced. Such behavior is peculiar of the system with feedback control and modifies the topological scenario of the basins erosion in the close proximity of the resonance frequency; for the controlled system, in fact, the disappearance of the nonresonant basin is caused by the erosion produced by the unbounded basin, while in the reference case ( Fig. 9(b)) it occurs by replacement with the resonant basin later on eroded by the unbounded one, without penetration of the latter, which enlarges starting from the outer edge of the potential well.
Remembering the bifurcation diagrams of Fig. 17, which indeed well samples the general behavior of the system, it is evident that the basin of attraction of the resonant solution corresponds to different solutions if referred to the controlled system or to the uncontrolled one. In particular, in the controlled case it represents the new P1' solution basin for low values of the forcing amplitude U (less than 0.1, i.e. for all red basins in Fig. 22(b), at ω u = 0.8), while it corresponds to the same P1H solution of the reference system, with a proper operation of the control technique, for high values of U where, however, the basin has been already meaningfully eroded.  This can be understood by referring to the 3D (x = x ref , y = y ref , z) phase space; here, the basins reported in Fig. 22(b) correspond to (x = x ref , y = y ref ) sections of the five-dimensional basin at z = 0. While for the nonresonant basin the attractor belongs to a parallel section at z = z s = 0.01, as expected, the attractor of the resonant basin fails to reach the same section, and settles onto a plane with different (and often negative, see Fig. 22) z, as shown in Fig. 23, which refers to the case U = 0.01. Note that the z = z s = 0.01 plane in the figure corresponds to the phase-plane of the uncontrolled reference system, so that the cross section of the 5D controlled basin of attraction in this plane coincides with the second picture of Fig. 9 To quantify the erosion process of the system safe basin illustrated in Fig. 22, the relevant erosion profiles obtained by means of GIM and IF measures, are reported in Fig. 24 and compared to those of the uncontrolled system (black curves). The reduction of the safe basin to the sole 1-period dominant basin causes a lowering of the profiles in the whole range of forcing amplitude, and for all the frequencies analyzed, however with meaningful differences. At ω u = 0.7 (orange profiles of Fig. 24(a)), the erosion develops smoothly as in the reference case, with not even the initial sudden fall due to the escape tongue therein separating the outer resonant basin from the nonresonant one. This feature is worthily kept, at this frequency, also at higher excitation amplitudes, as highlighted by the nearly constant intervals in-between the successive (orange) iso-GIM curves of Fig. 25. In this sense, despite its slightly lower integrity level, the system with feedback control turns out to be more robust than the uncontrolled one with respect to possible perturbations. Looking at the behavior at ω = 0.9, (orange profiles of Fig. 24(c)), on the contrary, the limit value U = 0.598, responsible for the birth of improper new solutions, causes the collapse of the profiles from 50% to 0% (consistent with the strong packing of the corresponding iso-GIM curves of Fig. 25), making the feedback controlled system particularly dangerous when working around these high amplitude values. Around the resonance, conversely, due to the fact that the safe basin of the controlled system is reduced to the sole, swiftly reduced, nonresonant one, the fall down of the erosion profiles relevant to the system with feedback control occurs for very low values of the forcing amplitude, and the controlled system results to be much more dangerous than the uncontrolled one. A further contribution to the robustness reduction is furnished by the way the erosion develops as the amplitude increases. In fact, it involves mostly the formerly rolled part of the nonresonant basin, which becomes very fractal, thus causing a sharp decrease of the GIM profile, while affecting the basin compact part (quantified by the IF measure) in a minor way. Figure 25. Residual iso-integrity curves of the system safe basin for the controlled (z = 0, z s = 0.01, orange solid lines) and reference (z = z s = 0.01, black dashed lines) systems.
As a summary of the obtained results, the iso-integrity curves in the ω u − U plane are shown in Fig. 25, i.e. frequency-dependent thresholds with constant residual integrity which actually govern the system practical safety, to be compared with those obtained for the system in the absence of control.
The overall outcomes highlight that around the resonance frequency the controlled system undergoes a severe worsening of its practical stability, with the residual integrity being reduced from 90% to 10% in a very narrow range of forcing amplitude values (∆U = 3.6·10 −3 at ω u = 0.82 ≈ ω 1 ). Moreover, analyzing the evolution of the iso-integrity curves a shift in their lowest peak can be recognized, with the minimum value moving right from the nonlinear resonance frequency to the natural one, as already observed for the escape threshold of Fig. 20. Both effects are associated with the system softening behavior, that prevents the response from attaining the higher amplitude resonant solution, which is characterized by too low values of residual integrity. This is similar to what observed also in other softening systems, e.g. the MEMS capacitor (Alsaleem et al., 2010), where it has been shown also responsible for the occurrence of practical escape (therein corresponding to pull-in phenomenon) at much lower excitation amplitudes and higher forcing frequencies in experimental conditions than in theoretical conditions (see also Rega and Lenci (2015)).
Note that the loss of stability due to the control presence is narrower in the left side of the chart and wider in the right, due to the fact that the P1H solution, governing the system response in the latter range, is no longer acceptable for the system.
Until now, the attractors robustness and the basin erosion of the system with external feedback control has been investigated with respect to variation of external action parameters (i.e. horizontal scan excitation and frequency). Yet, it is also important to verify the practical integrity of the system when varying some intrinsic parameter, like the distance between microcantilever tip and sample surface represented by the control parameter z s . Due to its dependence on the roughness of the sample to be scanned, this parameter is highly variable during the AFM scanning operation, and the analysis of its effect on the global behavior of the system is particularly important to assess the actual safety in operating conditions. For this reason, several bidimensional sections of the system basins of attraction in the (x = x ref , y = y ref ) plane are constructed for increasing values of z s , fixing the initial condition of Figure 26. Residual iso-integrity curves for the controlled system, at U = 0.01. the control variable z = 0 and the forcing amplitude at U = 0.01, which is indeed a low value, and for different forcing frequencies near the fundamental resonance. Note that varying z s in the theoretical model corresponds to implicitly accounting for the actual time-varying distance during the scan operation.
As a remark, it can be noted that in this case the phase-plane corresponding to the reference uncontrolled system is set to z = z s = 0, differently from what done in the previous Figs. 22-25 where the uncontrolled plane is defined as z = z s = 0.01. Since the reference system is obtained from the controlled one when equations (9) and (4) decouple, i.e. when z = z s for each value of z s , the different choice between the two integrity analyses is due to reasons of consistency. Indeed, while previously z s was a fixed parameter, equal to z s = 0.01, here its variation is investigated, so that it is natural to determine the reference plane as the one for which z s equals the fixed value of the variable z initial condition, i.e. z s = z = 0.
The outcomes highlight that when the tip-sample distance increases (i.e. for increasing z s ) the basin of the unbounded solution enlarges inside the potential well, to the main detriment of the (solely controllable) nonresonant solution, whose basin of attraction is dramatically reduced up to vanishing. The negative effect of the increasing tip-sample gap on the dynamical safety of the system is summarized by the iso-integrity curves reported in Fig. 26, where the resonance region proves to be critical also with respect to variation of the tip-sample distance. Around ω 1 , in fact, the dynamical integrity of the controlled solution is drastically reduced as the gap increases, decreasing up to ≈ 10% for z s = 0.01.
Finally, it is of interest to analyze the combined effect of the two main system intrinsic parameters, i.e. the z s parameter and the atomic interaction parameter represented by Γ 1 .
Analyzing and quantifying the basin evolution as a function of z s at different values of the atomic interaction parameter provides the erosion profiles of the safe basin in Fig. 27. They confirm that an increase of either parameters has an evident negative effect on the system robustness. This information has an important practical consequence, as indicates that a rough sample surface and/or a strong atomic interaction between tip and sample represent dangerous situations for the application of the external feedback control to an AFM, which could lead to errors in the topography. Moreover, for a given acceptable value of residual integrity, the z s value available for control in safe operating conditions becomes smaller as the atomic interaction increases. Overall, in a design perspective, useful hints are obtained to calibrate the tip-sample interaction (e.g., tip material choice) depending on the sample characteristics and roughness, in order to guarantee a proper safety level during the scan operation.

AFM with global control
The results presented in the previous Section 3 have shown that, while being generally effective in realizing the specific aim for which it is designed, the insertion of an external feedback control in a noncontact AFM model causes a generalized reduction of the stability region and a dangerous decrease of system safety with respect to the unwanted jump to contact phenomenon. Yet, from a practical stability perspective, an acceptable system-dependent residual integrity is needed to guarantee secure AFM operation since it is nowadays well known that the safety of a nonlinear system depends not only on stability of its solutions but also on the uncorrupted basins of attraction surrounding them.
Following this assumption, and focusing on the preservation of dynamical integrity, it is worth analyzing whether a non-feedback control technique, specifically aimed at favorably affecting a global bifurcation event which triggers basins erosion, can also work for the AFM system and reduce the basin erosion which leads to the loss of safety. The method (Lenci and Rega, 2004) is based on the optimal modification of the shape of the reference harmonic excitation in order to shift the occurrence of the global event (i.e., homo/heteroclinic bifurcation of some saddle) responsible for the sudden fall down of the dynamical integrity profiles, thus obtaining an overall control of the dynamics and an enlargement of the system safe region in parameters space. When the global bifurcation triggering erosion involves the manifolds of the sole hilltop saddle of the system, then its (approximate) occurrence can be determined by means of the Melnikov method, based on the interpretation of the nonlinear system as perturbation of the relevant Hamiltonian one, and the subsequent determination of the amplitude and phase of the controlling superharmonics by solving an optimization problem. This analytical control procedure has demonstrated its effectiveness in increasing the system safety of several different mechanical systems, such as the Helmholtz oscillator (Lenci and Rega, 2003a), the Duffing oscillator (Lenci and Rega, 2003b), the rigid block (Lenci and Rega, 2005) and a reduced model of MEMS microbeam (Lenci and Rega, 2006), which have allowed to verify the influence on the erosion of some main mechanical features, like smoothness, stiffness and symmetry. With reference to the latter, two different approaches have been proposed, the "one-side" control, which is very effective but acts on a single global (homoclinic) bifurcation thus being particularly useful for asymmetric systems, and the "global" control, which is able to increase simultaneously two homo/heteroclinic bifurcations in case of two (or more) potential wells or heteroclinic loop involving two saddles. In the examined cases, the analytical procedure has demonstrated its ability in delaying the erosion also in cases where the loss of safety is likely due to the presence of secondary global bifurcations involving other internal saddles. However, when other secondary saddles play a major role in the development of the basins erosion, and when the profile fall down occurs away from the hilltop bifurcation, the control technique has to be applied by developing a purely numerical procedure to define the optimal shape of the controlling superharmonics.

Analytical control of homoclinic bifurcation of hilltop saddle
In order to apply the analytical procedure, reference has to be made to the underlying Hamiltonian system (see Eq. (6) and Fig. 2 of Sect. 2.1), from which the expression of the homoclinic orbit y h (t) of the hilltop saddle S H , reported also in Fig. 28, can be derived (Settimi et al., 2016): where x S H is the coordinate of the hilltop saddle. The integral (10) cannot be solved in closed form but has to be computed numerically. The curve borders the single potential well, which includes the bounded periodic solutions and represents the safe domain for noncontact AFM operation, separating it from the outer region of unbounded solutions leading to the jump to contact phenomenon. Due to the unperturbed nature of the Hamiltonian system, the homoclinic orbit represents also the coinciding right unstable (W u (S H )) and stable (W s (S H )) manifolds of the hilltop saddle. When damping and excitation are added into the system (6), the separation of the stable and unstable manifolds of the Hamiltonian system occurs (see Fig. 28, where the (gray) homoclinic orbit is also plotted, for comparison), and their evolution can lead them to intersect or to remain disjoint depending on the value of the excitation parameters. The critical situation of manifolds tangency corresponds to the occurrence of a global homoclinic bifurcation, which can be analytically detected by the classical Melnikov method. The control procedure aims at shifting such bifurcation by optimally modifying the excitation shape, i.e. by adding controlling superharmonics to the harmonic excitation of the system (4) with V g = 0 and U g = U sin ω u t: Figure 28. Qualitative behavior of the stable (W S ) and unstable (W U ) manifolds of the hilltop saddle for the Hamiltonian system (gray) and the ε-perturbed, forced, damped system (black).
where U 1 is the horizontal parametric-like reference excitation while U j and Ψ j are the amplitudes and phases of the controlling superharmonics. In presence of weak excitation and damping, the system (11) can be expressed as an ε-perturbation of the Hamiltonian system (6): where ε is a smallness parameter. To the first order, the distance between the stable (W s (S H )) and unstable (W u (S H )) manifolds of the perturbed system is furnished by the Melnikov integral: where are integrals to be numerically computed, and is 2π-periodic with zero mean value (m = ω u t 0 ). The tangency of the stable and unstable manifolds, i.e. the occurrence of the homoclinic bifurcation, corresponds to a simple zero of the Melnikov integral M(t 0 ) = 0 for some t 0 : Note that in presence of the sole reference harmonic excitation (j = 1) h 1 = 1 and M = 1. To achieve the best control effect, the homoclinic tangency has to be shifted to the highest possible value of the forcing amplitude, which corresponds to the following optimization problem (Lenci and Rega, 2003a): Maximize min The optimization can be achieved by properly choosing h j and Ψ j , where j = 1, 2, 3, .., N is the number of superharmonics to be added to the system. To measure the improvement obtainable with respect to the reference harmonic excitation, the gain G is introduced: G = U 1,cr /U h 1,cr = 1/M , where U h 1,cr represents the critical amplitude for the reference harmonic system. Setting Ψ j = 0, focus is paid to the addition of one even superharmonic, i.e. j = 2, due to the asymmetric feature of the system which calls for the application of the so called "one-side" control. The solution of the optimization problem (15) furnishes the value h 2 = 0.353553, as its minimum reaches the maximum value in h(m) = −0.707107, corresponding to the maximal theoretical gain G = 1.4142 (Lenci and Rega, 2004). The homoclinic bifurcation thresholds detected by means of the Melnikov method are reported in Fig. 29 in the excitation parameters plane (ω u , U 1 ), for the reference harmonically forced system (black curve), and for the system with one optimal controlling superharmonic (green curve). The outcomes confirm that the control manages to shift the bifurcation threshold to higher values of the forcing amplitude with respect to the harmonic case, apart from the region around the "anti-resonant" frequency ω u ≈ 1.03 at which the threshold goes to infinity due to the first order nature of the Melnikov approximation. Figure 29. Analytical bifurcation thresholds in the ω u -U 1 plane. Black curve: reference harmonic system; Green thick curve: system with superharmonic h 2 = 0.353553 (optimal).
The numerical detection of the hilltop manifolds validates the analytical results, highlighting the ability of the optimal superharmonic in separating the stable (gray) and unstable (green) hilltop manifolds, as shown in Fig. 30 at frequency values before and after the fundamental resonance ω 1 = 0.8357. It is worth noting that the increase of the bifurcation threshold occurs for negative values of the amplitude of the controlling superharmonic (see Fig. 34 forward), while positive amplitudes produce a worsening effect of the global dynamical behavior. From an operational point of view, the same result can be achieved by applying positive superharmonics with a phase shift Ψ 2 = π.
With reference to the main practical objective of the applied control, i.e. the overall enlargement of the system safe region, the verification of the homoclinic hilltop bifurcation shift has to be complemented by the investigation of the safe basin evolution, the latter being carried out by the analysis of the system dynamical integrity. In this respect, the safe basin is defined as the union of the classical basins of attraction inside the potential well, and the integrity measures used to quantify the erosion are the Global Integrity Measure (GIM) and the Integrity Factor (IF) which work both well in this respect. Once the basins of attraction are numerically obtained for increasing values of the forcing amplitude, and after quantification of the basins integrity, the erosion profiles of Fig. 31 are built, for the system with harmonic excitation (black curves) and for the one with optimal superharmonic (green curves). The results highlight that the controlling superharmonic is not able to improve the system integrity, i.e to shift the profiles fall down to higher values of the forcing amplitude, different from what observed in several other systems with escape (Rega and Lenci, 2008). Moreover, right of the fundamental resonance the controlled profiles result to be sharper than the harmonic ones in their final part (see Fig. 31(b)), meaning that for high values of the forcing amplitude the control worsens the system robustness (this is anyway a well known behavior, detected in other literature systems (Rega and Lenci, 2008)). The causes of this behavior are to be sought in the significant distance, in terms of values of the forcing amplitude, between the occurrence of the homoclinic bifurcation and the profile fall down, as can be seen in Fig. 31. Focusing on the case ω u = 0.7 and observing the basins evolution (not reported), in fact, it can be noted that, once the manifolds intersect, the unbounded solution basin which surrounds the potential well starts penetrating inside the well through fractal tongues from the outer edge. However, the safe basin erosion develops very slowly and, when the amplitude approaches the value U 1 = 0.0067 (numerically detected) corresponding to the profiles fall down, the positive effects of the homoclinic bifurcation shift have been exhausted; in fact, at this value the stable and unstable Figure 30. Stable (gray) and unstable (green) manifolds of the hilltop saddle S H for the reference harmonic system (a,c) and for the system with optimal superharmonic h 2 = 0.353553 (b,d) at U h 1,cr = 0.001823 and ω u = 0.7 (a,b) and at U h 1,cr = 0.0027652 and ω u = 0.9 (c,d). Square point identifies one of the manifolds tangency points. manifolds of the hilltop saddle clearly intersect each other. Of course, the same considerations hold at frequencies higher than the fundamental resonance one, as well as, for both cases, the arise of competing basins of attraction before the profiles fall down. Their appearance suggests the presence of other internal saddles whose manifolds are responsible for the basins behavior inside the well, which, however, has to be investigated with solely numerical procedures.

Numerical control of bifurcations of secondary saddles
Before the fundamental resonance frequency An accurate numerical investigation of the global bifurcation scenario of the AFM model at ω u = 0.7 is required to properly identify the global event responsible for the reduction of system integrity. The results reported in Table 4.2 display the occurrence of four different global bifurcations, involving stable and unstable manifolds of a secondary saddle internal to the potential well, in addition to the hilltop manifolds. In particular, the arise of such secondary saddle S 1 is related to the appearance, at U = U S1 = 0.001697, of the basin of attraction of the 1-period (resonant) solution inside the former (nonresonant) basin, and its stable manifolds represent the smooth edges in between the two basins. This event occurs just before the homoclinic bifurcation of the hilltop saddle at U = U hom1 = 0.001823, already detected and controlled in the previous Sect. 4.1, and corresponding, from a phenomenological viewpoint, to the beginning of the fractalization of the safe basin outer edge. At considerably higher values of the forcing amplitude, at U = U hom2 = 0.006375, the tangency between the unstable left manifold W u l (S 1 ) and the stable left manifold W s l (S 1 ) of the in-well saddle points out the occurrence of the left-side homoclinic bifurcation of the saddle S 1 , which marks the beginning of the disruption of the smooth separation between the basins inside the well, representing the starting point for the (a) (b) Figure 31. Erosion profiles as a function of the forcing amplitude U 1 for the reference harmonic system (black) and for the system with optimal superharmonic (thick green), at ω u = 0.7 (a) and at ω u = 0.9 (b). Dashed lines: profiles obtained with Integrity Factor (IF); Continuous lines: profiles obtained with Global Integrity Measure (GIM).
fractalization of the resonant basin boundary inside the nonresonant one. As the amplitude slightly increases, the left unstable manifold of the in-well saddle crosses the left stable one and approaches the right stable manifolds of both the in-well saddle (purple line in the following Fig. 32) and the hilltop one (gray line), whose evolutions have developed keeping them very close to each other. At U = 0.006676, the two mixed-side tangencies (W u l (S 1 ) ∩ W s r (S 1 ) and W u l (S 1 ) ∩ W s r (S H )) occur almost simultaneously, corresponding to a homoclinic bifurcation of the in-well saddle and to a heteroclinic bifurcation between the in-well and hilltop saddles, respectively, reported in Fig. 32. They lead to the penetration of tongues of the unbounded solution basin inside the potential well, which causes the separation of the resonant and nonresonant basins (right panels of Fig. 32), and consequently the strong reduction of the safe basin integrity highlighted by the sharp fall down of the erosion profile. The described global bifurcations are marked on the erosion profiles of the system safe basin in Fig. 33, confirming that the two last bifurcations occur just before the profile sharp decrease, thus suggesting they are the global events to be controlled for increasing the system safety.
To this aim, a fully numerical procedure is developed (Settimi and Rega, 2016b), starting from the numerical detection of the saddles and of the relevant stable and unstable manifolds involved in the bifurcation to be controlled. Then, a proper region in the state plane including candidate, and graphically handy, bifurcation points is identified, and the computation of the manifolds distance along a proper (i.e., non parallel to one of the manifolds) direction is carried out via the arclength method. Once one manifold has been numerically discretized, distances are computed starting from each point of it and the minimum value is selected. The direction for calculating the distance is chosen to be that of the unstable eigenvalue of the hilltop saddle, although distances along other (rotated) directions have been also calculated to check the consistency of the obtained Table 1. Main global events at ω u = 0.7.

Harmonic Amplitude
Global Event Saddle involved  results. Apart from obvious quantitative differences, the outcomes have shown the same trend, thus validating the adopted choice. The procedure is repeatedly implemented by varying the amplitude of the harmonic excitation U 1 , until the calculated distance becomes equal to zero, warning that the tangency of the manifolds (thus the global bifurcation) occurs. The application of the same process is then carried out for different amplitudes of the controlling superharmonic, i.e. various U 2 /U 1 ratios in Eq. (11), assuming Ψ j = 0 and j = 2, as done in the analytical case illustrated in Sect. 4.1, thus considering the presence of one even superharmonic (2ω u ). The subsequent detection of the bifurcation threshold in the U 2 /U 1 , U 1,cr plane allows to finally determine the optimal superharmonic able to shift the bifurcation to the highest value of the forcing amplitude.
To test the numerical control, the homoclinic bifurcation of the hilltop saddle occurring at U 1 = U hom1 = 0.001823 has been investigated, since for it the Melnikov method already allowed to identify the value of the optimal controlling superharmonic (Sect. 4.1). The good agreement between the analytical and numerical results shown in Fig. 34(a) confirms the ability of the numerical procedure in detecting the value of the superharmonic to be added to the system for maximizing the delay of the global bifurcation.
Thus, the technique is applied to the bifurcations occurring just before the erosion profile fall down, focusing the attention on the hom3 bifurcation, since it requires the detection of the position and of the relevant manifolds of the sole internal saddle, simplifying the implementation of the numerical procedure. The selected region in the state plane is chosen to be x ∈ [−0.2, 0.2], x ∈ [−0.3, 0.3], since it includes a significant number of tangency points. The obtained bifurcation threshold for varying amplitude of the controlling superharmonic is reported in Fig. 34(b), and allows to determine the new value of the optimal superharmonic, which corresponds to the peak of the threshold occurring at U 2 /U 1 = −0.32. It is worth noting that, in this case, the whole bell of increased bifurcation thresholds is contained in a range of amplitude values of the added superharmonic which are meaningfully lower than the value of the reference harmonic amplitude, thus corresponding to a control suitability also in terms of relatively low cost.
Once the optimal amplitude of the controlling superharmonic is detected, the erosion profiles of the controlled system are realized and compared with those of the uncontrolled one, to verify the effectiveness of the control in delaying the sharp decrease of the dynamical integrity. Figure 35 confirms that, beyond delaying the occurrence of the targeted homoclinic bifurcation, the numerical control is able also to move the profiles fall down to higher values of the forcing amplitude, in this case increasing the system safety of about 4%. This beneficial effect can be easily justified by observing the influence of the control on the behavior of the basins of attraction. With reference to the uncontrolled system (Figs. 32(d,f)), the optimal superharmonic manages to avoid the manifolds intersection (see comparison between Figs. 32(d,e) without control and Figs. 36(a,b) with optimal superharmonic), and entails the reduction of the basins fractalization and the postponement of the penetration of the white tongues of the unbounded basin inside the potential well, as highlighted comparing Fig. 32(f) and Fig. 36(c). As a consequence, the safe basin, i.e. the union of the resonant Figure 35. Erosion profiles for the uncontrolled system (black line) and for the system with optimal superharmonic U 2 /U 1 = −0.32 (thick green line). and nonresonant basins taken into account by the IF measure, maintains its compactness for a wider range of forcing amplitudes.
With reference to the outcomes presented in literature concerning different dynamical systems, it is worth observing that the presented analysis furnishes coherent and acceptable results. In fact, for the AFM system investigated here, the improvement of about 4% of the dynamical robustness is obtained with relatively low cost in terms of controlling amplitude (U 2 /U 1 = 0.32). Other literature systems are in line with these values: for an archetypal Duffing equation, the addition of one odd superharmonic has demonstrated to succeed in reducing the scattered chaotic region of about 1% (with U 2 /U 1 = 0.8) (Lenci and Rega, 2003c), while one even superharmonic has been able to increase the saved region of about 5% (with U 2 /U 1 = 0.4); moreover, in a microelectro-mechanical system (Lenci and Rega, 2006) one controlling superharmonic has proved to shift the relevant erosion profile of about 8% for U 2 /U 1 = 1.6591 (which is indeed a major costing amplitude), and of about 4% with U 2 /U 1 = 0.5.
After the fundamental resonance frequency The numerical procedure previously described has a very general nature, and can be applied to any kind of global bifurcation with the addition of Figure 37. Erosion profiles of the reference system at ω u = 0.9. any number of superharmonics. To verify its versatility, the numerical control has been applied also at ω u = 0.9, a sample case of the system behavior after the fundamental resonance, for which the analytical control has shown to even deteriorate the safe basin robustness by shifting the smoother part of the relevant erosion backward (see Fig. 31(b) of Sect. 4.1).
Observing the erosion profiles at this frequency, reported in Fig. 37 for the sake of clarity, it is evident that the sharp loss of robustness occurs for quite high values of the harmonic forcing ẋ x S1 ẋ x S1 ẋ x S1 ẋ x  amplitude, highlighted by the red box in figure, after a long smooth decrease, and the final fall down of the profile is definitely stronger than the one detected at ω u = 0.7 (Fig. 31(a)) swiftly leading to the zeroing of the system dynamical integrity. Also in this case, attention has to be focused on the sharp part of the profiles, which marks the loss of safety for the system, so that the global bifurcation to be controlled has to be sought for amplitude values much higher than that associated with the homoclinic bifurcation of the hilltop saddle analytically determined (i.e. U h 1,cr = 0.0027652). A thorough investigation of the system global behavior allows to detect the arise of an internal saddle, coinciding with the birth of a new basin related to a competing 1-period solution, at a very large excitation amplitude (U = 0.604). At this value, the manifolds of the hilltop saddle are already fully intersected (not reported in the next figures for readability reasons), suggesting a crucial role of the internal saddle in the loss of robustness of the system. In fact, the relevant left manifolds W u l (S 1 ) and W s l (S 1 ) undergo a homoclinic bifurcation at U ∼ = 0.686, that from a phenomenological viewpoint corresponds to the basins separation inside the potential well, leading to a strong decrease of the system dynamical integrity. The overall behavior can be recognized in Fig. 38, where manifolds and basins of attraction are reported before, during, and after the homoclinic bifurcation W u l (S 1 ) ∩ W s l (S 1 ). Figure 39(a) is a zoom of the upper panel of Fig. 38(b) around S 1 . By applying the same procedure of the case ω u = 0.7, one even superharmonic is applied to the model. However, this time the ability of control in preventing the homoclinic bifurcation, which is apparent in the separation of the left manifolds of the S 1 saddle shown in Fig. 39(b) for a very low value of the controlling amplitude, is accompanied by a critical modification of the behavior of the corresponding right manifolds. As highlighted by the purple and orange lines in Fig. 39(b), in fact, they intersect each other causing the occurrence of another homoclinic bifurcation W u r (S 1 )∩W s r (S 1 ), which obviously neutralizes the beneficial effect achieved on the left manifolds. This result can be easily explained by observing that in this case, as already noticed in commenting Fig. 9(c), despite the overall asymmetric one-well potential of the system, the in-well basins of attraction for high values of the forcing amplitude are organized in the state plane so as to assume a local topology similar to that of a symmetric two-well potential, which is now playing the role of safe basin to be optimally controlled by adding odd superharmonics (so-called "global" or, herein better, "two-side" control (Lenci and Rega, 2004)). For this reason, to apply the global control at ω = 0.9, one odd superharmonic (3ω u ) is added to the harmonic forcing, unlike what done at ω = 0.7. The implementation of the numerical procedure furnishes the bifurcation threshold of the homoclinic bifurcation for different values of the amplitude ratio U 3 /U 1 (Fig. 40), whose peak occurs for U 3 /U 1 = −0.45 and represents the optimal value of the controlling amplitude to be applied (Settimi and Rega, 2017). The effects of the control application on the basins evolution are reported in Fig. 41, to be compared with those of the reference system in Fig. 38; the evolution of the manifolds is significantly smoother, for both the right and left manifolds of the internal saddle S 1 , confirming the ability of the "two-side" control in modifying the behavior of both the wells (basins in this case) governed by the considered saddle. Moreover, the robustness of the (green) dominant basin is increased by considerably delaying the arise of the competing 1-period solution (blue basin), and consequently the basin separation. As a consequence, the fall down of the "new" safe-basin profile is shifted to higher values of the harmonic forcing amplitude, enlarging the range of U values in which the system safety is guaranteed, while leaving the qualitative trend of the profiles substantially unchanged, as can be observed by comparing the green and black profiles in Fig. 42. From the results, furthermore, it appears that at ω = 0.9 the global control is even more effective than at ω = 0.7, likely due to being farther away from the left-shifted nonlinear resonance (A-vertex of Fig. 6), thus being less influenced by its amplifying effect.

Conclusions
Local and global dynamics of a reduced order model of noncontact Atomic Force Microscope have been addressed in this Chapter, by considering also the possible presence of two different control techniques, in order to comparatively point out their effects on the overall dynamical behavior of the system. The presented results allow to draw some general comments: • The analysis of the basin of attraction evolution, and the quantitative description of the relevant erosion through the evaluation of different integrity measures, prove to be a fundamental tools to obtain a comprehensive characterization of the system dynamics. In particular, the definition of residual integrity levels associated with the system global dynamics has highlighted a marked variability of the classical stability boundary obtained via numerical integration with prescribed initial conditions and, mostly, a meaningful lack of homogeneous safeness of the latter as regards robustness of the system periodic solutions. In contrast, resorting to the tools of dynamical integrity permits to detect thresholds of residual integrity able to ensure acceptable safety targets established a priori according to the required system performances. • The addition of a simple external feedback control into the model, aimed at keeping the single local response of the AFM cantilever to a suitable reference one and representing a simple and efficient procedure for reliably measuring the sample surface, causes a significant modification of both the local and global dynamics. Due to the increased number of degrees of freedom, in fact, the controlled model shows a richer bifurcation scenario than in the original (uncontrolled) system, with elimination of the region of coexistence of resonant and nonresonant solutions and extremely dangerous reduction of the domain of stability of system bounded responses just in the resonance regions where the system behavior has to be taken more strictly under control in practical applications. In these regions, moreover, the analysis of the basins dynamical integrity highlights a generalized detrimental effect of the control on the robustness of the basins of attraction, with small perturbations of parameters which can cause dramatic changes in the system safety. • In a practical perspective, the suppression of the region of coexistence of resonant/nonresonant responses (known as bistability, and investigated mainly in the tapping AFMs literature) can be also considered as a positive effect of control, as it avoids possible discontinuous transitions from low-amplitude to high-amplitude responses which could lead to distortions in the sample scan. However, such a potentially positive effect of control is vanished in parameter space by the arise of an instability region just near the resonance frequency and in phase space, even more dangerous from an operational point of view, by the penetration of tongues of unbounded response inside the potential well which cause a strong reduction in the system safety for very low values of the excitation amplitude. In this sense, charts with constant integrity obtained for several varying parameters represent useful tools not only to assess the effect of a parameter (or a combination of them) on the system robustness, but also to select an appropriate operating setting able to ensure the functionality of the control and thus the success of the AFM scan process. • It is worth underlining that, while dynamical integrity is reasonably easy to investigate in singledegree-of-freedom systems, basins of attraction of dimension larger than two are computationally onerous to implement, and interpreting the results becomes anyway a considerably difficult task.
That's the case of the system with feedback control, for which, to achieve readable results, the investigation is carried out by analyzing suitable bi-dimensional sections of the five-dimensional basins of attraction. Of course, this implies making assumptions about some of the system initial conditions, thus restricting the analysis to selected (although reasonable) cases. An alternative approach is to resort to parallel computing, which anyway requires the implementation of ad hoc routines and the possibility to use powerful digital devices (Kreuzer and Lagemann, 1996;Xiong et al., 2015;Belardinelli and Lenci, 2016). • The second control method applied to the AFM model, which is specifically based on exploiting the global properties of the dynamical system, demonstrates to be able to preserve the system dynamical integrity by reducing the safe basin erosion. The results of Sect. 4 show that the addition of one controlling superharmonic entails a shift of the selected global bifurcation and simultaneously a delay of the fall down of the erosion profiles, thus increasing the system overall safety. Of course, the main challenge in the application of such technique concerns the identification of the global event actual triggering the loss of safety for the system, which calls for the detection of the main saddles with the relevant manifolds. This step of the control procedure is certainly the most computationally onerous, especially in systems (like the one discussed here) with a quite involved topological scenario characterized by tightly rolling manifolds. This is mainly due to the low considered value of the damping parameter (ρ 1 = 10 −3 ) typical of the vacuum operating conditions which the analyses refer to. Nevertheless, the control technique has proved to be effective also in such tricky cases, confirming the robustness and the wide applicability of the method. • As a further aspect of the global control application, the behavior of the various manifolds in the state plane reveals to be crucial for the proper selection of controlling superharmonic to be added into the model, i.e. for the choice between "one-side" and "two-side" control. In fact, the system behavior after the resonance highlights the inadequacy of looking only at the Hamiltonian system as selector benchmark, suggesting to carefully investigate the in-well basins organization as the system parameters change. • From a methodological viewpoint, analysis of dynamical integrity can be used as a general framework for evaluating the effects of control techniques to be possibly employed for securing AFM response robustness and safety against jump to contact in different dynamic conditions. In this sense, some general observations can be pointed out: -Unlike the feedback control, which can be considered as a local control technique, the global one has the explicit objective to increase the overall robustness of the system, an issue which (a) (b) Figure 43. Erosion profiles for the uncontrolled system (black line), for the system with feedback control (orange line) and for the system with optimal superharmonic (thick green line), before the resonance (ω u = 0.7)(a) and after the resonance (ω u = 0.9)(b).
is successfully accomplished with also beneficial effects on the main solutions displayed by the system, whose basins of attraction are either made less fractal (before the resonance) or actually enlarged (after the resonance). -However, although the feedback control has some detrimental effect on the behavior of the system with respect to the response robustness, it is able, when properly working, to exactly reproduce the reference response. Conversely, in the global control, the addition of the superharmonic slightly modifies the amplitude and the shape of the system responses. -Finally, some comments can be drawn comparing the IF (dashed) and GIM (continuous) profiles of the safe basins relevant to the reference, locally controlled and globally controlled systems, respectively, reported in Fig. 43. At the left of the resonance frequency, the erosion develops from the outer edge of the basins, reducing their volume while preserving their compact core. As a consequence, IF is always grater then GIM for all three systems ( Fig. 43(a)). At the right of the resonance, the erosion has the aforementioned development for low values of the forcing amplitude (IF > GIM), whereas in the dangerous initial part of the sharper profiles of both the reference and the globally controlled systems the safe basin develops by stretching its overall shape, as previously shown in Fig. 38(b) and 41(b), thus maintaining the volume almost unchanged while reducing its core (IF < GIM) ( Fig. 43(b)). It is thus apparent that the choice of the integrity measure to be used for the analyses is of great importance, and it must be made on the basis of the system properties, of the features of the erosion involving the basins and of the level of safety to be guaranteed.