Advanced atomic force microscopies and their applications in two-dimensional materials: a review

Scanning probe microscopy (SPM) allows the spatial imaging, measurement, and manipulation of nano and atomic scale surfaces in real space. In the last two decades, numerous advanced and functional SPM methods, particularly atomic force microscopy (AFM), have been developed and applied in various research fields, from mapping sample morphology to measuring physical properties. Herein, we review the recent progress in functional AFM methods and their applications in studies of two-dimensional (2D) materials, particularly their interfacial physical properties on the substrates. This review can inspire more exciting application works using advanced AFM modes in the 2D and functional materials fields.


Future perspectives
Numerous advanced or functional atomic force microscopies have been developed and applied in various research fields, from mapping sample morphology to measuring physical properties. The research and development of AFM-based physical characterization methods require the researchers to have clear natural science concepts and electronic circuit knowledge. The future development direction of AFM-based modes includes the following points: (i) The single physical property measurement method has limitations because the electrical properties are coupled with other properties. Therefore, a method that can explore the mechanical-thermalelectrical-optical coupling is needed. (ii) The functional AFM technology under vacuum conditions needs further improvement. (iii) Optimizing device performance requires understanding the properties of materials under work conditions. A method that can measure physical properties under work conditions is crucial for future research.

Introduction
Investigations at the nano or atomic scale are necessary to understand material properties fundamentally. Scanning probe microscopy (SPM) allows direct characterization and construction of structures at the nano and atomic scales [1][2][3][4]. In addition, numerous advanced or functional SPM methods, including atomic force microscopy (AFM), have been developed and applied in various research fields [5][6][7][8]. SPM detects the near-field interactions between the local probe (or tip) and sample surfaces, such as current and forces, depending on the probe-sample distance and their physical properties (figure 1).

Scanning tunneling microscopy (STM)
The first kind of scanning probe microscopy, the scanning tunneling microscope, was invented in 1981/1982 by Binnig and Rohrer [3,4,9] who went on to receive the 1986 Nobel Prize in physics. STM works using the tunneling effect of quantum mechanics, which states that even if the electron energy is lower than the barrier height, it has a certain probability of passing through the barrier. If the voltage between the tip and sample surface is minimal, we can obtain the exponential dependence of the tunneling current (I) on the tip-sample distance (z), as follows.
where ρ s is the electronic density of states of sample, and k = √ 2mφ/ℏ, φ is the work function of the sample. The main characteristic of tunneling current is the exponentially decaying with the tip-sample distance. This exponential dependence causes the extremely high vertical resolution of the STM, which can reach the picometer regime.

Atomic force microscopy (AFM)
One disadvantage of STM is that it can be used only on conducting samples as the tunneling current is the measured quantity. Instead of the tunneling current in STM, the ubiquitous interactions/forces between the tip and sample are measured in AFM [10]. AFM can be used on conductor, insulators, semiconductors, polymers, and biology cells.
A qualitative sketch of the force between the tip and sample is given as a function of the tip-sample distance (figure 2(a)). Three different regimes can be distinguished. If the tip is far away from the surface, the force between the tip and sample is negligible. An attractive (negative) force between the tip and sample occurs for closer distances, marked by blue color. For very small distances, a strong repulsive force exists between the tip and sample, marked by the red line.
There are three basic AFM working modes, which are primarily used to obtain the surface topography. The first mode is contact mode (figure 2(b)), which is a static mode. In this mode, the tip contacts the sample surface during scanning, and the interaction force between them is the repulsive force, which is measured by the cantilever deflection signal. The zfeedback loop maintains a constant force between the tip and sample. The corresponding changes in the z-position for maintaining a constant tip-sample distance (i.e. constant force) correspond to the sample surface topography.
The second mode is the amplitude modulation (AM) mode, also known as the Tapping mode (figure 2(c)), which is different from the static mode. The AM mode is a dynamic mode, and the feedback loop is the oscillation amplitude of the AFM probe. Herein, the micro cantilever is forced to vibrate, and the frequency ω drive is its free resonance frequency. When the tip and sample approach each other, the oscillation amplitude at the fixed excitation frequency ω drive will change. The preamplifier converts the defection signal from the photodiode to a voltage signal proportional to the cantilever deflection. The AC voltage signal has a frequency ω drive and amplitude proportional to the oscillation amplitude of the vibrating cantilever. Using a lock-in amplifier, the amplitude of the AC signal at frequency ω drive is measured. Moreover, maintaining a constant oscillation amplitude corresponds to maintaining a constant tip-sample distance. Thus, the z-feedback signal is used as the height signal to map the topography during data acquisition.
The third mode is the frequency modulation (FM) mode, also called non-contact dynamic mode (figure 2(d)) [11]. In this mode, the resonance frequency shift is measured due to the tip-sample interactions. The z-feedback loop maintains a constant frequency shift of the cantilever resonance frequency as the setpoint. A schematic of the implementation (figure 2(d)) consists of an oscillator loop, to which the measured oscillation signal is fed back (after a phase shift) as the driving signal of the cantilever. Please see details in [1].

AFM probes
For sensing normal tip-sample forces, the force sensor should be rigid in two axes and relatively soft in the third axis. This property is fulfilled with a cantilever beam. The first cantilevers were made from a gold foil with a small diamond tip attached to it. Later, silicon micromachining technology was employed to build cantilevers in parallel production with welldefined mechanical properties. The first micromachined cantilevers were built at Stanford in the group of Calvin F. Quate. Initially, mass-produced cantilevers were built from SiO 2 and Si 3 N 4 . Later, cantilevers with integrated tips were machined from silicon-on-insulator wafers. The most common cantilevers in use today are built from all-silicon with integrated tips pointing in a [001] crystal direction. Figure 3(a) shows the type of cantilevers that are mainly used today: micromachined silicon cantilevers with integrated tips.
In dynamic atomic force microscopy, some requirements for the force sensor are similar to the desired properties of the time-keeping element in a watch: utmost frequency stability over time and temperature changes and little energy consumption. Around 1970, the watch industry was revolutionized with the introduction of quartz tuning forks as frequency standards in clocks. Experimental studies of using quartz-based force sensors were carried out soon after the invention of the AFM. Tuning forks were used as force sensors in acoustic near-field microscopy, scanning nearfield-optical microscope. Figure 3(b) shows a quartz cantilever based on a quartz tuning fork [12]. In contrast to micromachined silicon cantilevers, the quartz forks are large. Hence a wide selection of tips can be mounted on a tuning fork with the mere help of tweezers and a stereoscopic microscope-sophisticated micromachining equipment is not needed.
The tip apex, which is the key sensing element, must be small for the highest spatial resolution. Nowadays, with a COfunctionalized tip attached to a quartz tuning fork sensor, the high resolution structure of 2D materials have been obtained, which provides support for in-depth understanding physical properties of new materials [13][14][15].
Although CO-functionalized tip was ultrathin and has high resolution, it still not an ideal tip. An ideal tip would be a sharp spike, such as a column of single atom width and infinite length. However, even this is not ideal, considering the force would act on the entire column in addition to the apex atoms. An ideal tip would consist of an isolated atom or a cluster of atoms, which is not realistic.
A more typical tip for basic AFM imaging is usually a cone or pyramid (table 1 and figures 3(a) and (b)), with an apex ratio between 1:4 and 1:10. The tip is a rounded tip, usually approximated by a sphere with a radius of 5-20 nm. The probe is typically made of silicon (figures 3(a)-(c)) or silicon nitride (figure 3(d)) using methods typical for micro-fabrication, including lithography, sputtering, anisotropic etching, or focused ion beam milling. Probes can  [12]. © IOP Publishing Ltd. All rights reserved. be modified by doping and coating thin films to enhance the contrast of a particular analytical method. Figure 3 shows the SEM images of some typical commercial AFM tips. The high aspect ratio tip (figure 3(c)) is suitable for performing measurements on samples with a sidewall angle approaching 90 • . The tip has an overall height of 10-15 µm, which allows for measuring highly corrugated samples. At the last few micrometers, the tip exhibits a high aspect ratio portion that is symmetric when viewed from the side and along the cantilever axis. The tip radius is typically 10 nm. The pyrex-nitride probe (figure 3(d)) is hydrophobic and used for biological sample detection. Cantilevers have been made from silicon and silicon nitride. However, in some cases other materials, such as diamond with outstanding mechanical properties, have been used to obtain specific cantilever and tip properties (table 2).
Two important properties of the cantilevers are their free resonant frequency (f 0 ) and spring constant (k). Determining the cantilever resonant frequency is straightforward using resonance tuning. However, the determination of the cantilever spring constant is relatively complicated. Currently, the commercial AFM (Asylum Research) provides thermal tune method to quickly determine the f 0 and k of each probe. The resonance frequency range of common probes is listed figure 4, which provides a reference for probes selecting in the experiments. Different from the above basic modes, much more complex AFM probes are needed in the functional AFM modes, like scanning thermal microscope (SThM), scanning near field optical microscopy (SNOM), scanning microwave impedance microscopy (sMIM) and so on. In the following section 1.4, we will introduce probes appropriate for each advanced AFM methods.

