# Fileset

[s-info5.pdf](https://mdr.nims.go.jp/filesets/8209abb9-eda4-49dd-a298-6c4dda48c5dc/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

Supplementary information:An Open-Source Framework for Quantitative ElectrostaticSimulations in Kelvin Probe Force MicroscopyNobuyuki Ishida11National Institute for Materials Science,1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan1I. SUPPLEMENTARY NOTE I: NUMERICAL IMPLEMENTATION OF MET-ALTIP IN GENERALIZED PROLATE SPHEROIDAL COORDINATESThis Supplementary Note provides a complete description of the numerical implementa-tion of METALTIP, including the coordinate definitions, parameter relations, discretizationscheme, boundary conditions, and solver procedures. The overall computational strategylargely follows the SEMITIP framework developed by Feenstra [1–6].A. Generalized prolate spheroidal coordinatesIn METALTIP, following the SEMITIP framework, the vacuum region between the tipand the sample is described using a generalized prolate spheroidal coordinate system (ξ, η).A key advantage of this coordinate choice is that the tip boundary can be aligned with acoordinate surface, η = ηT , which allows the boundary condition on the tip to be imposedstraightforwardly while retaining a compact grid representation of a sharp apex and thevacuum gap.The generalized coordinates (ξ, η) are defined in the (r, z) cross-section of the cylindricalcoordinate system, where r is the radial coordinate and z is the axial coordinate, by theimplicit relationsr2ξ2 − 1+(z − ac η)2ξ2= a2, (1)r21− η2+(z − ac η)2η2= a2, (2)where a sets the characteristic length scale of the mapping (equal to the focal length in thestandard prolate spheroidal system) and c is a dimensionless parameter introduced in thegeneralized formulation. Introducing c adds geometric flexibility: it allows the tip surface tocoincide with a coordinate surface while enabling the tip–sample separation (vacuum gap)to be tuned independently, as detailed in Sec. I B. Surfaces of constant ξ and constant ηcorrespond to spheroidal (ellipsoidal) and hyperboloidal quadric surfaces, respectively (seeFig. S1).Solving Eqs. (1) and (2) for (r, z) yields an explicit mapping to the physical space in2FIG. S1: Example of the grid in the generalized prolate spheroidal coordinates (ξ, η) used inSEMITIP/METALTIP. The ξ = const. and η = const. lines form spheroidal and hyperboloidalfamilies, respectively. For c 6= 0 the coordinates are generally non-orthogonal.Alt text: Coordinate grid used in SEMITIP/METALTIP. Lines of constant ξ and η formspheroidal and hyperboloidal coordinate families, and the generalized coordinate system isnon-orthogonal when c 6= 0.cylindrical coordinates:r(ξ, η) = a√(ξ2 − 1)(1− η2), (3)z(ξ, η) = a(ξ + c)η. (4)With this definition, the sample surface coincides with η = 0 (i.e., z = 0). The tip boundaryis represented by choosing a coordinate surface, η = ηT , to coincide with the tip surface.The definition of ηT in terms of the tip geometric parameters is summarized in Sec. I B.Figure S1 shows a representative grid of the generalized prolate spheroidal coordinatesused in METALTIP. Unlike the standard (orthogonal) prolate spheroidal coordinates, thepresent generalized system with c 6= 0 is generally non-orthogonal; the coordinate lines ofconstant ξ and constant η are not mutually perpendicular in the (r, z) cross-section. Thisnon-orthogonality introduces mixed ξ–η terms in the Laplacian (see Sec. I E), while retaininga compact grid in the near-apex region where the electrostatic field is strongly localized.3B. Determination of the tip boundary parameter ηTThe tip is modeled as a rotationally symmetric hyperboloid, specified by the tip radiusR, the full opening angle α, and the tip–sample separation s. Note that ηT is determinedsolely by the opening angle α, whereas R sets the scale parameter a and s sets the shiftparameter c.We define the shank-slope parameter as the asymptotic slope in the (r, z) cross-section,b ≡ dzdr∣∣∣∣asym= cot(α2). (5)The tip surface is aligned with the coordinate surface η = ηT , whereηT =1√1 + b−2= cos(α2). (6)The scale parameter a is determined by requiring the local apex curvature to match R,a =Rb2ηT. (7)Defining s′ ≡ aηT , the vertical shift is chosen such that the apex position corresponds to theprescribed separation s:z0 = s− s′, c =z0aηT=saηT− 1. (8)These relations ensure that the coordinate surface η = ηT coincides with the hyperbolictip boundary and that the apex point (ξ, η) = (1, ηT ) is mapped to (r, z) = (0, s) with thedesired local curvature R.C. Grid construction and indexingThe η coordinate is discretized uniformly between the sample plane (η = 0) and the tipsurface (η = ηT ) with∆η =ηTNv. (9)The interior vacuum nodes are indexed as j = 0, 1 . . . , Nv − 1 and are placed atηj = (j + 1)∆η, (j = 0, 1, . . . , Nv − 1). (10)The sample surface at η = 0 is not included as a grid point; the first vacuum node is locatedat η = ∆η.4To maintain high resolution near the axis while keeping a sufficiently large lateral domain,we introduce a variable (non-uniform) grid in the radial direction using a tangent-typemapping,ri =2Nr∆rπtan[π2Nr(i+12)], (i = 0, 1, . . . , Nr − 1), (11)where ∆r controls the near-axis spacing. Note that the central axis at r = 0 (correspondingto ξ = 1) is not included as a grid point. The corresponding ξi values are determined fromthe sample-plane relation at η = 0,ξi =√1 +(ria)2. (12)Once {ξi} and {ηj} are defined, the physical-space grid points follow from Eqs. (3) and (4).In particular, the radial and vertical coordinates at a general grid point (i, j) are given byrij = a√(ξ2i − 1)(1− η2j ) = ri√1− η2j , (13)zij = a(ξi + c)ηj = aηj(√1 +(ria)2+ c). (14)For the three-dimensional (3D) solver, an azimuthal angle θ is introduced and discretizeduniformly,θk = k∆θ, ∆θ =2πNp, k = 0, 1, . . . , Np − 1, (15)with periodic boundary conditions identifying θ−1 ≡ θNp−1 and θNp ≡ θ0.D. Boundary conditionsThe boundary conditions are specified as follows:1. Sample surface (η = 0): the surface potential is prescribed as a Dirichlet boundarycondition,φ(ξ, η = 0, θ) = Vs(ξ, θ). (16)Here Vs(ξ, θ) denotes the surface potential defined on the η = 0 plane, i.e., Vs(ξ, θ) ≡Vs(r(ξ, 0), θ). For the 2D calculations, we assume a grounded sample (Vs = 0).2. Tip surface (η = ηT ): the tip is treated as an ideal conductor at the applied bias U ,φ(ξ, η = ηT , θ) = U. (17)53. Symmetry axis (r = 0, ξ = 1): axial symmetry yields the Neumann condition,∂φ∂r∣∣∣∣r=0= 0. (18)4. Outer boundary (ξ = ξmax): the outer boundary is placed sufficiently far from theapex, and a homogeneous Neumann-type condition is imposed in the ξ direction.∂φ∂ξ∣∣∣∣ξ=ξmax= 0. (19)As mentioned in Section I C, the plane η = 0 is not explicitly included as computationalnodes; the Dirichlet condition on the sample is therefore enforced using a ghost layer adjacentto the first vacuum layer. In particular, when updating the first vacuum layer (j = 0,η = ∆η), the surface Dirichlet condition is applied by assigning the prescribed value Vs as aghost value at j = −1.The symmetry axis ξ = 1 (r = 0) is a coordinate singularity and is not explicitly includedin the grid. Regularity at the axis requires that the potential remain finite. Accordingly, theaxis condition is enforced in a finite-difference sense by introducing an on-axis ghost valueobtained from a short extrapolation stencil. In practice, the on-axis potential is evaluatedfrom the nearest off-axis points (e.g., i = 0 and i = 1) asφax,j,k ≈98φ0,j,k −18φ1,j,k, (20)which is consistent with ∂φ/∂r|r=0 = 0.E. Discretization of the Laplace equation in the generalized (ξ, η) coordinatesIn the vacuum region, the electrostatic potential φ satisfies Laplace equation, ∇2φ = 0.In the generalized prolate spheroidal coordinates (ξ, η, θ), this equation can be written asthe following second-order partial differential equation:f1(ξ, η)∂2φ∂ξ2+ f2(ξ, η)∂2φ∂η2+ f3(ξ, η)∂2φ∂θ2+ f4(ξ, η)∂2φ∂ξ ∂η+ f5(ξ, η)∂φ∂ξ+ f6(ξ, η)∂φ∂η= 0, (21)6where the coefficient functions fn(ξ, η) are given byf1(ξ, η) =(ξ2 − 1) [(ξ + c)2 − η2(2cξ + c2 + 1)]ξ(ξ + c)− η2(cξ + 1), (22)f2(ξ, η) =(1− η2)(ξ2 − η2)ξ(ξ + c)− η2(cξ + 1), (23)f3(ξ, η) =ξ(ξ + c)− η2(cξ + 1)(ξ2 − 1)(1− η2), (24)f4(ξ, η) =−2cη(ξ2 − 1)(1− η2)ξ(ξ + c)− η2(cξ + 1), (25)f5(ξ, η) =g5(ξ, η)[ξ(ξ + c)− η2(cξ + 1)]2, (26)f6(ξ, η) =g6(ξ, η)[ξ(ξ + c)− η2(cξ + 1)]2, (27)withg5(ξ, η) = c3 + 3c2ξ + c(2 + c2)ξ2 + 3c2ξ3 + 4cξ4 + 2ξ5+ η4[c3 + (2 + 3c2)ξ + c(6 + c2)ξ2 + 3c2ξ3]− 2η2[c(c2 − 1) + 3c2ξ + c(6 + c2)ξ2 + (2 + 3c2)ξ3 + cξ4], (28)g6(ξ, η) = −η{c2 + 4cξ + c2ξ2 + 2ξ4 + η4[2 + c2 + 4cξ + c2ξ2]− 2η2[c2 + 4cξ + (2 + c2)ξ2]}. (29)Under axial symmetry, φ is independent of the azimuthal angle θ, so that ∂φ/∂θ = 0and the θ-derivative term vanishes, i.e., f3(ξ, η) ∂2φ/∂θ2 = 0. The vacuum problem is thenreduced to a two-variable equation in (ξ, η). In the generalized system (c 6= 0), the mixedξ–η term remains.Equation (21) is discretized on the (ξ, η, θ) grid using second-order finite differences. Theη and θ directions use uniform spacings ∆η and ∆θ, whereas the ξ direction follows thenon-uniform grid described in Sec. I C. We define the local ξ spacings as∆ξ−i ≡ ξi − ξi−1, ∆ξ+i ≡ ξi+1 − ξi, ∆ξi ≡ ∆ξ−i +∆ξ+i . (30)At the symmetry axis, ξ = 1 is treated as a ghost boundary; in the implementation, ∆ξ−0 istaken as the distance from ξ = 1 to ξ0.The discrete equation at an interior vacuum node (i, j, k) can be written as a weighted7sum of neighboring values (a stencil):A+ijφi+1,j,k + A−ijφi−1,j,k +B+ijφi,j+1,k +B−ijφi,j−1,k+ Cij(φi+1,j+1,k − φi+1,j−1,k − φi−1,j+1,k + φi−1,j−1,k)+ Pij (φi,j,k+1 + φi,j,k−1)−Dijφi,j,k = 0. (31)where the coefficients areA+ij = f1(ξi, ηj)2∆ξ+i ∆ξi+ f5(ξi, ηj)1∆ξi, (32)A−ij = f1(ξi, ηj)2∆ξ−i ∆ξi− f5(ξi, ηj)1∆ξi, (33)B+ij = f2(ξi, ηj)1∆η2+ f6(ξi, ηj)12∆η, (34)B−ij = f2(ξi, ηj)1∆η2− f6(ξi, ηj)12∆η, (35)Cij = f4(ξi, ηj)12∆η∆ξi, (36)Pij = f3(ξi, ηj)1∆θ2, (37)Dij =(A+ij + A−ij)+(B+ij +B−ij)+ 2Pij. (38)For numerical efficiency, these coefficients are precomputed and stored in the code as Ap=A+, Am=A−, Bp=B+, Bm=B−, C=C, P=P , and inv_den=1/D.Equation (31) is solved iteratively using Gauss–Seidel iteration with successive over-relaxation (SOR). At each iteration (sweep), the interior vacuum nodes are updated in place(i.e., array entries are overwritten during the sweep). We first compute the Gauss–Seidelestimate by rearranging Eq. (31) asφGSi,j,k = inv_denij[Apijφi+1,j,k +Amijφi−1,j,k + Bpijφi,j+1,k + Bmijφi,j−1,k+ Cij(φi+1,j+1,k − φi+1,j−1,k − φi−1,j+1,k + φi−1,j−1,k)+ Pij (φi,j,k+1 + φi,j,k−1)], (39)and then apply SOR,φi,j,k ← (1− ω)φi,j,k + ω φGSi,j,k, (40)where ω is the over-relaxation factor. Dirichlet boundary nodes are excluded from theupdates. The iteration is continued until the maximum change in φ between successive8sweeps falls below a prescribed tolerance,maxi,j,k∣∣∣φ(n+1)i,j,k − φ(n)i,j,k∣∣∣ < εφ. (41)To improve numerical accuracy while keeping the computational cost manageable, weemploy a staged grid-refinement procedure. Starting from an initial coarse grid, the Laplacesolver is executed through a sequence of refinement steps (typically five or six steps in total),in which the grid density is doubled at each step. The converged solution at a given step isinterpolated onto the refined grid and used as the initial guess for the subsequent step.F. Electric-field evaluation on the conductor surfacesWe next describe how the electric field is evaluated from the numerically computed poten-tial. After obtaining φ(ξ, η, θ), we compute the electric field in physical space as E = −∇φ.To evaluate the physical-space derivatives, we first compute ∂φ/∂ξ and ∂φ/∂η on the coor-dinate grid and then transform them to physical-space derivatives using the Jacobian of thecoordinate mapping.For the tip surface (η = ηT ), the normal field component is evaluated using a higher-accuracy one-sided boundary derivative to reduce discretization error. Specifically, insteadof using a single-step finite difference at the boundary, we sample the potential at several gridpoints (seven points) along the η coordinate line in the vacuum region immediately adjacentto the tip surface (toward decreasing η at fixed ξ), and evaluate ∂φ/∂η|η=ηTusing a one-sidedfinite-difference stencil (equivalently, polynomial interpolation followed by differentiation)[7]. Because the tip surface is an equipotential (Dirichlet boundary), the tangential derivativealong the surface vanishes; accordingly, ∂φ/∂ξ|η=ηT= 0.The derivative information is then converted to the physical-space gradient components(∂φ/∂r, ∂φ/∂z) using the Jacobian of the mapping r(ξ, η) and z(ξ, η). DefiningJ ≡rξ rηzξ zη , J−1 =ξr ξzηr ηz , (42)the chain rule gives∂φ/∂r∂φ/∂z =(J−1)T ∂φ/∂ξ∂φ/∂η =ξr ηrξz ηz∂φ/∂ξ∂φ/∂η . (43)9Here, rξ = ∂r/∂ξ, etc., follow from Eqs. (3) and (4). The resulting ∇φ is perpendicular tothe equipotential tip surface, and its magnitude gives the normal field component on thetip, |En| = |∇φ|.For the 2D calculations, the sample boundary at η = 0 is laterally uniform, so thatthe electric field at the sample plane has only the normal component. The normal fieldis obtained from a one-sided derivative along the η direction. Specifically, the potential issampled along the η line starting from the Dirichlet boundary value at the surface and severaladjacent vacuum nodes, and ∂φ/∂η|η=0 is evaluated using a one-sided finite-difference stencil[7]. Using the local relation z(ξ, η) = a(ξ + c)η, this derivative is converted to the verticalderivative at the sample plane as∂φ∂z∣∣∣∣z=0=1a(ξ + c)∂φ∂η∣∣∣∣η=0, (44)which gives Ez = − ∂φ/∂z|z=0.For the 3D calculations, in contrast, an arbitrary laterally varying surface potentialVs(ξ, θ) is prescribed at η = 0. The normal field component is evaluated in the same manneras in the 2D case, namely from a one-sided derivative along η followed by the transformation∂φ∂z∣∣∣∣z=0=1a(ξ + c)∂φ∂η∣∣∣∣η=0, (45)so that Ez = − ∂φ/∂z|z=0. However, because Vs varies laterally, tangential field componentsgenerally remain finite on the sample plane. In cylindrical coordinates, these are evaluatedasEr = −∂Vs∂r, Eθ = −1r∂Vs∂θ. (46)In the implementation, Er is evaluated from finite differences on the nonuniform radial grid,whereas Eθ is obtained from a periodic finite difference in θ.G. Maxwell-stress force evaluationThe electrostatic force is computed from the Maxwell stress on the conductor surfaces.For the tip side, the electrostatic pressure in vacuum is given byp =12ε0E2n, (47)10where En is the electric-field component normal to the surface. The vertical force is thenobtained by integrating the pressure over the tip surface,Fz =∫∫p nz dS, (48)where nz is the z component of the unit normal vector and dS is the surface area element.For the planar sample surface in the 2D calculations, no tangential field component ispresent, and the vertical traction reduces to the usual electrostatic pressure,Tzz =ε02E2z . (49)The corresponding sample-side force is thereforeF samplez = 2π∫ rmax0Tzz(r) r dr, (50)which is evaluated numerically on the discrete grid.For the 3D calculations, because laterally varying surface potentials generally produceboth normal and tangential electric-field components on the sample plane, the vertical com-ponent of the Maxwell stress is evaluated asTzz =ε02(E2z − E2r − E2θ). (51)The total sample-side vertical force is then obtained by integrating Tzz over the samplesurface,Fz =∫∫Tzz r dr dθ. (52)Thus, unlike in the 2D case, the 3D force evaluation explicitly includes the contributions ofboth the normal and lateral electric-field components on the sample surface.H. Internal consistency check of tip- and sample-side forcesAs an internal consistency check (action–reaction), the mismatch between the forcesobtained from the tip and sample surfaces is quantified asδF ≡2∣∣|F tipz | − |F samplez |∣∣|F tipz |+ |F samplez |. (53)We find that δF decreases upon refining the η-direction mesh, reflecting the sensitivity of thetip-side evaluation to the strongly localized field near the apex (note that p ∝ E2n amplifies11discretization-induced errors). As described above, the calculations use staged refinement;here we only summarize the practical choice of the baseline (coarsest-level) η-grid numberNv,in. For tips with relatively large opening angles (α ≥ 15◦), Nv,in = 4 or 8 is typicallysufficient to achieve agreement between the tip- and sample-side force evaluations within afew percent, whereas for sharper tips (α < 15◦), Nv,in ≥ 16 is generally required to reach acomparable level of agreement.Importantly, when Nv,in is varied, the tip-side force F tipz changes substantially more thanthe sample-side force F samplez , while F samplez remains comparatively stable. This indicatesthat the residual mismatch is dominated by discretization sensitivity on the tip side, andthat increasing the η resolution mainly reduces this tip-side contribution to the mismatch.In principle, the electrostatic force acting between the tip and sample can be evaluated fromeither the tip-side or sample-side stress integration alone. Therefore, in practice, it is notnecessary to increase the grid resolution until F tipz fully matches F samplez . In the presentcalculations, to reduce computational cost in extensive parameter sweeps, we use a modestNv,in and evaluate the force from the sample-side stress integration, while performing addi-tional refinements only when needed to confirm the stability of F samplez and the consistencyof the result.I. Conversion from static STM setpoint distance to mean oscillation distanceIn ultrahigh vacuum atomic force microscopy (AFM)/Kelvin probe force microscopy(KPFM) measurements, the tip–sample distance is often controlled by scanning tunnel-ing microscopy (STM) feedback, i.e., by using tunneling current. However, in AFM/KPFMmeasurements, the tip is mechanically oscillated, and the tip–sample distance is thereforeeffectively regulated through the tunneling current averaged over the oscillation cycle. Con-sequently, for a given STM setpoint, the mean tip position during oscillatory operation islarger than the corresponding position in a non-oscillating (static) configuration.In METALTIP, for convenience and practical parameterization, we use as the inputparameter a nominal tip–sample separation s defined for the static (non-oscillating) STMcase. This choice is convenient because the corresponding distance range is commonly usedas an experimental reference (typically on the order of ∼0.7–1.0 nm under usual STMconditions), whereas directly specifying the mean distance during oscillatory operation would12require explicit amplitude-dependent correction. The mean tip position 〈s〉 used in theelectrostatic potential and force calculations is then determined so that the oscillation-averaged tunneling current matches the tunneling current at the static reference distances.Assuming an exponential distance dependence of the tunneling current,I ∝ exp(−2κs), (54)and writing the instantaneous separation during FM operation ass(t) = 〈s〉+ A cos(2πf0t), (55)the matching condition for the tunneling current is〈exp[−2κ (〈s〉+ A cosωt)]〉 = exp(−2κs), (56)where ω = 2πf0. The left-hand side can be separated into the contribution from the meantip position and the oscillatory contribution:exp(−2κ〈s〉) 〈exp[−2κA cosωt]〉 = exp(−2κs). (57)Using the integral representation of the modified Bessel function,I0(x) =12π∫ 2π0exp(x cos θ) dθ, (58)and noting that I0(−x) = I0(x), the oscillation average is given by〈exp[−2κA cosωt]〉 = I0(2κA). (59)The current-matching condition therefore becomesexp(−2κ〈s〉)I0(2κA) = exp(−2κs), (60)which yields〈s〉 = s+12κln[I0(2κA)] , (61)where I0 is the modified Bessel function of the first kind. In the simulations, 〈s〉 is used as themean tip position for evaluating electrostatic quantities and oscillation-averaged observablessuch as ∆f .13II. SUPPLEMENTARY NOTE II: SAMPLE PREPARATIONA Cu(111) single-crystal substrate was cleaned by repeated cycles of Ar+-ion sputteringand annealing. After cleaning, the sample was transferred to the cryostat of the scanningprobe microscopy (SPM) chamber and cooled to ∼78 K. It was then transferred to theroom-temperature preparation chamber, where NaCl was deposited immediately by thermalevaporation from a glass crucible heated to 580 ◦C for 3 min. After deposition, the samplewas left in the preparation chamber for 20 min and was then transferred back to the SPMchamber. This procedure resulted in the formation of rectangular-shaped NaCl islands.14III. SUPPLEMENTARY FIGURES� ��� ��� ��� ���  �� ���������������������������� ����������������������������������������� �������� ������
�������������
����������FIG. S2: Radial dependence of the electrostatic pressure on a logarithmic vertical scale, evaluatedon the tip surface (blue circles) and on the sample surface (red circles) for a sharp hyperbolic tipgeometry (R = 10 nm, α = 10◦) at an applied bias of 1.0 V. Near the apex region, the pressureis larger on the tip side, whereas the weak long-range tail is more pronounced on the sample side.Although this tail is small in magnitude, it contributes to the total force through the area-weightedsurface integral (∝ r dr). Therefore, a sufficiently large radial integration range is required to obtainquantitative agreement between the tip-side and sample-side force integrations (typically within∼5%).Alt text: Log-scale radial profiles of electrostatic pressure on the tip and sample surfaces for asharp hyperbolic tip. The tip-side pressure is larger near the apex, whereas the sample-sidepressure shows a larger long-range tail at larger radial distances, indicating different spatialweightings of the two surface integrations.15� ��� ��� ��� ���  �������������������������������μ��������� ������� ������� ������������������������������
������������������
��������������������FIG. S3: Cumulative vertical electrostatic force on the tip side as a function of the vertical distancefrom the apex. The solid line shows the force obtained by integrating the calculated tip-sideelectrostatic pressure within the computational domain, which extends to approximately 18 µmfrom the apex for a full opening angle of 10◦. The dashed line shows an extrapolation of the outertail of the pressure distribution to 500 µm, chosen to cover the typical length scale of the actualtip used in self-sensing force sensors. The cumulative force is already close to saturation withinthe modeled hyperbolic-tip region, indicating that the additional contribution from more distantmetallic parts of the probe is expected to be small.Alt text: Cumulative vertical electrostatic force plotted as a function of distance along the zdirection from the tip apex. The curve approaches saturation within the modeled hyperbolic-tipregion, supporting the assumption that electrostatic contributions from more distant probe partsare small.16���� ��� ��� ��� ��� �� ���������"#����������������������������������� ��������#!�"#�#��� !�""$!����������
������"�� ���� ��!�� ������"�� ���FIG. S4: Logarithmic-scale comparison of the sample-side electrostatic pressure profile calculatedwith METALTIP (solid circles) and the analytical sphere–plane model (dashed line). The twoprofiles agree well in the near-apex region (small r), where the local apex curvature dominates theelectrostatic response. At larger radial distances, however, the METALTIP result increasinglydeviates from the sphere–plane prediction, indicating additional long-range contributions from thehyperbolic tip shank, which are included in METALTIP but absent in the idealized sphere–planegeometry.Alt text: Log-scale comparison of sample-side electrostatic pressure calculated by METALTIPand by the analytical sphere–plane model. The two profiles agree near the apex but deviate atlarger radial distances because METALTIP includes long-range contributions from thehyperbolic shank.17� �� �� � �� ����������� ����������� ����� �������������������������a���
��2� ����������� �!�� ������������������������FIG. S5: Dependence of the parabolic curvature coefficient a on the tip radius R. For each R,the coefficient a was obtained by fitting the calculated F–U curve with F (U) = a(U − U0)2 + F0.The resulting a(R) values (solid circles) were fitted using a linear function, a(R) = mR + b, overthe full range R = 2–100 nm (dashed line). The linear fit yielded a coefficient of determinationof R2lin = 0.99990. The maximum relative deviation from the fit was 16.6%, which occurred atthe smallest radius, R = 2 nm, while the root-mean-square relative deviation over the full rangewas 4.0%. This behavior is consistent with the near-field scaling expected for the tip–sampleelectrostatic response for fixed tip-geometry parameters other than R.Alt text: Parabolic curvature coefficient a plotted as a function of tip radius R. The coefficientincreases almost linearly with R over the range R = 2–100 nm, with small deviations from thelinear fit except at the smallest radius.18���� ���  ��� ��  ������ �#��!���������������������������������
���"���$� ���!���%�������������!���"��!���FIG. S6: Bias dependence of the frequency shift (∆f–U curve) calculated for the same tip–samplegeometry as in Fig. 3 of the main article. Solid circles denote the calculated values, and the dashedline shows a parabolic fit.Alt text: Simulated frequency-shift–bias curve for a laterally homogeneous conducting surface.The calculated ∆f–U data follow a parabolic dependence on bias voltage.19��� ��� � �� ����������������������������� ����������������������FIG. S7: Surface-potential profile used for the single-step model. The surface potential is given byVs(x) =CPDdiff2[1 + tanh(xwinterface)], which describes a smooth step from 0 to CPDdiff centered atx = 0. Here, CPDdiff = 0.5 V and the interfacial width parameter was set to winterface = 1.0 nm.Alt text: Smooth surface-potential step used in the single-step model. The potential changesfrom 0 to CPDdiff = 0.5 V across the interface with a smoothing width of winterface = 1.0 nm.20FIG. S8: Representative scanning electron microscopy images of as-provided tungsten (W) tips(Unisoku), together with single-hyperbola fits to the outer tip profile (orange dashed curves). (a)Example of a tip whose outer profile is well reproduced by a single hyperbolic curve over a broadaxial range from the apex; the fit shown uses R = 10 nm and α = 7◦. (b) Example of a tip forwhich a single hyperbolic curve cannot simultaneously reproduce the near-apex region and theshank profile within the fitting window; the fit shown uses R = 6 nm and α = 6◦. For this tip,further decreasing α does not substantially improve the fit to the shank profile within the fittingwindow.Alt text: Scanning electron microscopy images of tungsten tips with hyperbolic fits to the outerprofiles. One tip is well reproduced by a single hyperbola, whereas another shows a mismatchbetween the near-apex region and the shank profile.21��� � ����  ����� ����  ��� � ����  ���� ���  �� �
���� �����������������������������������������!�������Δf���"�α=5 ∘α=1∘ ∘α=2∘ ∘FIG. S9: Simulated ∆f–U curves calculated for different full shank opening angles, α = 5◦, 10◦,and 20◦. The CPD value and bias-voltage range were set to the same conditions as those used forthe comparison with the experimental ∆f–U curve in the main text. The frequency-shift responseshows only a weak dependence on α in this range, supporting the use of α = 10◦ as a representativefixed opening angle in the fitting simulations.Alt text: Simulated ∆f–U curves for full shank opening angles of α = 5◦, 10◦, and 20◦. Thecurves differ only weakly, supporting the use of α = 10◦ as a representative fixed opening angle inthe fitting simulations.22FIG. S10: STM-mode topography and corresponding line profile of the NaCl/Cu(111) interface re-gion used for acquiring the interfacial CPD profile for fitting. (a) Representative STM topographicimage of the measurement area containing a partially NaCl-covered Cu(111) surface. The orangedashed line indicates the position at which the line profile across the NaCl/Cu(111) interface wasacquired. The grayscale corresponds to a height range of 0–400 pm. (b) Line profile taken alongthe dashed line in (a), showing an apparent step height of approximately 320 pm between the bareCu(111) terrace and the NaCl-covered terrace, consistent with a two-monolayer NaCl film.Alt text: STM topographic image and line profile of the NaCl/Cu(111) interface region used forCPD-profile fitting. The line profile shows an apparent height difference of approximately 320 pmbetween the Cu(111) terrace and the two-monolayer NaCl island.23�� ��� � �  �� ���������������������� �����������������������������������FIG. S11: Simulated CPD profiles across a surface-potential boundary for different tip radii, R = 3,5, 7, and 15 nm. As R increases, the interfacial transition becomes broader, whereas sharper tipsyield a narrower CPD transition. This demonstrates that the width of the simulated CPD profileis strongly governed by the electrically effective tip radius.Alt text: Simulated CPD profiles across a surface-potential boundary for different tip radii.Larger tip radii produce broader interfacial transitions, showing that the CPD-profile width ismainly governed by the tip radius.24��  ��� �  �   �� � ������������������������������������������������������
�����������FIG. S12: Simulated CPD profiles across a surface-potential boundary for several tip–sampleseparations, s = 0.7, 0.8, 0.9, and 1.0 nm. Over this physically reasonable range for STM-moderegulation, the CPD-profile shape shows no appreciable dependence on s, indicating that theinterfacial broadening is governed primarily by the electrically effective tip radius rather than bythe operating separation.Alt text: Simulated CPD profiles for tip–sample separations from s = 0.7 to 1.0 nm. The profileshape changes very little over this range, indicating that the interfacial broadening is much lesssensitive to s than to the effective tip radius.25��  ��� �  �   �� � 
����������������������������
����α=5 ∘α=1∘ ∘α=2∘ ∘FIG. S13: Simulated CPD profiles across a lateral surface-potential step calculated for differentfull shank opening angles, α = 5◦, 10◦, and 20◦. The potential difference between the two regionswas set to the experimentally determined CPD difference, and the tip–sample geometry parameterswere fixed at R = 3.0 nm and s = 0.9 nm, except for α. The CPD-profile width changes only weaklyover this range of α, indicating that the interfacial broadening is much less sensitive to the shankopening angle than to the effective tip radius R under the present small-separation conditions.Alt text: Simulated CPD profiles for full shank opening angles of α = 5◦, 10◦, and 20◦. TheCPD-profile width changes only weakly with α, indicating that the broadening is less sensitive tothe opening angle than to the effective tip radius.26� � � �  �����������R�����������������������������������������FIG. S14: Root mean square error (RMSE) between the experimental and simulated CPD profilesplotted as a function of the tip radius R. The simulated CPD profiles were calculated for fixed s =0.9 nm and α = 10◦, and compared with the experimental profile after applying the same interface-position correction as used in the main-text fitting. The RMSE exhibits a shallow minimum aroundR = 2.5–3.0 nm. The dashed line marks the 20%-above-minimum tolerance level, correspondingto an acceptable range of approximately R = 2.0–4.0 nm.Alt text: Root-mean-square error between experimental and simulated CPD profiles plotted asa function of tip radius R. The RMSE has a shallow minimum around R = 2.5–3.0 nm, and the20%-above-minimum criterion gives an acceptable range of approximately R = 2.0–4.0 nm.27�� ��� � �  �� ���������������������� ������
����������������������������������������������������FIG. S15: Simulated CPD profiles calculated using the flat-surface model for R = 3.0 nm andα = 10◦, without and with a +320 pm offset in the tip–sample separation. The offset correspondsto the apparent height difference between the Cu(111) terrace and the 2 ML NaCl island in theexperiment. Because both profiles were calculated for a geometrically flat sample surface, theoffset calculation should be regarded as a sensitivity test for the neglected height difference ratherthan as a full simulation of the non-flat boundary geometry. Nevertheless, the small differencebetween the two profiles indicates that the flat-surface approximation does not substantially affectthe CPD-profile broadening used for the tip-radius extraction.Alt text: Simulated CPD profiles calculated with and without a +320 pm offset in thetip–sample separation. The two profiles differ only slightly, indicating that the neglectedapparent height difference between Cu(111) and two-monolayer NaCl does not substantiallyaffect the CPD-profile broadening.28���� ��
� ��
� ��������������������������������������� ��������������"��������!������� �FIG. S16: Root-mean-square error (RMSE) between the experimental and simulated ∆f–U curvesplotted as a function of the static STM-mode tip–sample separation s. The simulated curves werecalculated for fixed R = 3.0 nm and α = 10◦, and each curve was compared with the experimentaldata after interpolation onto the experimental bias points. The RMSE shows a clear minimumaround s = 0.91–0.92 nm. The dashed line indicates a tolerance level 20% above the minimumRMSE, giving an acceptable range of approximately s = 0.90–0.93 nm.Alt text: Root-mean-square error between experimental and simulated ∆f–U curves plotted asa function of static STM-mode tip–sample separation s. The RMSE is minimized arounds = 0.91–0.92 nm, with an acceptable range of approximately s = 0.90–0.93 nm.29�  �� � �� � ������������������ ������������ ������FIG. S17: Additional examples of CPD profiles measured across the NaCl/Cu(111) interface.Although the CPD changed sharply at the interface, some profiles did not become fully flat inregions far from the boundary. One possible interpretation is that these tips remained relativelysharp near the apex but became thicker in the outer region than expected from the ideal hyperbolicgeometry assumed in the model.Alt text: Additional experimental CPD profiles measured across NaCl/Cu(111) interfaces. Someprofiles show a sharp interfacial change but do not become fully flat far from the boundary,suggesting deviations of the actual tip shape from the ideal hyperbolic model.[1] R. M. Feenstra, Journal of Vacuum Science & Technology B 21, 2080 (2003).[2] R. M. Feenstra, S. Gaan, G. Meyer, and K.-H. Rieder, Physical Review B 71, 125316 (2005).[3] R. M. Feenstra, Y. Dong, M. P. Semtsiv, and W. T. Masselink, Nanotechnology 18, 044015(2007).[4] Y. Dong, R. M. Feenstra, M. P. Semtsiv, and W. T. Masselink, Journal of Applied Physics103, 073704 (2008).[5] N. Ishida, K. Sueoka, and R. M. Feenstra, Physical Review B 80, 075320 (2009).[6] S. Gaan, G. He, R. M. Feenstra, J. Walker, and E. Towe, Journal of Applied Physics 108,114315 (2010).30[7] B. Fornberg, SIAM Review 40, 685 (1998).31