# Fileset

[draft013.pdf](https://mdr.nims.go.jp/filesets/e29562d0-99a1-4b86-8588-eb79b57d0948/download)

## Creator

[Nobuyuki Ishida](https://orcid.org/0000-0003-0161-0583)

## Rights

This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in Nobuyuki Ishida; An open-source framework for quantitative electrostatic simulations in Kelvin probe force microscopy. J. Appl. Phys. 14 July 2026; 140 (2): 024301 and may be found at https://doi.org/10.1063/5.0336722.[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[An open-source framework for quantitative electrostatic simulations in Kelvin probe force microscopy](https://mdr.nims.go.jp/datasets/1a45a006-52b0-4c09-91c7-9661b1333176)

## Fulltext

An Open-Source Framework for Quantitative Electrostatic Simulations in KelvinProbe Force MicroscopyNobuyuki Ishida1, ∗1National Institute for Materials Science, Tsukuba, Ibaraki 305-0047, Japan(Dated: July 18, 2026)We developed an open-source numerical code for calculating electrostatic interactions betweena hyperbolic metallic tip and a conducting surface under Kelvin probe force microscopy (KPFM)-relevant conditions. The code solves the electrostatic potential in the tip–sample vacuum gapunder an applied bias by imposing Dirichlet boundary conditions on the tip and sample surfaces,and evaluates electrostatic force observables from the calculated potential. In addition to laterallyhomogeneous samples, the full three-dimensional implementation treats spatially inhomogeneoussurfaces and reproduces contact potential difference profiles. Quantitative comparison with exper-iments is achieved without arbitrary scale factors or other ad hoc proportionality constants. Themethod provides a practical framework for experimental-parameter extraction, quantitative KPFMdata analysis, and future investigations of contrast formation and other measurement aspects inKPFM.I. INTRODUCTIONKelvin probe force microscopy (KPFM) is a power-ful technique for visualizing nanoscale surface-potentialdistributions [1–3] across a wide range of materials un-der realistic device-operating conditions [4–9]. However,because the KPFM signal originates from electrostaticinteractions, it inevitably includes long-range contribu-tions. As a result, the measured potential does not neces-sarily correspond directly to the local value immediatelybeneath the tip apex, but rather reflects a weighted aver-age of the surrounding electrostatic environment. Conse-quently, a consistent framework that relates probe geom-etry and measurement conditions to the resulting electro-static field and force is essential not only for the quan-titative interpretation of KPFM data, but also for thecontinued improvement of KPFM methodology.To date, numerous studies have evaluated electro-static tip–sample interactions using numerical and semi-analytical approaches to discuss the mechanisms of signalformation in electrostatic force microscopy (EFM) andKPFM [10–14]. The nonlocal response of KPFM has alsomotivated complementary inverse-analysis approaches,including point-spread-function/deconvolution methods,Green’s-function methods, and finite-element/KPFM-hybrid approaches for reconstructing surface-potential orcharge-density distributions from measured KPFM im-ages [15–20]. Most of the forward-modelling studies fo-cused on cantilever-based force sensors and examined, forexample, the differences in contrast formation betweenamplitude-modulation (AM) and frequency-modulation(FM) EFM/KPFM [21–23], as well as how different partsof the probe assembly, including the tip apex, the shank,and the cantilever, contribute to the electrostatic inter-action [11–13, 23, 24]. By contrast, relatively few studieshave explicitly considered AFM configurations employ-∗ ishida.nobuyuki@nims.go.jping a sharp metallic tip mounted on a self-sensing forcesensor and operated under ultrahigh-vacuum (UHV) con-ditions, where very small tip radii and tip–sample sepa-rations are typically assumed for high-spatial-resolutionmeasurements [25, 26].Representative studies by Belaidi et al. [27] and Shenet al. [22] examined how the probe–sample geometry af-fects the apparent spatial resolution, showing that, in theregime s � R, it scales as (Rs)1/2, where R and s denotethe tip radius of curvature and the tip–sample separation,respectively. However, these studies were primarily con-cerned with forward analyses of contrast formation andresolution limits, rather than with rigorous comparison toexperimental observables or the establishment of a practi-cal framework for quantitative data analysis. As a result,a general approach for extracting physically meaningfulprobe parameters from measured data has not yet beenfully established. One likely reason is that the tip radius,opening angle, and tip–sample separation can be mutu-ally correlated, making it difficult to impose physicallymeaningful constraints on the fitting parameters. More-over, the accuracy of existing models and the precision ofexperimental calibration have not always been sufficientto support fully quantitative fitting capable of reliablyidentifying physical parameters from experimental data.In this study, we developed a numerical software pack-age, METALTIP, to compute the electrostatic responseof a hyperbolic metallic tip above a conducting surface.For user-defined parameters such as the tip radius, open-ing angle, and tip–sample separation, the code solves theelectrostatic potential in the tip–sample vacuum gap un-der an applied bias and computes the resulting electro-static force. This enables direct simulation of the bias(U) dependences of the electrostatic force, F (U), andfrequency shift, ∆f(U). The full three-dimensional (3D)implementation further treats laterally varying surfacepotentials and reproduces contact potential difference(CPD) profiles across inhomogeneous interfaces. We val-idated the numerical accuracy by benchmarking againstestablished theoretical models and demonstrated quanti-2tative fitting of experimental data without ad hoc tuningparameters. This quantitative accuracy enables in situextraction of key experimental parameters, including thetip–sample separation and the effective tip radius. Thesecapabilities establish METALTIP as a practical frame-work for the quantitative interpretation of KPFM dataand for the development of new measurement and anal-ysis approaches.II. NUMERICAL METHODSTo quantitatively calculate the electrostatic interactionbetween a metallic tip and a conducting sample surface,METALTIP builds on the potential-solver framework ofSEMITIP developed by Feenstra [28–33] and largely fol-lows its numerical implementation. SEMITIP was origi-nally developed to model scanning tunneling microscopy(STM) measurements on semiconductor surfaces by cal-culating the tunneling current, which requires solvingthe electrostatic potential and charge distribution insidesemiconductors. In contrast, METALTIP is designed foratomic force microscopy (AFM)/KPFM measurements,where the relevant observable is the electrostatic forceacting on a metallic tip. The code addresses the elec-trostatic boundary-value problem of an ideal metallic tipabove a sample surface whose potential is specified asa Dirichlet boundary condition. Thus, the electrostaticpotential is solved only in the tip–sample vacuum gap.This formulation does not describe electrostatic re-sponses inside the sample, such as finite screening, dielec-tric polarization, charge redistribution, or band bending.Therefore, it is not directly applicable to general semi-conductor surfaces where electric-field penetration intothe sample plays an essential role. However, semiconduc-tor surfaces that behave approximately as equipotentialboundaries owing to strong screening or a high densityof surface states may still be treated within the sameeffective-boundary framework.In the following, we outline only the essential computa-tional workflow and key equations; full numerical detailsand implementation-specific procedures are provided inSupplementary Note I.In the vacuum region, the electrostatic potential φsatisfies the Laplace equation, ∇2φ = 0. In the two-dimensional (2D) case, corresponding to a laterally ho-mogeneous sample, axial symmetry is assumed, and thepotential is written as φ(r, z) in cylindrical coordinates,where r is the radial coordinate and z is the surface-normal coordinate. The (r, z) domain is discretized us-ing a generalized prolate spheroidal coordinate system(ξ, η), which enables a compact grid that captures boththe apex curvature and the shank opening. The samplesurface (z = 0) corresponds to η = 0, and the tip surfaceis aligned with η = ηT , where ηT is determined by thefull shank opening angle α (see Supplementary Fig. S1).After solving the Laplace equation numerically to ob-tain the electrostatic potential φ, we evaluate the elec-trostatic force from the electrostatic pressure on the tipand sample surfaces. The surface pressure is obtainedfrom the electric field at the conductor surfaces using theMaxwell-stress expression. In vacuum, the pressure isgiven by [34]p =12ε0E2n, (1)where ε0 is the vacuum permittivity and En is thesurface-normal component of the electric field. Becausethe electric field at a conductor surface is normal to thesurface, En can be evaluated from E = −∇φ at theboundary.The vertical force component is then obtained by inte-grating the pressure over the surface,Fz =∫∫p nz dS, (2)where nz is the z component of the unit normal vectorand dS is the surface area element. In practice, we eval-uate Fz from the pressure on the planar sample surface,which is numerically more stable, especially for sharptips. We verified that the sample-side and tip-side forceevaluations agree to within a few percent after mesh re-finement; details of the implementation and convergencetests are provided in Supplementary Note I.The 3D solver introduces the azimuthal angle θ andsolves φ(ξ, η, θ) with periodic boundary conditions in θ.An arbitrary surface-potential distribution can be pre-scribed on the sample surface as φ(ξ, η = 0, θ) = Vs(ξ, θ),enabling laterally varying potentials (and correspondingCPD variations), such as interface steps and patch-likeinhomogeneities. In contrast to the laterally uniform 2Dcase, such spatially varying surface potentials generallyproduce not only a normal electric-field component butalso lateral field components on the sample surface. Ac-cordingly, the vertical force component in the 3D calcula-tions is evaluated from the Maxwell stress using the fullsurface traction, so that the contributions of both thenormal and lateral electric-field components are properlytaken into account. The explicit expression used for thisevaluation is given in Supplementary Note I.III. EXPERIMENTAL METHODSTo obtain experimental data suitable for direct com-parison with the simulations, we prepared a Cu(111)surface partially covered with an ultrathin NaCl over-layer [35, 36]. The coexistence of bare Cu(111) terracesand NaCl islands on the same substrate provides a well-defined lateral boundary across which the surface poten-tial changes abruptly [36]. Details of the sample prepa-ration are provided in Supplementary Note II.All KPFM measurements were performed at 78 K un-der ultrahigh-vacuum (UHV) conditions (base pressure< 1× 10−8 Pa) using a low-temperature scanning probe3microscopy (SPM) system (Unisoku USM-1400). qPlussensors [37, 38] equipped with electrochemically etchedtungsten (W) tips (P-100WS, Unisoku) were used asprobes. The resonance frequencies of the sensors rangedfrom 27 to 31 kHz, and the quality factors ranged from18,000 to 25,000 at 78 K. The W tips were subjected toAr+-ion sputtering for 20 min to remove the oxide layerand sharpen the tip apex. Tip–sample forces were mea-sured in FM mode [39] after calibration of the oscillationamplitude [40]. The tip–sample separation was regulatedin STM mode, that is, using the cycle-averaged tunnelingcurrent as the feedback signal. Bias voltages are definedas the sample voltage with respect to the tip. CPD sig-nals were acquired by Kelvin probe force spectroscopy(KPFS) [41–43], in which ∆f–U spectra were recordedat each measurement position and fitted to a parabola.IV. RESULTS AND DISCUSSIONA. Potential and electric-field distributions nearthe tip apexFigure 1(a,b) shows the calculated electrostatic po-tential distribution and the corresponding electric-field-magnitude map for a representative tip–sample geometrywith a tip radius R = 10 nm, a full shank opening an-gle α = 10◦, a tip–sample separation s = 1.5 nm, and abias voltage U = 1 V. The computational domain spans∼20 µm in the radial direction and ∼18 µm in the z direc-tion, ensuring that the outer boundaries are sufficientlyfar from the tip–sample junction. To clearly visualize thestrong-field region near the apex, only a magnified viewaround the tip apex is presented in Fig. 1. The elec-tric field was calculated from the potential distributionas E = −∇φ, and the color scale in Fig. 1(b) representsits magnitude, |E|. The |E| map indicates pronouncedlocalization of the electric field in the immediate vicinityof the apex. In the following section, this localization isanalyzed more quantitatively through the radial distri-bution of the electrostatic pressure.B. Radial distribution of electrostatic pressure andvalidation against the sphere–plane modelUsing the same numerical data as in Fig. 1, we evalu-ate the radial dependence of the electrostatic pressure onthe tip and sample surfaces, as shown in Fig. 2. In bothevaluations, the pressure decreases steeply with increas-ing r over a length scale comparable to the tip radiusR (10 nm), indicating that the dominant contributionarises from the near-apex region. A weak long-rangetail, which is not readily visible on a linear scale butbecomes evident on a logarithmic plot (SupplementaryFig. S2), nevertheless contributes to the total force be-cause the surface integral is area-weighted (∝ r dr). Ac-cordingly, a sufficiently large radial integration range isFIG. 1. Calculated maps of (a) the electrostatic potential φand (b) the electric-field magnitude |E| in the vacuum regionnear the tip apex for a representative tip–sample geometry(R = 10 nm, α = 10◦, and s = 1.5 nm). The tip is held atU = 1 V and the sample surface is grounded (0 V).� � �� �� �� ���������������������������������������������������������� ������� � ����������� ����������������������������������������������������FIG. 2. Radial profiles of the electrostatic pressure on the tipand sample surfaces, calculated for the same geometry andbias conditions as in Fig. 1. Solid circles show the MET-ALTIP results, evaluated on the tip surface (blue) and thesample surface (red). Dashed lines indicate the analyticalsphere–plane solution.required for the tip-side and sample-side force evalua-tions to achieve quantitative agreement (typically within∼5%). This agreement serves as an action–reaction con-sistency check.The convergence behavior of the cumulative verticalelectrostatic force also provides a basis for estimatingthe influence of probe parts not explicitly included in4the present model, such as the macroscopic sensor bodyand the transition region between the modeled hyper-bolic tip and the sensor body. In the present calculations,the hyperbolic metallic tip is included within a compu-tational domain extending in the z direction to approxi-mately 18 µm from the apex, where the cumulative ver-tical electrostatic force is already close to saturation, asshown in Supplementary Fig. S3. In typical self-sensingforce sensor configurations, the metallic parts outside thismodeled region are located much farther from the tip–sample junction, typically several hundred micrometersfrom the apex. Their additional contribution beyond theincluded computational domain is therefore expected tobe small. Moreover, for the FM-KPFM response primar-ily considered here, such distant contributions are furthersuppressed because the frequency shift is governed by theforce gradient rather than by the absolute force, and cantherefore be regarded as negligible.When the tip–sample separation is sufficiently smallcompared to the tip radius (s � R), the near-apex elec-trostatic response is expected to be governed primarily bylocal curvature and narrow-gap geometry. Accordingly,in the immediate vicinity of the apex, the pressure pro-file should approach that of the canonical sphere–planemodel [12], for which an analytical solution is availablevia the method of images [44].To test this expectation, we compared the pressuredistributions calculated by METALTIP with the corre-sponding analytical sphere–plane solution. For this com-parison, the sphere radius and sphere–plane gap werechosen to match the apex curvature and tip–sampleseparation of the hyperbolic tip geometry. As shownin Fig. 2, the METALTIP results closely follow thesphere–plane profiles on both the tip and sample sidesfor r . 3 nm, with deviations below 1.0%. This excel-lent agreement in the near-apex region supports the va-lidity of the electric-field evaluation and the subsequentMaxwell-stress-based force calculation. In contrast, thelogarithmic-scale plot in Supplementary Fig. S4 showsthat, at larger r, the sample-side pressure increasinglydeviates from the sphere–plane prediction. This discrep-ancy is naturally attributed to contributions from the tipshank.C. Bias dependence of electrostatic forceFigure 3 shows the bias dependence of the electrostaticforce (F–U curve) calculated with the same geometric pa-rameters as in Fig. 1. The calculated F (U) characteristicis well described by a parabola,F (U) = a (U − U0)2 + F0, (3)where U0 and F0 denote the bias and force offsets, respec-tively. In the present calculation, the tip–sample CPDis set to zero, and only the electrostatic force above alaterally homogeneous surface is considered; accordingly,���� ���� ��� ��� ��� �����������������������������������������������������������������������
���������FIG. 3. Bias dependence of the electrostatic force (F–Ucurve) calculated for the same tip–sample geometry as inFig. 1. Solid circles denote the calculated values, and thedashed line shows a parabolic fit using Eq. (3).U0 = 0 and F0 = 0. This parabolic behavior confirms theexpected dominant U2 scaling of the electrostatic force.As the tip radius R increases, the parabolic curvatureparameter a increases monotonically and follows an ap-proximately linear dependence on R over the range ex-amined here (2–100 nm), as quantified in SupplementaryFig. S5. This scaling is consistent with the near-field limit(s � R) of the semi-analytical sphere-plus-cone model ofHudlet et al. [12], which predicts that a ' πε0(R/s). Theobserved a(R) dependence therefore further supports thevalidity of the force calculation.D. Conversion from force to frequency shiftIn FM-AFM and KPFM, the experimentally accessibleobservable is the resonance-frequency shift of the oscil-lating force sensor rather than the force itself. The zcomponent of the tip–sample force, Fz, can be convertedinto the resonance-frequency shift ∆f using the standardexpression [45, 46],∆f(z0) = − f0kA2〈Fz(z0 + q(t)) q(t)〉 , (4)where A, f0, and k are the oscillation amplitude, reso-nance frequency, and spring constant of the force sen-sor, respectively, and z0 denotes the mean tip position.The tip deflection during one oscillation cycle is given byq(t) = A cos(2πf0t), and 〈· · · 〉 denotes the average overone oscillation period. In the implementation, this cy-cle average is evaluated by sampling Fz(z) at discrete tippositions along the oscillation trajectory and calculatingthe corresponding weighted average.One motivation for developing METALTIP is toquantitatively simulate KPFM under experimental con-ditions aimed at high spatial resolution. In this regime,5KPFM is often performed in UHV using sharp metal-lic tips mounted on self-sensing force sensors [37, 38].In such measurements, the tip–sample separation is fre-quently regulated in STM mode using the averaged tun-neling current (I) as the feedback signal. Accordingly,in simulations of FM-mode experiments, the tip–sampleseparation s determined by STM feedback under staticconditions (i.e., with the tip oscillation switched off) wasused as an input parameter. For electrostatic calcula-tions under oscillating conditions, s was converted intothe mean tip–sample distance 〈s〉 by taking into accountthe oscillation amplitude and the tunneling-current de-cay constant κ. Assuming I ∝ exp(−2κs), the conversionis given by 〈s〉 = s + (2κ)−1 ln[I0(2κA)], where I0 is themodified Bessel function of the first kind (see Supplemen-tary Note I for the detailed derivation). In the presentanalysis, we used κ = 11.2 nm−1, as experimentally de-termined from the distance dependence of the tunnelingcurrent. Hereafter, s refers to the nominal input separa-tion defined by the static STM setpoint and should bedistinguished from 〈s〉.Supplementary Fig. S6 shows the ∆f–U curve calcu-lated using the same geometric parameters as those inFig. 3. The parameters f0 = 30 kHz and A = 500 pmwere chosen to represent typical experimental conditions,while k was set to the nominal design value of the qPlussensor (1852 N/m for our sensors). As for the F–Ucurves, the calculated ∆f–U curve is well described bya parabola. The resulting frequency-shift values also fallwithin the range typically observed experimentally.E. 3D simulation of CPD profiles for a laterallyinhomogeneous surfaceThe 3D implementation of METALTIP allows arbi-trary surface-potential distributions to be prescribed andcan therefore be applied to laterally inhomogeneous sam-ples. As the simplest demonstration, we considered a sur-face consisting of two regions with different surface poten-tials, separated by a straight boundary passing throughthe origin. The potential difference between the two re-gions was set to 0.5 V. To avoid an unphysical disconti-nuity on the numerical grid, the potential contrast wasintroduced not as an ideal step function, but as a smoothinterpolation described by a hyperbolic tangent functionwith the interfacial width parameter winterface = 1.0 nm(see Supplementary Fig. S7). The tip–sample geometryparameters were set to R = 5.0 nm, s = 1.0 nm, andα = 10◦.The interfacial width parameter winterface was intro-duced as a numerical regularization parameter, ratherthan as a physical length representing the intrinsic transi-tion width of the surface potential at a particular bound-ary. A mathematically discontinuous surface-potentialstep produces an artificial singularity on the numericalgrid, whereas the actual electrostatic potential near aninterface is expected to vary over a finite atomic-scale��� ��� � �� �������������������������� �����
��������������������FIG. 4. Simulated CPD profiles across a lateral surface-potential boundary for a two-region model surface with a po-tential difference of 0.5 V. The calculations were performedwith R = 5.0 nm, s = 1.0 nm, and α = 10◦. The CPD pro-files were extracted by parabolic fitting of simulated F–U and∆f–U spectra at each lateral tip position, corresponding toAM- and FM-KPFM, respectively. Here, x = 0 correspondsto the interface.length because of charge redistribution and screening.We therefore used winterface = 1.0 nm as a practicalsmoothing width: it is larger than the expected atomic-scale screening length and the local grid spacing nearthe interface, but still small compared with the several-nanometer CPD broadening caused by the finite tip size.Using this prescribed surface-potential distribution,the tip was scanned laterally across the interface, and ateach position the corresponding F–U and ∆f–U curveswere calculated. CPD profiles were then obtained byparabolic fitting of the simulated spectra at each lateralposition. The profiles extracted from the F–U and ∆f–Ucurves correspond to the responses measured in AM- andFM-KPFM, respectively. Figure 4 presents the resultingCPD profiles. The calculated profiles clearly show thatthe FM-mode signal yields a sharper transition across theinterface. This reflects the stronger weighting of the tip-apex contribution in the frequency-shift signal [23] andis consistent with the higher spatial resolution generallyachieved in FM-KPFM [21].F. Fitting to experimental data and parameterextractionFitting the experimental data serves three purposes:(i) externally validating the numerical framework un-der realistic operating conditions, (ii) enabling quanti-tative interpretation of the experimental data, and (iii)extracting key parameters that are otherwise difficult toaccess, such as the electrically effective tip radius andtip–sample separation. In this context, the electricallyeffective tip radius denotes the radius parameter of theideal hyperbolic metallic tip that reproduces the mea-6sured electrostatic response, including the CPD-profilebroadening and the ∆f–U curvature. This value maydiffer from the geometric tip radius estimated from SEMimages, because the actual tip shape can deviate fromthe ideal hyperbolic geometry and the electrostatic in-teraction reflects the conducting surface contour over afinite region around the apex rather than only the localapex curvature.The third objective, parameter extraction, is not gen-erally straightforward because electrostatic observablesoften depend on multiple geometric parameters in astrongly correlated manner. Fitting a single datasetalone therefore does not necessarily yield a unique pa-rameter set. Here, we address this identifiability issue bycombining two complementary measurements acquiredwith the same tip: a CPD line profile across a hetero-geneous interface and a ∆f–U spectrum measured on alaterally homogeneous metallic surface.1. SEM-based evaluation of tip geometryIn this study, we used commercially available electro-chemically etched W tips (Unisoku). To characterize thestatistical variation in tip geometry, we extracted the ge-ometric parameters (R and α) from SEM images. Thetip radius spanned approximately 5–30 nm, and morethan half of the tips exhibited R < 10 nm. The fullshank opening angle ranged from 5◦ to 20◦, with mosttips showing values below 10◦.For many tips, the outer profile could be fitted wellby a single hyperbolic curve over an axial range of sev-eral hundred nanometers from the apex, as exemplifiedin Supplementary Fig. S8(a). In contrast, some SEMprofiles could not be reproduced over the entire fittingwindow (from the near-apex region to the shank) by asingle hyperbolic curve (Supplementary Fig. S8(b)). Inthe hyperbolic parametrization, R and α are not indepen-dent: for a fixed R, decreasing α substantially increasesthe axial length scale over which the profile approachesits asymptotic shank taper. Consequently, for tips whoseshank narrows more rapidly (i.e., becomes nearly conicalclose to the apex), the full SEM profile cannot be cap-tured within the fitting window by a single hyperboliccurve, even when α is reduced to the smallest practicalvalue.We expect this residual mismatch to have only a lim-ited impact on the electrostatic response, because thenear-apex region (which dominates the electrostatic in-teraction for sharp metallic tips) is still reproduced rea-sonably well, whereas the contribution from the moredistant shank region is comparatively weak. In addi-tion, over the tip-radius range considered here (2–20 nm),the frequency shift signals changed only marginally whenα was reduced below 10◦, as shown in SupplementaryFig. S9. Based on these findings, we fixed the open-ing angle at α = 10◦ in all simulations used for fittingthroughout this work.Another reason for fixing α at 10◦ is the computationalcost associated with smaller opening angles. For smalleropening angles, the hyperbolic coordinate grid becomesincreasingly elongated near the tip surface, and a sub-stantially finer grid is required to maintain good consis-tency between the tip-side and sample-side force evalu-ations. This requirement becomes particularly demand-ing in the 3D implementation, where the computationalcost increases rapidly with grid density. Therefore, usingα = 10◦ provides a practical compromise: it is consistentwith the SEM-based estimate and the weak α dependenceof the frequency-shift response, while avoiding the exces-sive computational cost associated with smaller openingangles in the 3D calculations.2. Fitting to CPD profileTo compare the 3D simulation with experiment, weconsidered the sharp surface-potential-step model intro-duced in Section IV E. As an experimental test sample,we employed a Cu(111) surface partially covered with anultrathin NaCl overlayer [35], for which the CPD changesabruptly at the NaCl/Cu(111) interface [36]. The mea-surements were performed across the interface betweena two-monolayer-thick NaCl film and the Cu(111) sub-strate [47]. The blue curve in Fig. 5 shows an experimen-tal CPD profile, in which the CPD changes by approxi-mately 0.85 V across the interface. The sensor resonancefrequency and oscillation amplitude were 30.0 kHz and510 pm, respectively. An STM topographic image of themeasured region is provided in Supplementary Fig. S10.In the simulation, the potential difference between thetwo regions was fixed to the experimentally determinedvalue. For simplicity, the topographic height differencebetween the Cu and NaCl regions was neglected, and thesample surface was treated as geometrically flat. The in-fluence of this approximation is discussed at the end ofthis section. The solid circles in Fig. 5 show the simulatedprofile fitted to the experimental data. In the fitting, theshank opening angle was fixed at α = 10◦, while R ands were varied. The dashed vertical line marks the inter-face position in the simulation, which coincides with theonset of the NaCl-island topographic rise (green curve).Reasonable agreement was obtained for R = 3.0 nm ands = 0.9 nm. A small deviation is visible near the on-set of the Cu-to-NaCl transition, which may reflect theneglect of the topographic height difference between theCu and NaCl regions in the present flat-surface model,a slight deviation of the actual tip shape from the idealhyperbolic shape assumed in the model, or both.The width of the Cu-to-NaCl transition region de-pended strongly on R, becoming narrower for sharpertips (Supplementary Fig. S11). In contrast, varying swithin the physically reasonable range for STM-moderegulation (0.7–1.0 nm) produced no appreciable changein the CPD-profile shape (Supplementary Fig. S12). Sim-ilarly, varying the shank opening angle over the range7FIG. 5. Experimental CPD profile across a NaCl/Cu(111) in-terface (blue curve), together with the fitted simulated profile(red circles) and the corresponding topographic line profile(green curve). The dashed vertical line indicates the interfaceposition used in the simulation. Reasonable agreement wasobtained for R = 3.0 nm and s = 0.9 nm.α = 5◦–20◦ produced only a minor change in the CPD-profile width (Supplementary Fig. S13). These resultsindicate that, within the small-separation regime con-sidered here, R is the primary parameter governing theinterfacial broadening.To estimate the robustness of the extracted tip radius,we performed a residual-based sensitivity analysis basedon the root-mean-square error (RMSE) between the ex-perimental and simulated CPD profiles. The RMSEshowed a shallow minimum around R = 2.5–3.0 nm,as shown in Supplementary Fig. S14. When parametersets with RMSE values within 20% of the minimum wereregarded as acceptable, the CPD profile was consistentwith R ≈ 2.0–4.0 nm.These results demonstrate that the transition widthof a CPD profile across a potential-step interface can beused to constrain the electrically effective tip radius insitu. The representative value of R = 3.0 nm extractedfrom this CPD-profile fitting was then used as a fixed in-put parameter in an independent validation of the modelby fitting the ∆f–U response measured on the Cu region.We also assessed the possible influence of the approx-imately 320 pm apparent topographic height differencebetween the Cu(111) terrace and the 2 ML NaCl is-land, which was neglected in the flat-surface model. Be-cause the present 3D implementation assumes a geomet-rically flat sample surface, explicitly treating the non-flatboundary geometry would require a substantial extensionof the code. As a practical sensitivity test, we thereforeperformed an additional flat-surface simulation in whichthe tip–sample separation was increased by 320 pm (Sup-plementary Fig. S15). This calculation does not explic-itly model the topographic step itself, but evaluates theseparation-related effect associated with the measuredheight difference. The resulting CPD profile differed onlyslightly from the original calculation, indicating that thiseffect is unlikely to substantially affect the CPD-profilebroadening or the extracted effective tip-radius range.3. Fitting to ∆f–U curveFigure 6 shows the experimental ∆f–U curve ana-lyzed here, together with a representative METALTIPfit. The spectrum was extracted from the same datasetused for the CPD-profile analysis, specifically from a Curegion where the CPD profile was flat and the influence ofthe NaCl/Cu(111) interface was negligible. Accordingly,this analysis was carried out using the 2D implementa-tion of METALTIP. In the simulation, the tip radiuswas fixed at R = 3.0 nm, as determined from the CPD-profile fitting, and only the tip–sample separation wasvaried. Good agreement was obtained for a representa-tive value of s = 0.91 nm, which lies within the physi-cally reasonable range expected for STM-mode operation(0.7–1.0 nm) [48].To quantify the sensitivity of the ∆f–U fitting to thetip–sample separation s, we evaluated the RMSE be-tween the experimental and simulated curves as a func-tion of s (Supplementary Fig. S16). The RMSE showeda clear minimum around s = 0.91–0.92 nm, with nearlyidentical values at these two separations. When param-eter sets with RMSE values within 20% of the minimumwere regarded as acceptable, the experimental ∆f–Ucurve was consistent with s ≈ 0.90–0.93 nm. This rangeshould be regarded as the sensitivity of the fit underthe assumed model parameters, rather than as an ab-solute uncertainty bound on s. The absolute uncertaintyrange is likely larger than this estimate, because it is af-fected by systematic factors that are not fully capturedin the present sensitivity analysis. These include the un-certainty in the electrically effective tip radius, possibledeviations of the actual tip shape from the ideal hyper-bolic model, and calibration uncertainties in the sensorparameters such as k and A.Importantly, the representative fitted separation lieswithin the physically reasonable range expected for STM-mode regulation, and the experimental magnitude of the∆f–U curve was reproduced without any arbitrary pro-portionality factor or unphysical mechanical parameters.The spring constant k was set to the nominal design valueof the sensor, while the resonance frequency f0 and oscil-lation amplitude A were taken directly from the exper-imental settings. The conversion from force to ∆f wasperformed consistently using the established FM-AFMexpression. This scale-free agreement supports not onlythe validity but also the high quantitative accuracy ofthe overall framework, including the electric-field evalua-tion, electrostatic-force calculation, and force-to-∆f con-version.We performed the same combined analysis for sev-eral different tips. In many cases, the CPD profileswere well reproduced with effective tip radii in the range8���  ���� ���  ��� �� ��� �#��!����������������������������������"���%� ���!���&�
$�������!���"��!���FIG. 6. Experimental ∆f–U curve measured on a Cu(111)terrace (blue line) and a representative METALTIP fit (redcircles). In the simulation, α = 10◦ and R = 3.0 nm (fromthe CPD-profile fitting) were fixed, and the agreement shownwas obtained at s = 0.91 nm.of R = 2–14 nm, and the corresponding ∆f–U fittingyielded tip–sample separations of s = 0.8–0.9 nm. Insome cases, however, keeping s within a physically rea-sonable range (s = 0.7–1.0 nm) required the tip radiusused in the ∆f–U fitting to differ by a few nanometersfrom the value obtained from the CPD-profile fitting.This tendency may indicate that the actual tip geometrydeviates from the ideal hyperbolic geometry assumed inthe model. We also found tips for which the CPD pro-file showed a sharp transition at the interface but didnot become fully flat even in regions sufficiently far fromthe boundary (Supplementary Fig. S17). Such behaviormay reflect an actual tip geometry that remains relativelysharp near the apex but becomes broader away from theapex than assumed in the ideal hyperbolic model.Finally, we discuss a simple strategy for tip-parameterextraction using METALTIP. Under STM-mode dis-tance regulation, the variation in s is expected to be lim-ited as long as the same setpoint is used. Moreover, fortips fabricated by the same procedure, the variation inthe shank opening angle is also expected to be relativelysmall. These considerations suggest that, by acquiring a∆f–U curve on a clean metallic surface at a fixed STMsetpoint and fitting it with the 2D METALTIP modelwhile treating R as the primary fitting parameter, theelectrically effective tip radius can be evaluated in situunder experimental conditions.V. SUMMARYWe developed a numerical software package forcomputing electrostatic interactions between a hy-perbolic metallic tip and a conducting surface un-der AFM/KPFM-relevant conditions. This code en-ables direct simulation of the bias-dependent force andfrequency-shift responses, F (U) and ∆f(U). In addi-tion, the full 3D implementation treats spatially varyingsurface potentials and reproduces CPD profiles acrossinhomogeneous interfaces. We validated the numericalframework by benchmarking against established theoret-ical models, confirming the expected scaling behavior andcorrect near-field limiting behavior. Using representativeexperimental datasets, we further demonstrated quanti-tative fitting without arbitrary scale factors, reproducingboth CPD profiles across a potential-step interface and∆f–U spectra on a homogeneous metal surface. Combin-ing these complementary observables provides strongerconstraints on key experimental parameters—includingthe tip–sample separation and electrically effective tipradius—than fitting either dataset alone, and therebyenables in situ parameter identification with high quan-titative accuracy. These results establish METALTIPas a practical platform for extracting key experimen-tal parameters, enabling quantitative interpretation ofKPFM data, analyzing contrast formation and measure-ment characteristics, and testing emerging KPFM mea-surement schemes under realistic tip–sample geometries.SUPPLEMENTARY MATERIALSee the supplementary material for additional nu-merical details, sample preparation, supplementary fig-ures, and supporting analyses, including the numeri-cal implementation of METALTIP, simulated surface-potential profiles, parameter-dependence plots, represen-tative scanning electron microscopy images of the tips,and additional CPD profiles.DATA AVAILABILITY STATEMENTThe METALTIP source code is publicly avail-able on GitHub at https://github.com/nobuyuki-ishida-nims/metaltip. The data that support the findings ofthis study are available from the corresponding authorupon reasonable request.ACKNOWLEDGMENTSThis work was partially supported by JSPS KAK-ENHI Grant Numbers JP17K06366, JP21H01818, andJP24K01367.9[1] M. Nonnenmacher, M. P. O’Boyle, and H. K. Wickra-masinghe, Appl. Phys. Lett. 58, 2921 (1991).[2] W. Melitz, J. Shen, A. C. Kummel, and S. Lee, Kelvinprobe force microscopy and its application, Surf. Sci.Rep. 66, 1 (2011).[3] S. Sadewasser and T. Glatzel, eds., Kelvin Probe ForceMicroscopy: Measuring and Compensating ElectrostaticForces (Springer, 2012).[4] R. Shikler, T. Meoded, N. Fried, and Y. Rosenwaks, Po-tential imaging of operating light-emitting devices usingkelvin force microscopy, Applied Physics Letters 74, 2972(1999).[5] H. Masuda, N. Ishida, Y. Ogata, D. Ito, and D. Fujita,Nanoscale 9, 893 (2017).[6] M. Cai, N. Ishida, X. Li, X. Yang, T. Noda, Y. Wu,F. Xie, H. Naito, D. Fujita, and L. Han, Control ofelectrical potential distribution for high-performance per-ovskite solar cells, Joule 2, 296 (2018).[7] H. Masuda, K. Matsushita, D. Ito, D. Fujita, andN. Ishida, Commun. Chem. 2, 140 (2019).[8] T. Nakamura, N. Ishida, K. Sagisaka, and Y. Koide, Sur-face potential imaging and characterizations of a gan p-njunction with kelvin probe force microscopy, AIP Ad-vances 10, 085010 (2020).[9] N. Ishida, Utilizing the surface potential of a solid elec-trolyte region as the potential reference in kelvin probeforce microscopy, Beilstein Journal of Nanotechnology13, 1558 (2022).[10] S. Belaidi, P. Girard, and G. Leveque, Electrostatic forcesacting on the tip in atomic force microscopy: Modeliza-tion and comparison with analytic expressions, Journalof Applied Physics 81, 1023 (1997).[11] H. O. Jacobs, P. Leuchtmann, O. J. Homan, and A. Stem-mer, Resolution and contrast in kelvin probe force mi-croscopy, Journal of Applied Physics 84, 1168 (1998).[12] S. Hudlet, M. Saint Jean, C. Guthmann, and J. Berger,Evaluation of the capacitive force between an atomicforce microscopy tip and a metallic surface, The Euro-pean Physical Journal B 2, 5 (1998).[13] J. Colchero, A. Gil, and A. M. Baró, Resolution enhance-ment and improved data interpretation in electrostaticforce microscopy, Physical Review B 64, 245403 (2001).[14] A. Sadeghi, A. Baratoff, S. A. Ghasemi, S. Goedecker,T. Glatzel, S. Kawai, and E. Meyer, Multiscale approachfor simulations of Kelvin probe force microscopy withatomic resolution, Phys. Rev. B 86, 075407 (2012).[15] G. Cohen, E. Halpern, S. U. Nanayakkara, J. M. Luther,C. Held, R. Bennewitz, A. Boag, and Y. Rosenwaks, Re-construction of surface potential from Kelvin probe forcemicroscopy images, Nanotechnology 24, 295702 (2013).[16] J. L. Neff and P. Rahe, Insights into Kelvin probe forcemicroscopy data of insulator-supported molecules, Phys-ical Review B 91, 085424 (2015).[17] J. F. Gonzalez, A. M. Somoza, and E. Palacios-Lidón,Charge distribution from SKPM images, Physical Chem-istry Chemical Physics 19, 27299 (2017).[18] M. F. Orihuela, A. M. Somoza, J. Colchero, M. Ortuño,and E. Palacios-Lidón, Localized charge imaging withscanning Kelvin probe microscopy, Nanotechnology 28,025703 (2017).[19] J. Xu, Y. Wu, W. Li, and J. Xu, Surface potential model-ing and reconstruction in Kelvin probe force microscopy,Nanotechnology 28, 365705 (2017).[20] F. Pertl, J. C. Sobarzo, L. Shafeek, T. Cramer, andS. Waitukaitis, Quantifying nanoscale charge density fea-tures of contact-charged surfaces with an FEM/KPFM-hybrid approach, Physical Review Materials 6, 125605(2022).[21] U. Zerweck, C. Loppacher, T. Otto, S. Grafström, andL. M. Eng, Accuracy and resolution limits of kelvin probeforce microscopy, Phys. Rev. B 71, 125424 (2005).[22] Y. Shen, M. Lee, W. Lee, D. M. Barnett, P. M. Pinsky,and F. B. Prinz, A resolution study for electrostatic forcemicroscopy on bimetallic samples using the boundary el-ement method, Nanotechnology 19, 035710 (2007).[23] T. Wagner, H. Beyer, P. Reissner, P. Mensch, H. Riel,B. Gotsmann, and A. Stemmer, Kelvin probe force mi-croscopy for local characterisation of active nanoelec-tronic devices, Beilstein Journal of Nanotechnology 6,2193 (2015).[24] G. Elias, T. Glatzel, E. Meyer, A. Schwarzman, A. Boag,and Y. Rosenwaks, The role of the cantilever in Kelvinprobe force microscopy measurements, Beilstein Journalof Nanotechnology 2, 252 (2011).[25] W. A. Behn, Z. J. Krebs, K. J. Smith, K. Watanabe,T. Taniguchi, and V. W. Brar, Measuring and tuningthe potential landscape of electrostatically defined quan-tum dots in graphene, Nano Letters 21, 5013 (2021),https://doi.org/10.1021/acs.nanolett.1c00791.[26] N. Ishida and T. Mano, Quantitative characterization ofbuilt-in potential profile across gaas p–n junctions us-ing kelvin probe force microscopy with qplus sensor afm,Nanotechnology 35, 065708 (2023).[27] S. Belaidi, F. Lebon, P. Girard, G. Leveque, andS. Pagano, Finite element simulations of the resolutionin electrostatic force microscopy, Applied Physics A 66,S239 (1998).[28] R. M. Feenstra, Electrostatic potential for a hyperbolicprobe tip near a semiconductor, Journal of Vacuum Sci-ence & Technology B 21, 2080 (2003).[29] R. M. Feenstra, S. Gaan, G. Meyer, and K.-H.Rieder, Low-temperature tunneling spectroscopy ofge(111)c(2×8) surfaces, Physical Review B 71, 125316(2005).[30] R. M. Feenstra, Y. Dong, M. P. Semtsiv, and W. T.Masselink, Influence of tip-induced band bending on tun-nelling spectra of semiconductor surfaces, Nanotechnol-ogy 18, 044015 (2007).[31] Y. Dong, R. M. Feenstra, M. P. Semtsiv, and W. T.Masselink, Band offsets of ingap/gaas heterojunctionsby scanning tunneling spectroscopy, Journal of AppliedPhysics 103, 073704 (2008).[32] N. Ishida, K. Sueoka, and R. M. Feenstra, Influence ofsurface states on tunneling spectra of n-type gaas(110)surfaces, Physical Review B 80, 075320 (2009).[33] S. Gaan, G. He, R. M. Feenstra, J. Walker, andE. Towe, Size, shape, composition, and electronic prop-erties of inas/gaas quantum dots by scanning tunnelingmicroscopy and spectroscopy, Journal of Applied Physics108, 114315 (2010).10[34] D. J. Griffiths, Introduction to Electrodynamics (Cam-bridge University Press, 2017).[35] R. Bennewitz, V. Barwich, M. Bammerlin, C. Loppacher,M. Guggisberg, A. Baratoff, E. Meyer, and H.-J. Gün-therodt, Ultrathin films of nacl on cu(111): a leed anddynamic force microscopy study, Surface Science 438,289 (1999).[36] T. Glatzel, L. Zimmerli, S. Koch, B. Such, S. Kawai, andE. Meyer, Determination of effective tip geometries inkelvin probe force microscopy on thin insulating films onmetals, Nanotechnology 20, 264016 (2009).[37] F. J. Giessibl, Atomic resolution on si(111)-(7×7) by non-contact atomic force microscopy with a force sensor basedon a quartz tuning fork, Applied Physics Letters 76, 1470(2000).[38] F. J. Giessibl, The qplus sensor, a powerful core for theatomic force microscope, Review of Scientific Instruments90, 011101 (2019).[39] T. R. Albrecht, P. Grütter, D. Horne, and D. Rugar, Fre-quency modulation detection using high-q cantilevers forenhanced force microscope sensitivity, Journal of AppliedPhysics 69, 668 (1991).[40] G. H. Simon, M. Heyde, and H.-P. Rust, Recipes for can-tilever parameter determination in dynamic force spec-troscopy: spring constant and amplitude, Nanotechnol-ogy 18, 255503 (2007).[41] T. Vančura, S. Kičin, T. Ihn, K. Ensslin, M. Bichler,and W. Wegscheider, Kelvin probe spectroscopy of a two-dimensional electron gas below 300 mk, Appl. Phys. Lett.83, 2602 (2003).[42] G. Münnich, A. Donarini, M. Wenderoth, and J. Repp,Fixing the energy scale in scanning tunneling microscopyon semiconductor surfaces, Phys. Rev. Lett. 111, 216802(2013).[43] F. Albrecht, M. Fleischmann, M. Scheer, L. Gross, andJ. Repp, Local tunneling decay length and kelvin probeforce spectroscopy, Phys. Rev. B 92, 235443 (2015).[44] W. R. Smythe, Static and Dynamic Electricity (McGraw-Hill, New York, 1968).[45] F. J. Giessibl, A direct method to calculate tip–sampleforces from frequency shifts in frequency-modulationatomic force microscopy, Applied Physics Letters 78, 123(2001).[46] F. J. Giessibl, Advances in atomic force microscopy, Re-views of Modern Physics 75, 949 (2003).[47] T. Kakudate, M. Nakaya, and T. Nakayama, Local mod-ification of nacl thin films on cu(111) under different biasvoltages, Thin Solid Films 520, 2004 (2012).[48] F. Albrecht, J. Repp, M. Fleischmann, M. Scheer, M. On-dráček, and P. Jelínek, Probing charges on the atomicscale by means of atomic force microscopy, Phys. Rev.Lett. 115, 076101 (2015).