Multi-functional AFM
Numerous multi-functional AFM-based techniques have been developed to characterize the local mechanical properties, electrical properties, optical characteristics and thermal properties of the sample. The methods include the scanning Kelvin probe microscopy (SKPM) and electrostatic force microscopy (EFM) for the surface potential (SP) and electrostatic potential measurement, dual AC mode and contact resonance (CR) mode for the contact stiffness and viscous response detection, scanning thermal microscopy (SThM) for the thermal response, and SNOM for the photonic optical characteristics detection. Several typical AFM-based techniques and their applications are presented in table 3.

Friction force microscopy and transverse shear
microscopy. FFM and TSM are the derivatives of the contact mode of AFM. They maintain the normal force constant and record the lateral signal between the tip and sample surface (figure 5) [18,19,[66][67][68]. The lateral bending or twisting arises from forces parallel to the sample surface plane and acting on the AFM-tip. In TSM, the scanning direction of a force microscope probe tip is parallel to the cantilever axis, and the lateral deflection or twist of the cantilever is recorded. This mode of operation differs from the better-known FFM  [20] Multi-harmonic -AFM Higher harmonic components introduced by the nonlinearity of the tip-sample interaction forces Interlayer coupling [7] and local variations of the elastic modulus [7,20] CR mode Contact resonance frequency and quality factor Contact stiffness and viscous response detection [21] C-AFM Current Local conductivity [22][23][24][25] SKPM/KPFM Electric potential Surface potential (SP) [26], work function [27], contact potential difference [28], charge transfer [29,30], surface point defect/adsorbate [31][32][33], voltage drop [34], and capacitance coefficients [35], EFM Electrostatic forces Capacitance coefficients [35], SP [36], and dielectric response [37] MH-EFM Electrostatic forces SP [38], work function [39], and mobile charge carriers (MCC) [40] sMIM Microwave reflection Dielectric constant [41,42], conductivity and permittivity variation [43,44], charge carrier variations [45], and doping concentration [  or lateral force microscopy (LFM) technique in that the scanning direction is perpendicular to the cantilever axis in the FFM/LFM technique [68]. Comparing to FFM/LFM, TSM has enhanced sensitivity to elastic deformation properties and can reveal crystallographic orientation. Owing to the flexibility of 2D materials, friction can easily drive stretch deformation, leading to in-plane mechanical (shear stiffness) symmetry breaking. In another word, the shear stiffness with uniaxial strain has significant lattice dependent anisotropy, especially for 2D atomic crystals with centrosymmetric structures. Based on normal tension straininduced shear anisotropy, we have reported that the visualization of crystallographic-orientation-dependent shear deformations in flexible, hexagonal 2D atomic crystals by using TSM. (See details in figure 4 in [18].)

Multi-harmonic and dual AC mode.
Conventional dynamic AFM methods focus on the excitation and detection of a single frequency component of the tip motion, usually at the fundamental resonant frequency. Information on sample properties included in other frequency components, such as harmonics, is irreversibly lost, which limits the capabilities of the AFM [5]. High harmonics have been recognized, but their utilization has been limited due to theoretical and experimental complexity [69,70]. Hence, developing multi-frequency AFM (MF-AFM) technology is imperative, based on simultaneous excitation and/or detection of several frequencies of the probe's oscillation [71,72].
In conventional AM AFM, the cantilever deflection is given as where A, z 0 , and ϕ are the amplitude, static deflection component, and phase shift with respect to the driving force, respectively. Owing to the presence of multifrequency components (harmonics) introduced by the nonlinear interaction between the tip and sample, the deflection of the cantilever is as follows.
where A n is the amplitude of the harmonic with angular frequency nω. Equations (2) and (3) are compatible with a pointmass description of the cantilever. A more precise description is achieved by considering the extended characteristics of the cantilever, with the probe deflection containing contributions from all its eigenmodes (q j ), In equation (4) the eigenmodes have been expressed in terms of the different harmonics (A n ). For multi-harmonic mode, the amplitude of higher harmonics can be expressed as For Q values greater than 5 and ω = ω 0 , where F ts represents the tip-surface force; d is the instantaneous tip-surface separation; and ω 0 is the fundamental resonant frequency. For dual AC mode, the instantaneous tip position z when interacting is described approximately by equation (7), which neglects the contribution of the other modes and harmonics where A 1 , ϕ 1 , and A 2 , ϕ 2 are the oscillation amplitudes and phase shifts at ω 1 and ω 2 , respectively. Thus, the total energy dissipated per cycle E dis can be approximated as the sum of the two eigenmodes where Q i , k i , and A 0i are the quality factor, force constant, and free amplitude of the ith eigenmode, respectively. The second eigenmode frequency is 6.27 times the first eigenmode frequency, i.e. n = 6.27. Equation (8) links the second modal observables, such as A 2 , with nonconservative interactions (E dis ). The schematics of the multi-harmonic AFM and dual AC modes are shown in figures 6(a) and (b), respectively. The typical AFM probes used in dual AC and multi-harmonic AFM mode measurements are shown in figures 6(c)-(f). In multi-harmonic mode, a special harmonic cantilever is fabricated for tuning of its second eigenmode, to six times the first eigenmode to increase the signal to noise level, as shown in figure 6(c). In dual AC mode, the silicon AFM probes (AC240, Asylum Research) (figure 6(e)), having first resonance mode f r1 of ∼76.3 kHz and second resonance mode f r2 of ∼451.4 kHz are used.
High harmonic supports the research on cantilever's dynamic behaviors with high compositional sensitivity and atmospheric resolution. Expressing the tip-sample interaction as a function of high harmonic amplitude or frequency has several advantages. First, the high harmonic method can provide more details and intrinsic characteristics of the sample because it depends on the nonlinear interaction and is measured with zero references. Second, it allows a high spatial resolution, particularly with small oscillation amplitude and low bandwidth reduction [74,75]. Thus, scanning with a relatively appropriate speed is possible by combining the small-amplitude higher-harmonic AFM with the resonanceenhancement technique. Third, high harmonic can be joined with other MF-AFM methods to develop interesting applications. Measurement innovations facilitate developing AFM and nanoscience. The high harmonic dynamic AFM will be attractive to nano-and pico-analytics in future physics, chemistry, biology, and materials science. 1.4.3. Contact resonance (CR) mode. The CR mode is an advanced AFM technique that measures the resonance peak of the probe and sample system. It offers improved nanomechanical characterization of layered materials compared to the existing techniques [76]. CR-AFM is achieved by introducing a slight vertical modulation to either the cantilever base [77] or sample [78,79]. The cantilever CR frequency (CRF) and quality factor (Q) change in response to the viscoelastic properties of the sample when the tip is in contact with the sample. It can quantitatively characterize the viscoelastic response of materials with high sensitivity and large dynamic detection range because of the high-frequency operation and resonance enhancement effect [80,81]. In this study, the detection limit of AFM used was approximately 10 pm; therefore, the excitation amplitude of the tip-sample contact was reduced to the picometer level, i.e. considerably smaller than the interlayer distance of two dimensional (2D) layered materials. Figure 7 shows the schematic and mechanical model of CR-AFM and related definitions and characterizations for quantitative analysis. CR-AFM is based on the typical contact mode imaging method; the key is the vertical modulation with a slight amplitude introduced to the sample, as illustrated in figure 7(a). The modulation is operated at a high frequency and hence does not affect the contact-feedback imaging; however, the modulation can be coupled to the cantilever deflection that can be extracted for detection. The high-frequency signal is continuously changing with the sample mechanical properties as the tip scans the sample in contact mode, detected by a lockin amplifier and referred to as CRF and CR amplitude (CRA). Figure 7(b) depicts the mechanical model of the tip-sample dynamic contact in CR-AFM. Because the tip location on the cantilever directly affects the final measurement, the cantilever is modeled as a distributed mass with spring constant k L rather than a point-mass approximation. The tip location is considered as a relative position γ = L 1 /L (with 0 ⩽ γ ⩽ 1) from the clamped end. The tip-sample interaction is considered as normal elastic and dissipative (damping) forces, modeled as a spring and dashpot in parallel (i.e. Kelvin-Voigt model). The variables k * and c indicate the contact stiffness and damping, respectively. The tip-sample stiffness and damping variations are reflected in the CRF and Q of the coupled system, which can be converted to the sample elastic and loss moduli with right contact mechanic models, respectively. Higher contact stiffness indicates high CRF, and larger viscosity (dissipation) at the tip-sample contact indicates low Q.
Consideration of the cantilever dynamic model is primarily based on the Euler-Bernoulli beam theory. The contact stiffness k * is calculated using the cantilever stiffness, relative tip position, and resonance wavenumber [21], as described in equations (9)- (12). The relative tip position γ is characterized by SEM, as shown in figure 7(d), The parameter A n , related to the wavenumber x n L of the cantilever resonance in free space, is defined by The wavenumber is denoted as y n L under the tip-sample coupling condition: Based on the equations (9)-(11), the contact stiffness k * can be determined as: k L (y n Lγ) 3 1 + cos y n L cosh y n L D , where The wavenumber x n L for the flexural mode n is the solution of the characteristic equation 1 + cos x n L cosh x n L = 0 for cantilever vibration in free space; the first two roots x 1 L and x 2 L are 1.8751 and 4.6941, respectively. Using the contact stiffness, the sample modulus can be extracted quantitatively based on contact mechanic models [21], such as the Hertzian (Hertz) contact, Derjaguin-Muller-Toporov (DMT), and Johnson-Kendall-Robert (JKR) models. The two frequently used models in CR-AFM are the Hertz and punch contact models, which depend on the shape of the tip (sphere or punch) [21]. The discussion here is restricted to the Hertz and DMT models as a spherical tip is used in all experiments. The main difference between the two models is the consideration of the adhesion force F ad . Figure 7(c) depicts the Hertz contact of the sphere tip and sample and defines the related quantities (F N , R, a c , and δ). The normal contact stiffness of the tip-sample system can be expressed by equation (13), which is related to the contact radius a c and reduced Young's modulus E * .
where E s and E t are Young's moduli of the sample surface and AFM tip, respectively. When F ad is considered in the measurement, the DMT model is more accurate. The contact radius a DMT is modified as Finally, the elastic modulus of the test sample is determined from the experimental values of the contact stiffness using equations (16) and (17).
The indentation can be deduced from the DMT model and is described as equation (18), In regard with the CR-AFM, dual AC resonance tracking, invented by Asylum Research [82], is an exclusive implementation of the CR. It tracks the CR ( f 0 ) by adjusting the two drive frequencies ( f 1 , f 2 ) of the cantilever to zero difference amplitudes (A 2 −A 1 = 0) by reducing the crosstalk, rather than using the phase as frequency feedback.
The work function φ and V CPD are significant parameters for understanding the working principle of SKPM. Herein, φ is defined as the energy difference between the Fermi energy and vacuum level (figure 8), and V CPD is defined as the difference in the φ of two materials, as: where e is the electric charge; and φ sample , and φ tip are the work function of the tip and sample, respectively. Figure 8 shows the energy level of the tip and sample surface with no electrical connection. The vacuum levels are related, but Fermi levels are dissimilar. Equilibrium needs Fermi levels to line up at a stable state if the tip-sample is sufficiently close for electrical contact. The Fermi energy level will align by electron current flow and the system will accomplish an equilibrium state ( figure 8(b)). If the applied V DC has a similar magnitude as V CPD but with a reverse sign, the applied voltage eliminates F el (figure 8(c)). Thus, when the tip work function is known, the sample work function can be obtained, and SKPM can measure V CPD . Figure 9 reveals that SKPM operates in the nap mode to measure the sample SP. During the first pass mode, the morphological image is obtained the same as in standard tapping mode. In the second pass mode, the tip is raised on the sample surface by a desired lift height. The mechanical excitation is switched off, and V AC and V DC are applied between the tip-sample [93]. The electrostatic force F el between the tipsample is given by: where ∂C ∂z and ∆V are the capacitance coupling and voltage difference between the AFM tip and sample, respectively. The total voltage difference between the tip and sample is where ω e is the angular frequency of the applied AC voltage. The F el applied to the AFM tip is: where During SKPM scanning, ω e is generally set to ω r (i.e. F ωr = F ωe ) to achieve a higher signal-to-noise ratio. The V CPD acquired at minimum F ωr shows the sample's local [89]. In the first scan, the morphological image is captured using mechanical excitation of the cantilever. The surface potential mapping is acquired in the second scan. The electronic circuit diagram of the first and second scans is given by the black and purple lines. Reproduced from [91], with permission from Springer Nature.

Conventional EFM.
Conventional EFM operates in a two-pass process. During the first scan, the topography image is mapped the same as the standard AC mode. During the second scan, the tip is raised by lift height (usually z is approximately 20-50 nm) to a controllable height z above the sample surface at a specific distance. DC voltage is applied between the tip-sample in lift mode, and the phase shift of the vibrating cantilever is monitored to explore the electrostatic gradient of the sample [91, 95, 100].

Multi-harmonic EFM (MH-EFM).
In MH-EFM (figure 10), the dual-pass mode detects the topography and electrical properties of the sample. In the first pass, the AFM works similar to the AC mode. In the second pass (lift mode), the short-range repulsive forces and van der Waals forces can be ignored because the tip-sample distance is sufficiently large; meanwhile, an AC voltage of frequency f ω is applied between the tip and sample, thereby resulting in multi-harmonic electrostatic force F, as follows.
Accordingly, the second and third harmonic components are given by Thus, the information corresponding to ∂C ∂z and ∂C ∂V can be obtained by detecting the 2 ω and 3 ω components, respectively. Further, A ω is proportional to SP, A 2ω is related to the mobile charge carrier density (MCD), whereas A 3ω reflects the capacitance (or dielectric constant) of the sample [38,101,102]. This technology provides an opportunity to explore the electrical response of samples in the kHz range.
During the operation of sMIM, a microwave signal at a high frequency (few GHz) was applied to the tip apex and subsequently transmitted to the sample (figure 10) [115]. Some of the microwave signals pass through, whereas the remaining are reflected off the surface. The reflection coefficient Γ of microwave signals from the contact point measuring the tipsample impedance mismatch is as follows.
where Z 0 and Z tip-sample represents the characteristic impedance of the transmission line system and tip-sample, respectively. A cancellation signal suppresses the background to amplify small changes without output saturation. Then, the amplified signal is demodulated by mixer M1. The effective tip-sample impedance can be divided into real and imaginary parts. Standard dielectric sample (Al 2 O 3 @SiO 2 ) calibrates the  [105]. The TiW/Au metal tip on the free end of the cantilever is connected to the wire bond pad by a conducting path buried inside two Si 3 N 4 layers. Both the front (tip side) and backsides of the cantilever are covered by shield metals, which are electrically grounded in the microwave measurements. The SEM images in figures 10(e) and (f) show the bond pad, the cantilever and the pyramidal metal tip.
We have review the basic principles and instrumentation of sMIM and related scanning microwave microscopy (SMM) in [118], and also discuss its widespread applications in electrical imaging of a number of novel materials and biological systems. Recently, a nice review [119] about sMIM were published, and outlines future opportunities in expanding the capabilities of sMIM. Please check these reviews for further information.

Scanning thermal microscope (SThM).
Understanding energy dissipation at the nanoscale level requires probing the temperature fields with micron/nanometer resolution. Scanning thermal microscopy (SThM) [120,121] can quantitatively map temperature fields or thermal properties by scanning a sharp SThM probe with a temperature sensor at the tip. This technique has been applied in diverse areas, including microelectronics [122][123][124], optoelectronics [125,126], carbon nanotubes [127,128], and 2D material [53], since the 1990s.
The setup of an AFM-based SThM system is illustrated in figure 11. The combination of the X-Y scan position data along with the force feedback and thermal signals measured by the sensor located either at the tip or on the cantilever provides the raw data of the surface topography and thermal images. The thermal image contrast reflects the change in the amount of heat locally exchanged between the tip and sample. The force feedback control system operates simultaneously but independent of the thermal measurement. The thermal control unit performs real-time thermal signal analysis [129].
Since 1990s, various thermal methods using different thermosensitive sensors or phenomena have been developed. According to the temperature-dependent mechanism, they are classified as thermovoltage [130][131][132], change in electrical resistance [133][134][135][136], fluorescence [137,138], or thermal expansion [139]. These methods use functional SThM probes, such as metallic or doped silicon. Herein, we introduce two kinds of probes and their corresponding measurement methods. Null-point (NP) SThM [140][141][142] is based on the thermodynamic principle that the heat flux through the tipsample thermal contact is zero when the probe tip and sample surface temperatures are identical. Hence, NP SThM measures the unperturbed temperature by the heat flux through the tipsample thermal contact. The experimental setup for the NP method is shown in figure 12. The signal-to-noise ratio is maximized by amplifying the thermoelectric voltage extracted from the Wheatstone bridge using a preamplifier and reducing the 60 Hz harmonic noise via a noise filter. The DC thermoelectric voltage is fed into the signal access module and is available with the topographic signal from the atomic force microscope. A DC power supply is used to Joule-heat the sample [141]. The SEM images of the SThM probe used in the experiment are shown in figure 12(b). The cantilever made of silicon oxide, with the lowest thermal conductivity among the materials suitable for nanobatch-fabrication, maximizes the SThM probe sensitivity. Figure 12(c) shows the principle of quantitative thermal profiling. The SThM probe with a nano-thermocouple junction integrated at the end of the tip, as shown in the inset, is mounted on an atomic force microscope and scans on the same line in the contact and nonthermal contact modes. Thus, the difference between the signals obtained in the two modes is owing to the heat flux through the tip-sample thermal contact. Considering the governing equation for the temperature distribution and boundary condition [143], the local temperature of the sample T s is given as where k slope is a dimensionless constant determined through the measurement process in figures 12(d) and (e). Moreover, T c and T nc are obtained from the temperature sensor integrated at the tip when scanning in the contact and nonthermal contact modes, respectively.  Figure 13(a) shows a typical thermal probe made of a high-doped silicon cantilever with a low-doped silicon heater region. As the heater region is thermosensitive, an electric current to the probe controls and monitors the temperature of the integrated heater region. Changes in the ambient thermal resistance (TR) can change the heater region temperature (∆T probe ) and are monitored through the voltage value of the probe (V probe ). In the tip-sample system, the total TR (TR total ) measured by SThM between the tip and sample consists of  (e) Temperature jumps at different sample temperatures. The slope of the graph, k slope , is 10.9 K K −1 for this particular probe. Reprinted from [141], with the permission of AIP Publishing. Reprinted with permission from [143]. Copyright (2011) American Chemical Society.
the tip TR (TR tip ), contact TR (TR contact ), and spreading TR of the sample (TR sample ), as shown in figure 13(b). The circuit diagram of the local thermal analysis (ZThermal) module is shown in figure 13(c). The electric resistance of the reference electrical resistor is known (R ref = 4000 Ω), and the total voltage (V total ) and probe voltage (V probe ) can be exerted or detected. The electric resistance of the probe is thermosensitive, and the relationship of the electric resistance (R probe ) versus temperature (∆T probe = T probe -T air ) is shown in figure 13(d). The maximum electric resistance is observed at 550 • C. Below 550 • C, the electric resistance of the probe corresponds to a specific temperature, as follows [53]: The thermal power generated by the thermal probe can be calculated as follows: We introduce the design idea of NP SThM [144] and the quantized thermal transport characterization [145] and use the thermal tip-sample approach curve to calculate TR. The model assumes that when the tip abruptly contacts the sample, the thermal circumstance and heat transfer coefficient (h) of the thermal tip in the air do not change. However, before contact, the heat (Q air,off ) generated by the thermal tip dissipates into the air, whereas after contact, a new thermal flow channel is opened and a part of the heat (Q sample ) diffuses into the sample. Thus, the heat transfer equation can be given as On contact: where A is the superficial area of the tip; and P off and P on are the heat power generated by the thermal probe before and after contact, respectively. Thus, the total TR can be deduced as: Based on the thermal tip-sample approach curve (figure 13(e)), the total thermal resistance can be solved using the tip voltages of points A and B, combined with equations (29)- (33). Although the SThM method has been calibrated and quantified in air conditions, two main limitations exist. One is the spatial resolution due to the liquid film existing at the tip-sample interface, and the other is the difficulty in quantitative measurements due to the parasitic air conduction between the sample and SThM probe. However, the ultra-high vacuum (UHV)-based SThM technique enables the quantitative measurement of thermal maps with ∼15 mK temperature resolution and ∼10 nm spatial resolution [146].
SThM opens new avenues for investigating physical phenomena in the quantum regime. In this review, we will mainly focus on the applications of SThM at the interface and surface of 2D materials.

Applications of AFM in surface/interface of 2D materials
Following the discovery of graphene [147], considerable research has been focused on 2D materials [148,149]. Transition metal dichalcogenides (TMDCs) have emerged as an important 2D-layered material owing to their superior characteristics [149]. Several desirable properties emerge at monolayer limits, the most notable of which is the presence of a direct bandgap [149,150]. From the perspective of spatial resolution, AFM has obvious advantages in studying the ultrafine structure of 2D materials. Moreover, the multi-functional AFM technology, which was extended from AFM, can realize high-resolution and physical and chemical properties quantitative characterization of materials. It is playing an increasingly important role in characterizing the local electrical, mechanical and thermal properties of 2D materials at the nanoscale level [17-19, 53, 151-153]. This section discusses recent studies on surface/interface of 2D materials using AFM.

Mechanical engineering of surface/interface
Tuning the material band structure by subjecting it to strain is an important strategy to enhance the electronic device performance. Diverse methods have been proposed for introducing strain in 2D materials, including the bending of films on elastic substrates [10,11,16], stretching of films using atomic force microscope probes [17][18][19], and thermal expansion mismatch [20][21][22][23]. Our review mainly focuses on the stress application methods discussed in sections 2.1.1 and 2.1.2.

Thermal strain.
The thermal expansion coefficient (TEC) mismatch is a convenient method to produce controllable strain in TMDCs [155,156]. Owing to the high temperatures when synthesizing TMDCs, the TEC-mismatch between the substrate and 2D semiconductor can control the strain in the synthesized 2D materials, as shown in figure 14. The tensile (compressive) strains accumulate on the TMDCs when the substrate TEC is less (greater) than that of the 2D material. Moreover, relaxed samples are obtained if the TEC of the substrate and 2D material match. The thermal strain-engineered behaviors of chemical vapor deposition (CVD)-grown triangular MoS 2 and WS 2 flakes were reported by our group and will be summarized in this section.

Edge dependent strain behaviors.
We consider a micrometer or larger-sized 2D layer as a continuum where the atomic details are averaged, and the classical elasticity theory governs, as illustrated using the six-fold hexagonal 2D material in figure 15(a). However, a close examination  Process schematic through which 2D materials realize strain via the TEC-mismatch between the substrate and TMDC. (a) Tensile strain is achieved when the substrate TEC is less than that of the 2D material, (b) relaxed samples are achieved when the TEC of the substrate and 2D material match, and (c) compressive strain is achieved when the substrate TEC is greater than that of the 2D material. Reprinted with permission from [154]. Copyright (2022) American Chemical Society. of the hexagonal flake at the nanoscale level, e.g. monolayer WS 2 , under strain ( figure 15(b)), shows that the two inequivalent edges can make the local elastic modulus nonuniform. The hexagonal flake is divided into six regions because of its C 3 symmetry and the two different edges. We denote these regions using the edge type as Zigzag (ZR) and Klein regions (KR). The CR mode characterizes the out-of-plane modulus by recording the resonance frequency shift of the probesample system ( figure 17). Based on the model explained in and (c)) show uniform contrast but different elastic modulus ( figure 17(e)). The out-of-plane elastic modulus of ZR was significantly smaller than that of KR, which is consistent with the Density Function Theory (DFT) results ( figure 17(a)). Furthermore, the ZR and KR energy dissipation differences are resolvable through CR ( figure 17(f)).
The friction force between the tip and sample is often affected by the out-of-plane, in-plane, and buckling modulus of a sample ( figure 17(g)), which are difficult to measure. Therefore, FFM was employed here to measure the friction force by imaging the lateral torsion of the cantilever through contact scanning along the perpendicular direction. Figures 17(h) and (i) show the friction force images of the WS 2 flake. The friction force of KR is larger than that of ZR, probably because of the smaller in-plane and larger out-ofplane moduli of KR. Thus, the micron-sized triangular KR and ZR explicitly show different measured mechanical properties in terms of elastic modulus, energy dissipation, and friction force. Advanced AFM techniques further investigate the strainengineered structure and property of MoS 2 flakes. The AFM topography image of the triangular sharp-corner MoS 2 monolayer shows no specific features within the flake, as shown in figure 19(a). In contrast, FFM (figure 19(b)) reveals three sharp corners (marked II) with distinctive strain-induced features from the central region (marked I). Moreover, the strain at the corners (II) is different from that at the central region (I).
According to the previous theoretical works, the tensile (compressive) strain can decrease (increase) the band gaps of the MoS 2 monolayer [161]. The DH-EFM qualitatively evaluates the strain-induced local MCD, which depends on the band gaps of MoS 2 flakes. A higher MCD is observed at the sharp corners than in central regions within the triangular flakes ( figure 19(c)), indicating a smaller bandgap at the corners under a larger tensile strain. The sMIM mode measures the electronic properties of the sharp-corner MoS 2 flakes ( figure 19(d)) and confirms a higher electronic conductivity at the sharp corners. The results indicate that the sharp corners suffer larger tensile strain than the central regions.
The strain-engineered structure and properties of the veinlike MoS 2 flakes were further investigated in comparison with the sharp-corner flakes (figures 19(e)-(h)). The distinctive strain-induced features were observed using FFM (figures 19(f) and (g)), with three straight lines originating from the center to the corner apex. Figure 19(g) shows the enlarged central region, with a hierarchical hexagonal ripple pattern around the central triangle. In DH-EFM (figure 19(h)), the highest mobile charge densities were observed at the three straight lines starting from the flake center, in contrast to the medium and lowest charge density at the ripple areas and corners, respectively. Thus, the three straight lines suffer larger tensile strain than the ripple areas.
Kim et al [158] presented nanoscale photoluminescence (PL) spectroscopy images of triangular CVD-grown WS 2 monolayers of different sizes under different temperatures and excitation power ( figure 20(a)). Intense PL emissions were observed around the edges of individual WS 2 grains. Meng et al [162] reported the SPs and work functions of strained WS 2 flakes of 20 µm size. In the interior regions, work functions are much larger than those of the unstrained regions (figures 20(b) and (c)). The strain domain distribution displays a three-fold symmetry in figures 18-20. However, the strain domains were more complex when the TMDC flakes increased to ∼90 µm.
The monolayer WS 2 flakes were also prepared via CVD and underwent fast-cooling to introduce thermal strain ( figure 21(a)). The magnitude of the tensile thermal strain within the WS 2 /SiO 2 flakes depends on the TEC-mismatch and flake size. Here, the WS 2 monolayer (size ∼100 µm) was split into different forms by several cracks, and the cause of the crack formation was studied using finite element analysis (FEA). The von Mises stress is a critical principle to characterize the yielding and fracture of materials; a point subject to a higher von Mises stress is the starting point of failure. The von Mises stress distribution of a triangular WS 2 flake with a hexagon nucleation center is shown in figure 21(b). Strong stress concentrations are generated at the vertices of the hexagon nucleation center (marked by arrows), which initiate cracks that extend outward to release the stress. Because the stress in the corner regions of the WS 2 layers is much larger than that in the edge centers, the cracks tend to extend toward the corner regions. Several cracks (broken lines) were found in the strained WS 2 flakes, and these cracks split the WS 2 flakes into six-splitting forms ( figure 21(c)).  Reproduced from [153]. © IOP Publishing Ltd. All rights reserved.
The electrical properties of the triangular-shaped sixsplitting WS 2 monolayers were examined. Figure 21   each part of the six-splitting WS 2 monolayer (figure 21(e)), thus suggesting nonuniform local electrical properties. Additionally, the band gap, frictional, viscosity, and elasticity characteristics of the different strain regions were investigated. The nanopattern should enable the flexible designing of more sophisticated devices based on 2D materials.

Geometry dependent strain behaviors.
The straininduced hierarchical ripple nanostructures, modified by the MoS 2 shapes, were also observed in MoS 2 flakes. Therefore, the structural evolution of these hierarchical nanoripples is further discussed based on the geometry and thickness of MoS 2 flakes.
The shape of the MoS 2 flakes can be modified by controlling the S/Mo ratio and growth temperature, as shown in figure 22(a). The most common triangular flakes are obtained at a high S/Mo ratio and low growth temperature. When the S/Mo ratio is reduced, and growth temperature is increased and triangular flakes with multi-apex corners are obtained, as shown in figures 22(b)-(g). The optical image of each flake shows a uniform contrast. Furthermore, a second layer can grow at the center of these triangular flakes and exhibit a perfect triangular shape. Herein, both 2H-and 3R-MoS 2 bilayer polytypes were obtained, as shown in figures 22(h) and (i). After CVD at high temperature, the isotropic tensile strain was applied on the triangular MoS 2 flakes through the substrate during fast-cooling. The tensile strain is due to the TEC mismatch between the MoS 2 flakes (10 −5 k −1 ) and SiO 2 /Si substrate (10 −7 k −1 ) and may induce the out-of-plane ripples due to the Poisson instability of MoS 2 layers. Considering the geometry of MoS 2 flakes, the strain at the three corners is larger than in other regions.   The universal self-similar hierarchy of ripples can be observed within the thin sheets or films when under tensile strain and boundary edge confinement [164][165][166][167]. This curtain effect is demonstrated in our experiment, shown in figure 23(h). These ripples are primarily along the longitudinal direction (same as the direction of gravity) and perpendicular to the boundary edges. Wrinklon is the localized transition zone where two ripples merge, which is a building block for these hierarchical ripple patterns. To further understand their structural evolution, we gradually bent the boundary edge, as shown in figures 23(e)-(g). When bent to a '∧' shape, the ripples fail to orient along the longitudinal direction but are affected by the bent edges. The ripples near the edges orient perpendicular to the bent edges. In contrast, the ripples farther from the edges gradually bend along the longitudinal tensile direction of gravity while retaining hierarchy. The hierarchical ripples in the bent curtains are relatively similar to those in the triangular MoS 2 flakes. Based on their inherent similarities, a simple diagram is proposed to elucidate their formation mechanism (figures 23(i)-(l)).

Strain-induced in-plane anisotropic shear behaviors
(crystallographic orientation imaging). Xu et al [18] reported in-plane shear stiffness anisotropy under uniaxial normal tension strain in monolayer MoS 2 and experimentally verified the shear characteristic based on friction-driven stretch deformation during the TSM contact scan.
For TSM, the scan direction of the AFM tip is parallel to the cantilever axis, and the lateral torsion of the cantilever is recorded. For a flexible film weakly bound onto a rigid substrate, the film puckers locally during contact scan by the AFM tip. When the tip moves the puckered region forward, the puckered geometry causes the film to relax at the front edge and stretch at the rear area of the tip simultaneously ( figure 24(a)), which breaks the original in-plane mechanical isotropy of the hexagonal 2D atomic crystal. FFM can characterize anisotropic stretch deformation and friction energy dissipation for an elastically anisotropic film. Such frictiondriven stretch deformation will cause the stretch force direction to deviate from the deformation direction. The intersection angle creates an additional transverse shear component that generates lateral torsion of the cantilever and is the source of TSM signals. On an isotropic surface, the transverse shear signal is zero.
The star-shaped monolayer MoS 2 flake was characterized using AFM (figures 24(c)-(e)). We can observe remarkable contrasts of the star-shaped domains in the TSM image (figure 24(e)) but not in the corresponding FFM (figure 24(d)) and topography (figures 24(b) and (c)) images. FFM fails to distinguish the frictional and elastic anisotropies of the domains probably because of its low sensitivity.
For TSM characterization with a series of clockwise rotations, the periods of 60 • of the shear signals are presented ( figure 24(g)), thereby implying consistency with the hexagonal symmetry of MoS 2 . More precisely, the scatterplot diagrams of shear signal vs. crystallographic orientation for the seven domains (figure 24(h)) exhibit the same distribution regularity, which confirms that the anisotropic shear deformations accurately correspond to crystallographic orientation.
TSM characterization is useful for crystallographic orientation imaging of monolayer samples and characterization of bilayer samples, such as bilayer MoS 2 and graphene [19].  Figure 25(e) shows the FFM images of the bilayer MoS 2 flake. No crystalline differentiated contrast is observed on top and bottom layer grains. The friction signal on the bilayer area is slightly lower than that on the monolayer area. Figure 25(f) shows the TSM images of this bilayer MoS 2 flake. The top (P 21 and P 22 ) and bottom layer grains (P 11 and P 12 ) are clearly resolved in the TSM image. The P 21 and P 11 grains are in a near-pristine 2H stacking. P 22 shows an almost identical shear signal as P 12 due to their pristine 2H stacking with the same but inversion-asymmetry crystallographic orientation. As the vdW interaction between the top and bottom MoS 2 layer is less than the bonding between the bottom MoS 2 layer and substrate, the difference between shear signals (trace and retrace) of the top layer is larger than that in the bottom layer. The top layer grains are stacked over the bottom layer grains. Whereas no shear contrast is observed within the top layer grains. This result demonstrates that the puckering effect occurs primarily on the top layer, and the corresponding puckering-induced shear signal in the bottom grains should be small.
The top layer-dependent crystallographic orientation imaging of 2D materials enables the crystallographic orientation imaging of bilayer films with TSM. Our research will be beneficial in understanding the nanomechanical behaviors of 2D systems and providing a convenient and powerful approach to facilitate nondestructive crystallographic orientation characterization of 2D atomic crystal systems.   [156]. Copyright (2021) American Chemical Society.

Strain-induced rippling and manipulation.
The strain-engineered rippling structures were reported in TMDCs by several groups. Ripples can strongly influence electronic properties by inducing effective magnetic fields and changing local potentials [168][169][170][171][172]. Here, the nanoscale ripples introduced by the in-plane thermal strain are reported. Using the TEC mismatch between the TMDCs and growth substrate, we can apply built-in strain in TMDCs layers, as shown in figure 26(a).  The ripple orientation is illustrated in the inset, and the direction of rippling domain I is horizontal (0 • ). The AFM manipulation was systematically performed as follows: first, repeated (forward and backward) scanning in contact mode was performed with a certain manipulation scan direction and large loading force (∼400 nN) on the entire flake. Next, the TSM image was obtained in the horizontal scan direction with a small loading force (∼0.1 nN) to check the result of the previous manipulation. By incrementally changing the manipulation direction from 0 • to 180 • , three rippling domains with different contrasts were obtained, labeled type I, II, and III. The ripple direction of type I (II and III) domain was further aligned along the y (b and p) zigzag edge of the WS 2 flake using the angle-dependent TSM images. In our experiments, the rippling domain of type I, II and III can only form by the manipulation directions within the angle interval of (−30 • , 30 • ), (30 • , 90 • ) and (90 • , 150 • ), respectively. Figure 27(e) describes the conclusion of manipulations: the manipulation is controlled primarily by the angle interval of the target zigzag-orientated ripple lines direction and manipulation directions. The types II or III rippling domain could be transformed to type I using AFM manipulation when the angle interval between manipulation direction and type I ripple direction range from −30 • to 30 • . Decreasing the angle interval allows successful manipulation using relatively fewer manipulation repeating times and smaller loading force. The direction of the rippling domain is the same if the angle interval between manipulation and pristine ripple direction is less than 30 • .
The AFM manipulation based on the above manipulation principle (figures 27(f) and (g)) has fabricated numerous distinctive rippling domain patterns in the monolayer WS 2 flake. Figure 27(f) shows the procedure of writing a matryoshkalike rippling domain pattern. More complicated patterns were 'painted' using the AFM lithography method ( figure 27(g)).
The single rippling domain was first prepared on the entire WS 2 flake and then formed within the flake using contact AFM scanning several times in the home-modified lithography AFM mode.
The artificial rippling domain patterns were stable in the ambient condition and can be imaged after several weeks. Furthermore, the rippling strongly influences the electronic properties of 2D materials by introducing effective periodical local pseudo-magnetic fields and electrical potentials [168,173]. This strain-engineered ripple and artificially-manipulated rippling domain could further investigate the exotic electron behaviors in the 2D limit.

Nanoscratching using AFM probe.
It is important to explore the fracture mechanics properties of 2D materials deposited on the substrate. To date, experimental research on fracture mechanics has been conducted on 2D materials, such as the nanoindentation test on exfoliated graphene [174] and MoS 2 [175], and nanoscratch test on graphene, MoS 2 , and h-BN [176].
The comparison of Young's moduli and breaking strengths of several materials are summarized in table 4. The strength of monolayer MoS 2 is exceeded only by carbon nanotubes and graphene. For a complete cognition of fracture behavior and to deepen the understanding of different fracture stages, nanoscratch with progressive and constant force is conducted, as in figure 28.
Ye et al [177] applied a progressive series of forces to perform a scratch along a linear routine on the single-layer MoS 2 , as shown in figure 28(a). The normal load increased from 50 to 70 µN in 20 µm. The experiment result is presented in figures 28(b) and (c), which correspond to the optical image and morphology of the sample after scratch test, respectively.
There is no distinct phenomenon in the sample when the normal load is less than 50 µN, and crack modes in the front and middle exhibit distinctive differences. Two novel crack forms can be triggered according to the magnitude of normal loads. Figure 28(d) shows the middle of the fracture where periodical serrated cracks can be found, thereby implying an anisotropic brittle fracture of the single-layer MoS 2 .    Figures 28(e) and (f) correspond to nanoscratch tests along the high-symmetric crystallographic orientations (armchair and zigzag). The size and distance between adjacent semi-circular cracks of the two crystallographic orientations are similar. However, the head shape of cracks along the armchair is flat, whereas that along the zigzag is angular with 120 • at the top. Thus, all the directions of the head cracks are consistent with zigzag orientation in this mode. Such a unique crack form originated from the special local stress distribution around the AFM tip [177].
In summary, semi-circular and periodical zigzag cracks in the fracture were generated according to the different loading force, and both present anisotropy in the generation and propagation progress that is helpful for the design of reliable and protective nanoelectronic devices.

Interfacial charge transfer
Understanding charge generation, transfer, and diffusion between 2D materials and their supporting substrates is important for their potential applications. Lu et al [178] observed significant N-doping in thin MoS 2 films on SiO 2 , dominated by charge traps at the sample − substrate interface by STM. Further, the charge trapping and its effects on the characteristics of the MoS 2 field-effect transistors and photodetection devices were investigated [179,180]. Liu et al [181] claimed to generate anomalously high DC using a sliding Schottky nano-contact in a sliding semiconductor-metal contact (on MoS 2 multilayers/Ag).
Kim et al [182] investigated the triboelectric charges of a graphene/SiO 2 system and found trapped triboelectric charges at the air-SiO 2 interface underneath the graphene that act as ghost floating gates. KPFM was used to measure the SP before and after triboelectrification (figures 29(a)-(c)). Figure 29(b) shows the initial equipotential state of the CVD graphene before rubbing. Figure 29(c) shows the SP after rubbing the central square; the rubbed region has a ∼50 mV higher SP than the unrubbed region. Such potential variation may not be attributed to charges stored in the graphene as electric charges are localized for long times only in insulating materials. Some generated charges tunnel through the monolayer graphene and are locally trapped on the underlying insulator. The trapped charges act as an immaterial (made of charges-only and not of a conductor) bottom floating gate and locally change the polarity and density of charges in the graphene and its work function. The ghost floating gates effectively control the carriers within graphene because of the thin graphene and comparably thin air gap.
∆V TT is defined as the SP taken with reference to the average SP of the unrubbed region. Figure 29(d) shows the averaged ∆V TT along the black dashed line of figure 29(c) taken after 0 and 72 h of rubbing. It is observed that the potential variation is well preserved even after 72 h. Figure 29(e) shows the averaged ∆V TT as a function of time: In practice, in addition to a small term of shorter time constant (τ short ∼ 3 h and 21 min), there is a dominant term with higher initial amplitude and exceptionally long time constant (τ long ∼ 278 h). Herein, τ long is more than two orders of magnitude than that of the standard triboelectrification. The immaterial ghost gates effectively control the charges inside the graphene owing to the thin 2D materials and sub-nanometer insulating air gap.
Recently, the charges transfer and diffusion of the MoS 2 /SiO 2 interface through contact and frictional electrification were investigated systematically in-situ using SKPM and DH-EFM (figure 30). In the contact electrification process, the contact voltages applied by the tip and the positive or negative charges can be transferred to MoS 2 /SiO 2 with the PtSi-coated AFM tip. The micro-process of contact electrification in the MoS 2 /SiO 2 system is shown in figure 30(a). This process can be divided into two steps: first, the biased tip is in contact with the center of MoS 2 and forms an equipotential MoS 2 flake. Second, the generated interface dipole moments act as floating gates. The image charges on SiO 2 can be induced by the charges on MoS 2 flakes by considering the air gap between them. The joint action of charges on MoS 2 and their image charges forms the dipole moments at the MoS 2 /SiO 2 interface. The effect of these interface dipole moments is evident in the DH-EFM images.
The decay process of the charges on MoS 2 /SiO 2 is monitored by the SP evolution, as shown in figures 30(e)-(i). Figure 30(e) shows the SP of the uncharged sample. We simultaneously applied −4 V on MoS 2 and SiO 2 to compare the decay process of charges on MoS 2 /SiO 2 with that on SiO 2 . Figures 30(f)-(i) show a series of SP images taken in the same areas after tip contact sample. The apparent area of charged regions becomes larger and the intensity becomes weaker. Figure 30(b) shows the evolution of SP value with time. The discharge process of MoS 2 /SiO 2 is much longer than that of the SiO 2 surface. The decay time τ is obtained using where ∆ SP is defined as the SP taken with reference to the average SP of the intact SiO 2 region. V 0 is a constant, which equals the initial SP difference ∆ SP0 [182,183]. The decay time of −4 and +4 V charged MoS 2 /SiO 2 is ∼1.98 and 2.42 h, respectively. Overall, the τ for the MoS 2 /SiO 2 areas is one order of magnitude larger than that of the contact electrification SiO 2 surface.
The longer τ on MoS 2 /SiO 2 confirms that the transferred charges are trapped at the MoS 2 /SiO 2 interface because the discharge process is related to the water molecules and atmospheric ions [184]. When the charge is trapped at MoS 2 /SiO 2 interface, the water molecules cannot effectively contribute to deelectrification due to the MoS 2 layers. The discharge process of interfacial charges depends primarily on the charge transfer in SiO 2 bulk. Hence, τ is larger in the MoS 2 /SiO 2 system. Furthermore, τ depends on the sign of charged voltages. Thus, the τ of +4 V charged flake is larger than that of the −4 V charged flake owing to the different charge diffusion coefficients of positive and negative charges in SiO 2 /Si.
Additionally, the characterization of frictional electrification on MoS 2 /SiO 2 was studied (figures 30(j)-(m)). Figure 30(j) shows the topography after rubbing. Figures 30(k) and (l) show the SP image of the sample before and after rubbing. The contrast between R and I areas can be observed even after 80 h. The rubbed areas are negatively charged. The Fermi surface of tip is higher than that of SiO 2 and MoS 2 , thus resulting in electron transfer from tip to sample. The nonuniform SP on SiO 2 substrate may originate from the charged impurities from the air deposited on the substrate.
The corresponding schematic of charge transfer is shown in figure 30(c). The ∆ SP of rubbed MoS 2 area vs. time is plotted in figure 33(d). The τ is ∼14.10 h for the rubbed MoS 2 /SiO 2 . The charge decay time of the MoS 2 /SiO 2 interface is one (or two) order of magnitudes larger than that of the SiO 2 surface. The longer τ confirms the transferred charges are trapped at MoS 2 /SiO 2 interface. Due to the screen effect of MoS 2 layers, water molecules or other atmospheric ions cannot effectively contribute to deelectrification, thereby resulting in a longer τ of interface charges [30].
In addition to the semiconductor/insulator interface, the metal/insulator (M/I) interface was also investigated (figure 31). The charges transfer at M/I interface ( figure 31(a)) is as follows. Electron flows from the material with higher Fermi level E F (smaller work function φ) to the one with lower E F (larger φ). Electron is trapped at the insulator surface when the work function satisfies φ I > φ M . The interfacial dipole moments (indicated by red arrows) collect at the M/I interface owing to the air gap between the insulator (d) Schematic of post-exposed h-BN surface preparation. Post-exposed h-BN surface is the h-BN surface covered by NbS 2 in the growth process but exposed after peeling the NbS 2 using AFM tips. NbS 2 layers framed by the dashed box stand for the stripped part before peeling. Charge distribution at the 1L-NbS 2 /BN interface and post-exposed h-BN surface is also shown. (e) AFM topography after partly peeling NbS 2 , the dash hexagon shows the position of the NbS 2 layer before peeling. (f) SP image of the post-exposed h-BN surface. Inset: SP profiles along the dash lines in (f). (g) Schematic energy band diagram of NbS 2 /BN showing electron transfer from 1L-NbS 2 to h-BN. Scale bars are 1 µm. Reproduced from [92]. © IOP Publishing Ltd. All rights reserved. Reproduced from [185], with permission from Springer Nature. and 2D metal. The NbS 2 /BN heterostructure was obtained using CVD on SiO 2 /Si substrate to identify the M/I interface charge transfer. Figure 31(b) shows the topography of CVD-grown NbS 2 /BN heterostructures. The NbS 2 layers grow on the h-BN substrate in a layered fashion with atomic smoothness. We call the initial growth layer the first layer (1L-NbS 2 ) and the following layers the second (2L-) and third (3L-) NbS 2 layer [92]. Figure 31(c) shows that the 1L-NbS 2 SP is larger than that of h-BN and 2L-NbS 2 by ∼130 and 10 mV, respectively. The SP difference between 1L-NbS 2 and h-BN indicates charge transfer between NbS 2 and h-BN interface, as demonstrated in figure 31(g).
An AFM probe was used as a manipulator to peel part of NbS 2 flakes to obtain post-exposed h-BN surface and 'see' the interfacial charge transfer ( figure 31(d)). A vertical knocking was applied on the NbS 2 layers with AFM tip in AC (tapping) mode. The NbS 2 layer is broken into several small lamellas with sufficiently strong vertical tapping, and the postexposed h-BN surface is obtained, as shown in figure 31(e). Figure 31(f) shows the post-exposed h-BN surface SP as ∼236 mV, lower than that of pristine h-BN. The negative charges accumulate at the post-exposed h-BN surface considering the neutral pristine h-BN surface.
Ding et al [185] used SKPM to characterize the SP of MoS 2 /PbI 2 heterostructures. In figure 32, a significant change in SP in the interlayer between MoS 2 and PbI 2 suggest interface charge transfer and electron transfer from MoS 2 to PbI 2 .
To summarize, SKPM and advanced EFM aid the study of interfacial electrons transfer and diffusion. The material properties can be adjusted by interface engineering and motivate further directive studies on the heterointerface of heterostructures.

Interface intercalation
Interfacial engineering, such as molecule or ion intercalations, can modify properties and optimize heterostructures and their device performances. For example, Bediako et al reported that lithium-ion intercalation at individual atomic interfaces motivates using vdW heteroepitaxy to realize new engineered functional interfaces for energy conversion and storage by manipulating the ion storage modes and 'job-sharing' characteristics of hybrid electrodes [186].
Moreover, the intercalation of atoms and molecules has been studied for h-BN monolayers [187,188] and other 2D materials, such as graphene [189], to modify the materials' properties. Yamasue et al synthesized hydrogen-intercalated graphene on a 4H-SiC(0001) substrate. Hydrogen intercalation at the interface eliminates covalent bonds and the original quasi-(6 × 6) corrugation, which indicates the conversion of the buffer layer into a second graphene layer by terminating Si bonds at the interface.
As water is an anomalous liquid, the water intercalation at the interface of heterostructures has attracted considerable attention. Lee et al [190] reported that H 2 O intercalating between graphene and mica increases the friction between the tip and substrate depending on the water and graphene layer thickness. Figure 33 illustrates the water intercalated singlelayer graphene (SLG) and bilayer graphene (BLG) on mica. The topography and friction images in figures 33(c) and (d) show the regions of mica, SLG + zero water layers (0 W), SLG + first water layers (1 W), and SLG + second water layers (2 W). The brighter contrast in figures 33(b) and (d) suggests high friction. The FFM images show increased friction by the water layer intercalation between graphene and mica. The friction enhancement in BLG is less visible than that in SLG but is higher than that on the BLG + 0 W on mica, as shown in figures 33(d) and (f). Figure 33(g) shows a schematic of varied friction over graphene layers with and without intercalated water using an AFM tip. Additionally, they reported that D 2 O intercalation decreases friction at the H 2 O-intercalated graphene on mica because the low rate of frictional kinetic energy dissipation at the interface affects the phonon contribution.
Hong et al [191] investigated the effects of the trapped interfacial ice-like water layer on the charge transfer between graphene and SiO 2 /Si substrate by recording the SP changes induced by partial removal of the water layer upon in-situ heating. Figures 34(a)-(d) illustrate the in-situ AFM height images of the graphene-interfacial water layer-SiO 2 /Si sandwich system at different temperatures. No noticeable change is observed below 40 • C ( figure 34(b)), whereas the interfacial water layer is partially removed at 80 • C (figure 34(c)). A further increase in temperature removes the interfacial water layer ( figure 34(d)). Interfacial water molecules fail to penetrate the impermeable graphene sheet, allowing the water molecules at the graphene edge to easily escape upon heating ( figure 34(e)). The same height of 0.37 ± 0.02 nm is measured for the interfacial water layer (figure 34(g)), which indicates an ice-like single water layer underneath the graphene. The contrast in the SP map between the 1L 1WL and 1L 0WL graphene regions shows the effects of the interfacial water layer on graphene SP ( figure 34(f)). The work function of graphene is calculated by φ G = φ tip − eV CPD , where φ tip and V CPD are the work functions of the AFM tip and contact potential difference between the AFM tip and graphene measured by SKPM. The SKPM mapping shows electronically modified graphene by the ice-like water layer as the electron density transfers from graphene to the water layer, resulting in graphene hole-doping.  [190], with permission from Springer Nature.
The water layers confined between hydrophobic/hydrophobic interfaces under ambient conditions have been studied by our group [73]. A special folding few-layers graphene film was elaborately prepared (figure 35(a)) to form a hydrophobic/hydrophobic interface. We precool the sample with constant temperature and humidity equipment for several hours before placing the sample in a warm and high humidity condition to realize H 2 O intercalation at the whole interface. Because the MF-AFM amplitude (A MF ) can reflect the local mechanical properties of the sample and has subsurface detection capability ( figure 35(b)), the internal interfacial intercalation water phases were further visualized with two MF-AFM modes: multi-harmonic and dual AC AFM.
The subsurface information obtained through nanomechanical coupling between the tip and sample helps distinguish between the two types of water phases via the higher harmonic A 6th detection ( figure 35(d)). Essentially, the elastic modulus increases for higher A 6th values. The corresponding histogram in figure 35(e) is extracted from figure 35(d). The larger (smaller) elastic modulus areas, corresponding to the bright (dark) areas in A 6th , represent the solid (liquid)-like phase of the intercalation water layers. Reproduced from [191] with permission from the Royal Society of Chemistry. Figure 35(g) shows the A 2 images of the water intercalated sample and the corresponding histogram distribution in figure 35(h). Two kinds of water structures coexist at the folding graphene interface. The higher A 2 (brighter contrast) indicates the strong repulsive (rigid) and less dissipated tip-sample forces, which reflects the solid-like phase. In contrast, the dark areas with lower A 2 reflect a liquid-like phase considering their soft and dissipative nature.
The experimental data of A 6th and A 2 demonstrate coexisting solid-and liquid-like water phases at the interface. The internal interfacial structures affect the solid and liquid layers, which confirm intercalated water layers at the interface. The distribution rates of the two phases are almost the same in the A 6th image. However, in the A 2 image, the distribution rates are ∼60% for solid-and ∼40% for liquid-like phases. This discrepancy arises because the tip-sample force is weaker in multi-harmonic mode (passive MF-AFM mode) than in dual AC mode (active MF-AFM mode).
Depending on the density, the liquid phase is divided into low-density liquid (LDL) and high-density liquid (HDL). Considering the fragility indices, LDL (m = 14) is a 'superstrong' liquid whereas HDL (m = 20-25) is relatively weak. Combing the MF-AFM data analysis and above discussion, phases of liquid water were determined using elastic modulus. LDL has a larger elastic modulus (solid-like phase) than the HDL (liquid-like phase). The water patches and films appear as one/two puckered bilayers of ice-I h in the section lines of the same area (figure 36(e)) [92,192]. The 2L-NbS 2 height remains unchanged before and after water adsorption. We conclude that except at the intercalation at NbS 2 /BN interface, the interface between 1L-and 2L-NbS 2 is free from water owing to the absence of charge transfer and dipole moments between the two layers.
The electrical properties of water-intercalated samples were further investigated using DH-EFM and sMIM. Compared to the mobile charge carrier (MCC) images of the pristine sample (figure 36(h)), 1L-NbS 2 is 'invisible' and has an identical signal with the insulating H 2 O/BN surface after water intercalation, as highlighted by the dotted lines in figure 36(i). This indicates the MIT of the 1L-NbS 2 layer when intercalating the H 2 O layer at the 1L-NbS 2 /BN heterointerface. The 2L-NbS 2 layer can be observed in MCC images, which indicates the preserved metallic state of 2L-NbS 2 by the buffering effect of the 1L-NbS 2 layer. Additionally, figure 36(j) shows that the SP of water-intercalated NbS 2 /BN heterostructures is ∼120 mV lower than that at  Dollekamp et al [193] tuned the friction of graphene by changing the size and vibrational modes of the intercalated molecules. They replaced water with larger alcohol molecules of different vibrational modes, e.g. stretching the C-C and C-O bonds. The extra vibrational modes contribute to the higher friction in methanol, ethanol, and 2-propanol compared to an intercalated double water layer, which scales with their size. In addition, the single-layer graphene has higher friction than bilayer and trilayer graphene.
To summarize, interfacial intercalation can tune the local physical properties of the sample by changing the environmental conditions and intercalated interface molecules. These results suggest manipulating the electronic/mechanical properties of 2D materials for diverse applications.

Interfacial thermology
Thermal measurements require nanoscale spatial resolution and high sensitivity. SThM can successfully determine the thermal conductivity, k, of samples [194]. For example, the in-plane k of a residue-free suspended graphene bridge was measured using NP-SThM with spatial resolution reaching 50 nm [144]. The temperature profile of suspended graphene bridges was measured and fitted with the 1D heat transfer equation to obtain k. Hwang and Kwon [195] used NP-SThM to obtain the in-plane k of suspended graphene disks of radius 50-3680 nm.
In addition, the thermal transport of heterostructures could be studied by SThM. The heater-sample thermal contact resistance (R X ) of graphene/MoS 2 heterostructure is given in figure 37(a). The sample has areas of bare SiO 2 /Si (SOS), graphene on SiO 2 (GS), MoS 2 on SiO 2 (MS), and MoS 2 on GS (MGS). The contrast in the R X image demonstrates the following trend from low to high thermal resistance areas. R X (GS) < R X (SOS) < R X (MGS) < R X (MS). This is because the interfacial thermal resistance plays a role in thermal transport in addition to the layer thermal conductivity [196]. The MS interface exhibited an R X between 4 × 10 −8 and 2.27 × 10 −6 m 2 KW −1 [197,198], which is higher than that of GS (between 5.6 × 10 −9 and 2 × 10 −8 m 2 KW −1 ) [199]. Thus, the thermal resistance of GS is increased by one order of magnitude by adding MoS 2 top layer, thereby providing pathways for increasing the efficiency of thermoelectric applications using vdW materials.
Additionally, the spatial distribution of temperature rise within the lateral heterostructure was studied [201]. The optical and AFM images (figure 37(c)) of a MoS 2 -WS 2 lateral heterostructure confirm a single monolayer sample. The compositional heterogeneity of this device is illustrated by the Raman map in figure 37(d). SThM maps the temperature rise spatial distribution within the monolayer TMDC devices during a high electrical power dissipation through a lateral interface (figures 37(e) and (f)). The results directly demonstrate that lateral heterojunctions between MoS 2 and WS 2 fail to impact the heat dissipation distribution, whereas GBs of MoS 2 considerably localize heating in the device.
Special physical phenomena could be observed when applying SThM. Xu et al [53] observed a novel mechanicalthermal coupling effect in monolayer/bilayer MoS 2 , and WS 2 films; essentially, puckering deformation can enhance interfacial TR (figure 38). Figure 38(a) shows the FFM images of monolayer/bilayer MoS 2 with increasing normal forces. Figure 38(b) shows the corresponding friction signals vs. normal forces. The AFM tip's normal force applied to the sample surface can controllably modulate the magnitude of puckering deformations, which facilitates corresponding frictional and thermal response detection. Corresponding to the FFM images, the in-situ dynamic TR images of the monolayer/bilayer MoS 2 sample are shown in figure 38(c). The lighter contrast stands for a larger TR. Under normal forces, the dynamic TRs of the monolayer/bilayer MoS 2 sample is enhanced. The corresponding relative dynamic/static TR vs. normal force (relative to the substrate with no puckering effect) are drawn in figure 38(d). Compared with the static TRs, the dynamic  TRs of 1L-and 2L-MoS 2 rapidly increases with increasing normal forces. The results suggest that puckering deformations modulated by normal forces contribute to increasing dynamic TRs.
The principle of static and dynamic TR is illustrated in figure 38(e). After contact, the thermal contact area increases with increasing normal force owing to the elastic compression, which reduces the static and dynamic TR. However, as the dynamic contact scan can induce puckering deformation, a slight interface gap exists between the top layer and substrate, which impacts longitudinal heat transfer with additional dynamic interfacial TR. Therefore, dynamic TR is higher than static TR (with no puckering deformation) under the same normal force.
Furthermore, the crystallographic orientation-dependent anisotropy of the puckering effect in the atomically thin 2D crystal was demonstrated using SThM (figures 38(f) and (g)).
The puckering deformation causes in-plane stress redistribution, thus breaking the original in-plane stiffness. These findings are significant to optimize the nanoscale tribological/thermal design and dynamic mechanical-thermal management of 2D materials in nanoelectronics.

Conclusion and perspectives
This review highlighted the principles of functional AFM methods and their applications in 2D materials. AFM-based physical characterization methods require the researchers to have clear natural science concepts and electronic circuit knowledge. In TSM, the enhanced sensitivity to elastic deformation properties of crystals can reveal the crystallographic orientation of samples. The nonlinear sample information can be obtained using multifrequency force microscopy. KPFM can detect the sample work function, whereas EFM can explore the carrier concentration, electrostatic gradient, and MCC. In sMIM, the dielectric constant, doping density, conductivity, and permittivity of the sample surface and subsurface can be obtained (semi-) quantitatively when combined with FEA analysis.
The AFM modes require further development and improvement in the applied research on the surface/interface of 2D materials. For example, (i) the single physical property measurement method has limitations because the electrical properties are coupled with other properties. Therefore, a method that can explore the mechanical-thermal-electricaloptical coupling is needed. (ii) The functional AFM technology under vacuum conditions needs further improvement.
Obtaining quantitative data is difficult as AFM detection in the atmospheric environment is inevitably affected by the water molecules. In addition, the functional AFM technology involves complex signal transmission of microwave or light waves, which results in unsatisfactory manufacturing costs and stability of functional AFM technology under vacuum conditions. (iii) Improving the application of 2D materials and optimizing device performance requires understanding the properties under work conditions. A method that can measure physical properties under work conditions is crucial for future research. With the development of AFM technology and 2D material preparation, we expect increased achievements in this field.