The theory of scanning quantum dot microscopy

Electrostatic forces are among the most common interactions in nature and omnipresent at the nanoscale. Scanning probe methods represent a formidable approach to study these interactions locally. The lateral resolution of such images is, however, often limited as they are based on measuring the force (gradient) due to the entire tip interacting with the entire surface. Recently, we developed scanning quantum dot microscopy (SQDM), a new technique for the imaging and quantification of surface potentials which is based on the gating of a nanometer-size tip-attached quantum dot by the local surface potential and the detection of charge state changes via non-contact atomic force microscopy. Here, we present a rigorous formalism in the framework of which SQDM can be understood and interpreted quantitatively. In particular, we present a general theory of SQDM based on the classical boundary value problem of electrostatics, which is applicable to the full range of sample properties (conductive versus insulating, nanostructured versus homogeneously covered). We elaborate the general theory into a formalism suited for the quantitative analysis of images of nanostructured but predominantly flat and conductive samples.


Electrostatic potentials at the nanoscale and their measurement
Electrostatic forces are among the most common interactions in nature. While they appear in the macroscopic world only when excess charges are present, they are omnipresent at the nanoscale because the constituents of matter, electrons and nuclei, carry discrete charges. These fields significantly influence microelectromechanical systems [1] as well as nanoelectronic components, for example as built-in interface potentials [2,3] or via unwanted background charges [4]. A common method to image and quantify electric potentials on surfaces is Kelvin probe force microscopy (KPFM) which measures the potential difference between a surface and a probe [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The measurement principle of KPFM is derived from the classical Kelvin probe insofar as the electrostatic interaction between two objects, tip and surface, is minimized by application of a compensating bias. An alternative technique, electrostatic force microscopy (EFM) images the tip-surface force (gradient) directly without the compensation mechanism. The lateral resolution of both techniques arises from the distance dependence of the electrostatic force [7,10,19]: The forces at the closest distance between sample and probe, i.e. close to the tip apex, outweigh the forces between the probe and more remote surface areas. Therefore, the potential difference in the vicinity of the tip apex is measured when nulling the overall electrostatic force (KPFM) and likewise the respective forces dominate the EFM contrast. Consequently, surface potential variations can be imaged when scanning the probe over the surface.
However, the lateral resolution of KPFM and EFM is intrinsically limited by this working principle in which the measurement signal is determined by the entire tip interacting with the entire surface. Different tip shapes have been explored but even ultimately sharp tips are not optimal because a shrinking (relevant) interaction area reduces the signal-tonoise ratio [7]. To obtain maximal resolution, the tip-surface distance has to be decreased to a point where chemical forces between tip and sample start acting [12,20]. A rather recent development in this direction is the passivation of a standard metal probe with a single CO molecule [15,16]. However, at these small distances additional effects besides the contact potential difference appear and hamper the quantitative interpretation of the data. Overcoming these limitations is an active field of research [12,16,[20][21][22].
There are alternative approaches which are not based on the detection of forces between the entireties of tip and sample. For classical KPFM, Brown et al suggest the use of a coaxial probe in which the core of a coaxial wire represents the tip to which a potential is applied [23]. The shield strongly reduces the field around this tip and screens long-range forces. Hapala et al derive local electric fields from the distortion of AFM images taken with a CO molecule attached to the tip [21,24]. A more direct technique to improve high-resolution imaging at small tip sample distances was introduced by Lee et al [22]. It utilizes tip-enhanced Raman scattering to measure changes in the vibrational frequency of the stretch mode of a tipattached CO molecule. Due to the vibrational Stark effect this change is directly correlated to the electric field at the position of the CO molecule. Like scanning quantum dot microscopy (SQDM) this technique reduces the sensitive region of the probe to a single molecule. The field on other parts of the tip has no influence on the measurement. However, since the transduction mechanism is based on the magnitude of (any) forces influencing the vibration of CO [25], the quantitative interpretation of high-resolution images recorded at close distances is similarly problematic as discussed for KPFM above.
Recently we developed SQDM, a new technique for the imaging and quantification of surface potentials [26][27][28]. SQDM is based on the gating of a nanometer-size tip-attached quantum dot (QD) by the local surface potential. The technique detects changes in the charge state of the QD via noncontact atomic force microscopy (NC-AFM) and (like KPFM) applies and adjusts a compensation voltage V b across the tipsurface junction to maintain the gating condition in which the QD charge state changes. Since local surface potential variations contribute to the gating, their effect has to be compensated by adjustments of V b . The introduction of a dedicated controller for SQDM imaging [29] reduces the effort of SQDM and permits the imaging of large surface areas [28]. Finally, an efficient image deconvolution technique enables the interpretation of SQDM images in terms of electric potentials in the sample surface plane [28].
The conceptional difference between SQDM one one side and alternative techniques for the imaging of electrostatic fields on the other side is the detection process. In SQDM, we aim to measure electric quantities (surface potential, surface dipole density, surface dipole). But we do not quantify them by a direct force measurement that is intrinsically unspecific, in the sense that the electric forces are overlaid by other forces such as chemical or van der Waals forces. Rather, we sense the electric potential via its influence on a QD, whose charging is then registered by the force measurement. Other methods utilizing a single quantum state for imaging are, e.g., NV centres [30]. In effect, this detour via the quantum dot charging provides a filter which makes sure that only interactions of purely electric origin are registered (chemical or van der Waals forces do not lead to the charging of the quantum dot). This filtering is one of the reasons why SQDM is much easier to make quantitative than other approaches, without the need to resort to models of any kind [21].
Here, we derive a rigorous formalism in the framework of which SQDM can be understood and interpreted quantitatively. In particular, we present a general theory of SQDM based on the classical boundary value problem of electrostatics, which is applicable to the full range of sample properties (conductive versus insulating, nanostructured versus homogeneously covered). Moreover, we elaborate the general theory into a formalism suited for the quantitative analysis of images of nanostructured but predominantly flat and conductive samples. This formalism is applicable to, e.g., atomically ordered metal surfaces with point or line defects, including adsorbates such as isolated molecules, molecular films or 2D materials.

The primary and secondary measurands of SQDM
Before presenting an extensive and rigorous discussion of the theoretical foundation of SQDM, we introduce the basic principles and quantities relevant for SQDM. The working principle of SQDM is illustrated in figure 1. SQDM requires a QD firmly attached to the tip apex of an NC-AFM ( figure 1(A)). The levels of the QD are electronically decoupled from tip, and at least one of them is positioned near the Fermi energy of the tip E F . Then, the charge state of that level can be controlled by gating via the sample bias V b . Changes in the charge state of the QD can be detected as abrupt steps in the tip-sample force and thus as dips in the measured frequency shift signal ∆f (figure 1(B)).
One experimental realization of SQDM is the case of a single molecule as the QD with an initially singly occupied level which can either be emptied or doubly occupied by gating [26]. The respective experiments were performed at 5 K and under ultra-high vacuum conditions. The bias voltages (centre of each dip) at which one electron is added (V + ) or removed (V − ) are specific properties of the tip + QD system (referred to as sensor below) and generally increase with decreasing QD size (figure 1(B)) [27,31]. For a given sensor the change in charge state happens at a specific welldefined potential Φ ± QD of the QD, measured with respect to the grounded tip. If the tip is scanned across the surface, the laterally varying surface potential leads to variations in Φ QD . This shifts the charging voltages V ± (figure 1(C)) [26]. Moreover, lateral variations of the gating efficiency which stem from the surface topography will affect the relation between Φ QD and V b and thus also cause V ± to change. V + and V − are the primary measurands of SQDM.
Next, we derive expressions which connect these primary measurands with properties of the sample surface that are of interest. The general expression for the potential at the QD is where Φ * denotes the contribution from the electrostatic potential distribution on the sample surface, αV b the contribution from the applied bias and Φ T the contribution from local potentials on the tip surface. The latter is constant for a given sensor. The gating efficiency is defined as Equation (1) will be rigorously derived in section 3.4 (equation (48)). Similar to KPFM, SQDM also measures changes in surface potential and no absolute surface potentials (or work functions). Hence, one has to select a reference position r 0 of the sensor above a region of the sample surface in which the surface potential is constant. At this reference position we define the surface potential as zero, which likewise defines Φ * (r 0 ) ≡ 0. We denote the primary measurands at r 0 as V ± 0 ≡ V ± (r 0 ), and the gating efficiency at r 0 as α 0 ≡ α(r 0 ). It is moreover helpful to define a relative gating efficiency as With these definitions we now equate two versions of equation (1), one at r 0 and one at any other sensor position r. Moreover, we do so for the positive and negative charging events, which yields respectively. We can derive the following two expressions via a series of straightforward algebraic transformation from equations (3)-(5), leaving out the variable r for simplicity With these equations we have quantified two main parameters of equation (1) in terms of measured V ± and V ± 0 . It is intuitively clear that the two quantities α rel and Φ * are related to the surface potential and topography. A substantial part of this paper is dedicated to the derivation of this explicit relation. Looking at equation (7), we observe that it only specifies V * ≡ Φ * α . The reason is that Φ * is measured by compensation and the applied compensation voltage acts on Φ QD via α.
We will discuss in the following sections how the secondary measurands of SQDM, V * and α rel , can be derived in a rigorous theoretical framework and how they are related to actually relevant surface properties like the local surface potential Φ s (r).

The imaging formalism of SQDM
SQDM is a scanning probe microscopy. The probe, in our case the sensor consisting of tip and quantum dot, is scanned at a certain height z across the sample surface, with the aim to learn something about the surface. Conceptually, it is opportune to distinguish between the imaging plane at z and the object surface, which is the sample surface itself. Note that since the imaging plane is a plane, while the object surface is an arbitrary surface, the vertical distance between the two is not the same at every lateral position of the probe.
Abstractly, the problem of SQDM can be formulated as follows: Knowing the imaged quantity I in the imaging plane from the measurement, we want to learn about the object quantity O in the object surface. Specifically, the following questions arise: What are the image and object quantities, what is the mathematical relation between I and O, and how can O be determined once I has been measured? We thus have to identify I and O and their relation where Ô is a functional and O(r ) and I(r) are functions defined in the object surface and imaging planes, respectively. In section 3.4 we prove that Φ QD (r), the electrostatic potential at the position of the QD, can be written as is the electrostatic potential in the object surface at r , V b is the bias voltage applied to the sample (object), Φ T is the constant contribution from the potential distribution on the tip surface, and γ is an integral kernel to be determined below. The integral is carried out over the entire object surface. Equation (9) corresponds to an electrostatic boundary value problem B , which abstractly can be written as with the object function O(r ) = Φ s (r ) and the parameter I = V b . Within the measurement protocol of SQDM, B is in fact executed, in the sense that the total surface potential O(r ) − I is adjusted by the parameter I in order to realize a specific constant value of Φ QD (corresponding to the charging potential Φ QD ± ) for all r. From the condition the image function I(r) is generated in the experiment. Mathematically, I(r) = V * (r) is given by The charging events can be detected in ∆f (V b ) curves as sharp dips superimposed on the Kelvin parabola. Shown are ∆f (V b ) curves for three different sensor QD molecules. (C) Series of ∆f (V b ) spectra around V − for a decreasing surface potential beneath the tip + QD sensor during scanning. As a result, V − shifts to more positive values such that the sum of surface potential and applied potential V b remains constant [26].
Equation (12), which can be derived straightforwardly from equation (9) as we will show later (equation (56)), has the form Hence, the inverse of equation (12) is the sought-after relation equation (8) and thus provides a solution to the problem of obtaining the object function from the measured image function of SQDM. In other words, SQDM corresponds to the inversion of a boundary value problem. For a certain class of samples this inversion corresponds to a deconvolution with a kernel γ * . The situation is different for the secondary measurand α rel . With respect to α rel as image function I(r), equation (13) reads (equation (51)) The corresponding object function O(r ) turns out to be the surface topography which determines the function γ(r, r ). However, because the relation between topography and γ is complex, a series of approximations must be made to retrieve the topography from α rel (section 4.4).

Poisson equation and Green's function formalism
As an essential abstraction for the forthcoming discussion we assume that the QD is point-like and that the tip is conductive and electrically connected to the sample at infinity such that tip and sample surface enclose a volume V (figure 2). Under these assumptions the SQDM setup and the imaging problem outlined above (equation (10)) can be formulated as a boundary value problem of electrostatics. The electrostatic potential Φ always and everywhere fulfills the Poisson equation With the help of Green's second identitẙ where S is a closed surface and V the enclosed volume, setting ψ(r ) = 1/|r − r | and φ(r ) = Φ(r ), as well as using the mathematical identity equation (15) can be turned into an integral equation for the potential Φ(r), namely The normal derivatives are defined by where n is the unit vector normal to S and pointing outward, i.e. away from V . Although equation (18) expresses the potential Φ(r) at r in the volume V as the integral over the charge density ρ(r ) in V and the potential Φ(r ) and its normal derivative ∂Φ(r ) ∂n on the boundary S of V (figure 2), this equation does not define a valid boundary value problem for Φ, because specifying arbitrary Φ and ∂Φ ∂n on S does not yield a solution to equations (15) or (18), as solutions are already uniquely defined by either specifying Φ (Dirichlet boundary conditions) or ∂Φ ∂n (Neumann boundary conditions) on S [32]. Specifying both independently overdetermines Φ.
Equation (18) is instructive, however, because it shows that the potential Φ at point r is determined by three contributions (figure 2): (1) the charge density ρ in V , i.e. each point charge at r that contributes to ρ being the source of a Coulomb potential, plus (2) the potential due to a surface charge density σ = 0 ∂Φ ∂n on S and (3) the potential due to a dipole moment density Π = − 0 Φ of a double layer on S . It is easy to show that the potential of a double layer of equal and opposite surface charge densities gives rise to a potential [32] where Π(r ) is the dipole density. Equation (18) applied to the potential at the position of the QD thus shows that there is a relation of either charge density or dipole density and SQDM measurands. Valid electrostatic boundary value problems for solving Laplace or Poisson equations can be formulated by the Green's function formalism [32]. Thereby, either Dirichlet or Neumann boundary conditions can be used. In a Dirichlet problem, the potential Φ is given on the boundary S of volume V (which corresponds to specifying a surface dipole density), while in a Neumann problem the normal derivative of the potential (surface charge density) is given on S . Since the charge density defines an electric field, which is the gradient of the potential, an integration constant in the form of a global potential offset is additionally required in the Neumann problem. In the following we will discuss the implications arising from the two types of boundary conditions. The Green's function specifies the contribution of each point on the surface S to the potential at any point in the volume V . A Green's function for the Laplace or Poisson equations in the volume V bordered by the surface S is a function G S (r, r ) that satisfies the equation Comparing to equation (17), one sees that e 4π 0 |r−r | , i.e. the electric potential at r of a point charge +e at r , is a possible Green's function. However, the most general expression for G S is where F S satisfies the Laplace equation The freedom which F S offers can now be used to implement Dirichlet or Neumann boundary conditions as well as specify the shape of the boundary. With ψ(r ) = G S (r, r ) and φ(r ) = Φ(r ), Green's second identity equation (16) becomes With equation (21) this becomes

Dirichlet problem for a metallic surface
It is possible to make the first term in the surface integral in equation (25) zero by choosing F S such that [32] This is the Dirichlet boundary condition, the superscript D stands for 'Dirichlet'. The physical meaning of Dirichlet boundary conditions equation (26) is clear: The potential G D S (r, r ) at r of a point test charge located at r on the surface S is forced to be zero everywhere in V and on S . This is achieved by external charges outside V which provide an additional potential F(r, r ) at r such that the total potential, including the potential of the test charge, is zero in V and on S . Since the external charge distribution which makes the potential vanish must depend on the position of the test charge, F(r, r ) depends parametrically on r . If S is the surface of a metal, then F can be interpreted as arising from a surface charge density located on S and induced by the test charge at r .
With equation (26), equation (25) becomes If moreover no charges are located in V , we obtain The essence of equation (28) is as follows: If one puts a test charge on a closed and grounded metal surface S , the freely movable charge of the metal will redistribute to zero the potential everywhere on the surface and in the enclosed volume V . This means that G D S (r, r ) = 0 for all r ∈ S . If now the test charge is moved away from the surface along the surface normal into the volume V , an image charge in the metal is created just below the test charge. This image charge comes on top of the charge distribution outside V that creates the potential F. Evidently, the now nonzero potential G D S (r, r ) at r is the potential of a point dipole at the position r . For this reason, the integrand in equation (28), which gives the total potential at r of an arbitrary surface potential distribution on S , is given by 1 of the point dipole per unit charge, times the actual size 0 Φ of the point dipole at r (see equation (20)). Note that any surface potential can be understood as the areal dipole density of a double layer, because at a double layer the potential always has a step (Φ 2 − Φ 1 = Π 0 ). The minus sign in equation (28) follows from the fact that the derivative is taken in the direction of the normal that points away from V , while we have moved the test charge at r into V .
We note that for Dirichlet boundary conditions the Green's function is symmetric, i.e. G D S (r, r ) = G D S (r , r) [32]. This means that source point and field point of the test charge can be exchanged. Hence, an alternative interpretation of Dirichlet boundary conditions is possible (figure 3): If a test charge is placed at r in V , then the total potential, i.e. the potential due to the point charge itself and the potential F due to the external charges, is zero for all r on S . − ∂G D S (r,r ) ∂n then gives the normal component of the electric field on S , which is small near surface depressions, and large near protrusions (figure 3). In this view, the parametric dependence of F on r arises because when the test charge changes its position in V , the external charges have to adjust in order to provide zero potential on S . If S is the surface of a metal, F can be interpreted as the potential due to image charges.

Neumann problem for a dielectric surface
Next, we analyze the situation of a dielectric rather than a metallic surface. In this case, Neumann boundary conditions, formulated in terms of surface charge densities instead of surface potentials, are appropriate. In the present section we recall the case of Neumann conditions on a closed surface S .
For Neumann conditions, F in equation (22) is adjusted such that [32] ∀ r ∈ V , ∀ r ∈ S : where S is the total surface area of S . While it might seem more natural to assume that Here we have used equations (21) and (19). If ∂G N S (r,r ) ∂n is constant on S , then equation (29) follows from equation (30). If equation (29) is inserted into equation (25), then one obtains If there are no charges in V , one has with Neumann boundary conditions The physical meaning of Neumann boundary conditions equation (29) can be understood in the following way: One puts a point test charge on S at r and then adjusts a distribution of additional charges outside V such that the total potential G N V (r, r ) = e 4π 0|r−r | + F(r, r ) at r inside V fulfils equation (29), i.e. changes with a fixed (independent of r ) slope as the test charge is moved along the surface normal n . Moreover, this fixed slope tends to zero as S becomes large. Equation (29) is evidently fulfilled if S is a dielectric surface that has, except at r , a specific areal charge density σ sitting below it, just outside V . The σ is constructed such that moving the test charge from its position on S to just outside V , that is into line with σ, creates the situation of constant potential (zero electric field) everywhere inside V (σ is chosen accordingly; if, e.g. S is a sphere, σ is a constant surface charge density). Before, with the test charge at its original position on S (i.e. slightly inside σ), the electric field was not precisely zero in V , and hence there are small changes ∂G N S (r,r ) ∂n (equation (29)) of the potential at r as the test charge is moved along n . However, since for increasing area of S the contribution of the one slightly 'misaligned' charge at r to the total potential in V becomes negligible, the potential change ∂G N S (r,r ) ∂n at r also converges to zero as S grows. With this interpretation of equation (29), it is straightforward to understand the physical essence of equation (33). The integral equation (33) gives the total potential at r of an arbitrary surface charge distribution applied on S (not to be confused with the charge density outside V that produces F). Its integrand is therefore given by 1 e G N S , which is the potential of the surface charge density element at r per unit charge, times 0 ∂Φ ∂n , the actual size of the applied surface charge density at r . The surface charge density is given by the latter expression, because − ∂Φ ∂n is the normal component of the electric field at S , and generally steps in the normal electric field at surfaces are related to surface charge densities (E 2 − E 1 = σ 0 ). The integral in equation (33) is positive, because negative signs due to the definition of the electric field and due to taking the normal derivative in −n direction cancel. In addition to the integral, the term Φ S appears in equation (33). It arises because, even in absence of the applied charge density 0 ∂Φ ∂n , the potential at r, given by G N S (r, r ), does not remain constant as the test charge moves along −n ; as an immediate consequence, the surface integral in equation (25) with in the integrand is not zero. As we have seen in the Dirichlet case, this integral describes the effect of a double layer at S on the potential at r. In the present case of Neumann boundary conditions, this double layer arises from the polarization of the dielectric at its surface S by the charge density that sets up F. Because by construction (equation (29)) is constant, only the average of the double layer potential on S contributes to the potential at r, although the double layer may well be inhomogeneous on S .
With this, we conclude our discussion of the formalism for the description of Dirichlet and Neumann boundary conditions [32]. We now turn to the application of this formalism to the specific case of SQDM.

Introduction
As stated earlier, the central task of SQDM image interpretation is finding the relation between the potential Φ QD at the QD and the properties of the boundary, i.e. tip and sample surface. In this respect, the goal is to actually determine the Green's function for a given measurement situation and invert the corresponding boundary value problem. We will present a formal solution below. However, it is practically always necessary to introduce simplifying assumptions about the nature of the boundary, particularly at the object surface, to solve this task. Clearly, there is not one single preferable approach, but rather a hierarchy of simplifications which can be only partially explored in this paper.
The boundary at the object surface is fully specified by the boundary shape and the material of which the boundary consists. Depending on the material, appropriate boundary conditions can be chosen: Dirichlet conditions for a conductive surface, Neumann conditions for a dielectric surface. We outline the general case of a boundary which has metallic and dielectric parts in section 3.2. In section 3.3 we introduce the first simplification, namely eliminating Neumann conditions at partially dielectric surfaces via the introduction of an effective dielectric topography of metallic nature. In section 4 we then apply various restrictions to the shape of the boundary, the strongest one being the assumption of a completely flat surface. Subsequently, we relax this assumption and sketch a possible approach to obtain the Green's functions for nonplanar surfaces.

Formal solution
In this section we consider the most general situation that may be encountered in SQDM. While the tip is always metallic, the sample surface may consist of several elements: 1. an arbitrary geometric topography of a metallic sample surface, 2. a dielectric layer of arbitrary dielectric properties and topography (not necessarily closed) on top of the metal surface; this includes the case of a thick dielectric on a metallic sample plate.
An example is given in figure 4. We split the closed surface S into two parts, the dielectric portion D and the metallic portion M . The latter also comprises the tip. In figure 4, we also could consider the surface S which is fully metallic. However, then in equation (27) the volume integral is not zero. Since in SQDM we must solve a problem of the kind defined in equation (13) for 2D image and object functions, it is therefore more advantageous to consider the inner surface S which consists of D and M . This situation can be handled by choosing Dirichlet boundary conditions on the metal, and Neumann boundary conditions on the dielectric. Mixed Dirichlet/Neumann boundary conditions have a unique solution [32], such that the potential at point r is given by The superscript M stands for mixed boundary conditions. Formally, we can define Then, we can write equation (34) compactly as where we have defined The average in equation (36) is carried out over the surface potential on all dielectric surface parts D . This equation has the form of equation (13), if we consider Φ(r) − Φ D as the image function. Its inversion would provide us, at least conceptually, with the object function of SQDM, which turns out to be a surface potential on the metallic part of the surface and a charge density on the dielectric part. The constant Φ D appearing in the object function is the average potential on the dielectric due to the polarization at its surface, resulting from the Neumann boundary conditions. While equation (36) establishes, for the most general sample, a formal relationship between the SQDM image and object functions, i.e. corresponds to equation (13), its practicality is hampered by the following: 1. the partitioning of the surface into metallic and dielectric parts is not known a priori, 2. γ M S (r, r ) is difficult to determine, and 3. when we adjust the bias voltage applied to the sample to compensate local surface potential variations (section 1.2), we execute a boundary value problem equation (10) on the metallic surface S . For this reason, our primary measurands V ± embody information on S , not on S . As a consequence, when inverting equation (13) to get equation (8), we obtain information on S , not S .
For these reasons, we will now introduce an approximation regarding the type of boundary conditions and boundary shape. We will use this assumption throughout the rest of the paper for the interpretation of SQDM images. Rather than pursuing equation (36), we go back to Dirichlet boundary conditions on the metallic surface S in figure 4 and regard the dielectric parts between D and S as a perturbation.

Dielectric topography
From the solution to the Dirichlet boundary value problem in section 2.2 it is clear that the topography of the metal surface influences the ∂G D S (r,r ) ∂n function, and hence also determines the effect of Φ s and V b on the potential in V . Since any di electric above the metal additionally influences how Φ s and V b on the metal surface affect the potential in V , the influence of the dielectric can be lumped together with the metallic topography into an effective function on the surface T (figure 5). We refer to this effective surface T as di electric topography. T is defined such that the gating efficiency resulting from the hypothetical metallic surface T is identical to the original situation of a mixed surface S , consisting of dielectric surface D and metal surface M . Note that generally T = S . When determining V * in the experiment by executing the boundary value problem equation (10) (i.e. nulling the effect of Φ s on the potential at the QD, as in equation (11)) we cannot distinguish between the metallic topography and the contribution of the dielectric. Hence, the inversion of equation (14) will necessarily yield an effective object function which combines metallic topography and dielectric contrib utions. Moreover, the object function Φ s (equation (12)) is given on the dielectric topography T (see figure 5), and not on S directly.
To illustrate the influence of the dielectric on the potential at the QD, we consider the case of a parallel plate capacitor that is partially filled with a dielectric. Specifically, we calculate how the partial filling with dielectric influences the gating efficiency α (equation (2)) which is the integral over r of the effective ∂G D T (r,r ) ∂n as we will show later. Splitting the capacitor of width z t into two capacitors in series, with widths t and z t − t, and filling the capacitor t by a dielectric with permittivity r , we obtain for the total capacity The gating efficiency is given by [31] such that Noting that C 0 QD/sample = 0Aeff zt , we obtain [31], we can expand the above equation (44) This shows that the dielectric will increase the gating efficiency, in lowest order linearly with its thickness t. Since t and | −1 r − 1| enter as a product, a given increase in gating efficiency can be explained by any combination of (increasing) thickness and (decreasing) permittivity. A metal with r = ∞ yields a minimum thickness. When we apply Dirichlet boundary conditions to the mixed surface S in figure 5, we implicitly assume an effective topography T of metallic nature. Therefore, in the dielectric regions T will lie below S .
Since we will use the assumption of Dirichlet boundary conditions on an effective surface T throughout the rest of the paper, we drop the indices D and T from now on. The respective Green's function and its gradient will hence be simply denoted as G and γ but should be interpreted in the context of this assumption.

Dirichlet solution at the surface given by the dielectric topography
In SQDM we essentially measure the electric potential Φ = Φ QD in V , the vacuum space in the NC-AFM junction where ρ = 0. The surface T consists of the tip and sample surfaces, which are metallic (tip) or effectively metallic (sample) and which we assume connected infinitely far away from the tip apex. Thus T is metallic and closed, and we can apply equation (28) to the SQDM situation. Specifically, the boundary values of the potential are Φ s (r + R) on the surface of the grounded tip (figure 6), and Φ s (r ) + V b on the surface of the sample. Local variations of Φ s arise from the local atomic structure of the tip and nanostructures on the surface. V b is the bias voltage applied to sample. The integration in equation (28) is over the tip and sample surfaces such that we obtain We assume that tip and sample are sufficiently well separated such that nanostructures on the latter do not influence the Green's function on the former and vice versa. Then, the integral over the tip surface yields a constant contribution Φ T which is an integral over the relative positions R between all points on the tip surface and the QD (figure 6), and is thus independent of the lateral or vertical positions of the sensor.
Thus we obtain (46) Here, 'surface' refers to the sample part of T . Defining With this, we have derived, in the context of a boundary value problem via the Green's function formalism, the fundamental equation of SQDM (equation (9) of section 1.3). We briefly analyze the sign of γ. Since n points outward from T , i.e. away from V , and since G(r, r ), as the potential at r of the test charge +e at r (see equation (21)) with G(r, r ) = 0 for r ∈ T , becomes less positive towards the surface, G(r,r ) ∂n must be negative. γ(r, r ), given by equation (47), is therefore positive. For a negative test charge −e at r the sign on the right hand sides of equations (21), (28) and (47) changes, and thus, because also G(r,r ) ∂n changes sign, γ(r, r ) remains positive. γ(r, r ) is therefore always positive.
A further notable property is the decay of γ with distance |r − r | which is generally expected because of the 1/|r − r | term in equation (22), notwithstanding the contribution of F. In fact, as we will see below, Dirichlet boundary conditions and the associated F in equation (22) lead to an even faster (exponential) decay of γ(r, r ) with the lateral distance between r and r . This effectively limits the area of the surface around r for which γ(r, r ) is nonzero to a few hundred nm 2 . For example, if we consider a flat surface, the flatness of the surface needs to be fulfilled only within a ≈10 nm radius around the lateral QD position, since the integrand in equation (48) is practically zero outside this region.
We now return to equation (1), Subtracting equation (48) for the case of Φ s = 0 from equation (49), we find Since this has to be valid for any V b , we can conclude that Φ * = 0 and which, with the definition of α rel (equation (3)) is identical to equation (14). Thus, the electrode geometry T that determines the boundary value problem also determines the gating efficiency α(r) of the quantum dot, through the common int egral kernel γ. Equation (51) allows the definition of a kernel γ * which is normalized if integrated over r . When V b = 0, the comparison of equation (1) with equation (48) shows that Dividing this by α(r), we obtain, because of equation (52), We define and thus finally arrive at This is equation (12). Applying −V * (r) to the sample exactly cancels the effect of the surface potential Φ s on the QD as can be seen by the following argument: We assume that at the sensor position r, the potential is Φ * (r). If then we apply a bias voltage V b = V ± 0 − V * (r) (given by equation (56)) we obtain In other words, additionally applying −V * keeps the potential at the QD at exactly the same charging threshold Φ QD ± (r) = α(r)V ± 0 as if Φ * was not present, irrespective of the local gating efficiency α(r). Equation (57) thus demonstrates that V * , as defined in equation (55) (or (12)), fulfills equation (11) and can thus be identified with the experimentally determined image function I(r) of SQDM. Since, for translationally invariant (i.e. flat) samples γ * (r, r ) = γ * (|r − r |), the inversion of equation (56) to determine the object function Φ s (r) becomes a deconvolution problem.

The role of non-local screening
In section 1.2 we found that SQDM provides two independent secondary measurands, V * and α rel . Within the formalism outlined in section 1.3, these two quantities can be interpreted as image functions. In section 3.4 we have identified the surface potential Φ s as the object function related to V * and derived a way to obtain the former via inversion of equation (56). While we have stated in section 1.3 that the topography of the sample surface is the object function of the second image function α rel , we have not formally derived the respective relation yet. In fact, this relation is not straightforward. The reason is nonlocal screening of the potential by the sample topography, as we will show now.
In section 2.2, we have seen that with Dirichlet boundary conditions the Green's function G D T (r, r ) is the sum of a point charge potential and the function F T (r, r ) (equation (22)). Here we have replaced S by T , following the discussion in section 3.3. F T can be understood as the potential of a charge distribution just outside T which screens the influence of the charge at r on the potential at r. Evidently, this charge distribution and therefore F T must depend on the shape of T , because the screening charges are distributed over the entire region just outside the surface. In other words, the shape of T non-locally determines F T and therefore, via equations (22) and (47), also γ. Consequently, γ and via equations (14) and (51), also α rel is therefore a functional of the surface topography T . This non-locality poses a problem for retrieving T from α rel . To approach this problem, we first express the shape of the surface T by a scalar function t d which then becomes the object function associated to α rel . Then, we derive relations between γ and t d for a series of approximations. These relations will also be important for the recovery of Φ s , since knowledge of γ * is required for the inversion of equation (56).

Infinite planes
In the simplest case for which we can solve the boundary value problem explicitly and thus obtain γ ≡ − 0 /e × ∂G/∂n , tip and sample are approximated by infinitely extended parallel planes. For the sample surface, this will become our general assumption for the rest of the paper, because we will treat topographic features and thus the object function t d as perturbations of this plane. This simplification allows us to recast the spatial coordinates r and r , separating them in vertical and horizontal components. We first define a positive z axis vertical to the sample plane which has its origin at the sample and points towards the tip. Details of the dielectric surface topography can then be described via scalar heights t d (r || ) along z at corresponding lateral positions r || , such that r ≡ r || t d .
We denote the height of the QD above the sample plane as z, 1. Parallel planar surfaces. As the simplest possible approximation we assume a planar tip and a sample topography with t d (r || ) = 0 ∀r || . An important consequence of this approximation is a simplification in the functional dependency of γ pp (the index pp stands for parallel plates) on r || and r || , because the Green's function becomes invariant under lateral translations: Here we have also assumed an isotropic surface, i.e. neglected the atomic lattice structure. Equation (58) makes equation (56) a straightforward 2D convolution with an axially symmetric kernel function. In this case, the kernel γ * pp can be easily calculated. To this end we place a point charge just above the surface at (r || , z c ). This point charge produces an image charge just below the surface. Charge and image charge together produce a minimal deformation Φ s (r || ) = φδ(r || − r || ) of the surface potential. According to equation (56), where z t denotes the tip height and d the tip-QD distance, to obtain, after suitable normalization, the function γ * pp (|r || − r || |, z). The potential can be calculated via an infinite series of image dipoles ( figure  7(A)). Figure 7(B) displays numerically calculated potentials Φ * = Φ QD of this dipole for two z t values, together with an asymptotic expression for large |r || − r || | which clearly reveals a (faster than) exponential decay [33], Note that the symmetry of the Green's function with respect to the positions of test charge and QD which we mentioned in section 2.2 is also present in equation (59), since it contains a product of the two sine functions with interchangeable z and z c . The numerical calculation ( figure 7(B)) and its asymptotic behaviour (equation (59)) provides us with the sought-after γ * pp (|r || − r || |, z). This is, in fact, the point spread function (PSF) of SQDM (for the planar electrode configuration). Although we sense long-range electrostatic fields in order to measure the object function, the joint screening of tip and sample leads to an exponential decay of the PSF and thus puts SQDM in line with scanning tunneling microscopy, where the tunneling probability also decays exponentially with distance. In both cases, the result is a superior lateral resolution, because the influence of objects which are not right beneath the probe is strongly attenuated.
Equation (59) shows that there is a considerable influence of the tip-surface separation z t on the shape of γ pp . However, this is not a problem, because the experimental tip height is precisely measured while we attach the molecular QD to the tip by molecular manipulation [26]. Hence, we can calculate the respective γ pp function for each V * image.

Beyond parallel planar surfaces.
There exists a hierarchy of approximations for the boundary shape which add increasing levels of complexity beyond the model of parallel planes for tip and surface. First, we will discuss the implications of a non-planar tip. If we initially retain the restriction of an axially symmetric tip, equation (58) is left intact. The only difference to the planar tip is that γ axial (|r || − r || |, z) will have a different shape than γ pp . γ axial , which denotes an entire class of functions depending on the precise tip shape, can in principle be obtained for any given axially symmetric tip surface by a finite element simulation. We have done finite element simulations for two tip models. The results are displayed in figure 8. A sharper tip results in a weaker screening of the point dipole at r || (figure 8(C)) and thus yields a weaker decay of γ axial with |r || − r || |. In contrast, a hypothetical tip with a concave apex into which the QD is embedded could enhance the decay compared to γ pp even further and thus increase the lateral resolution of SQDM beyond equation (59).
If we drop the assumption of an axially symmetric tip, equation (58) looses its validity and the dependency of γ on |r || − r || | is replaced by a dependency on r || − r || which still implies translational symmetry, Without axial symmetry of γ, V * images of a point-like object would loose their axial symmetry. However, on a larger scale the exponential decay of γ would still dominate. In fact, we have never observed any significant distortions of V * images of circular objects such as single adatoms which could have been attributed to an irregular tip shape. Nevertheless, distortions in an image of a highly symmetric object could, in principle, be used to gain some information about γ in equation (60). With equation (60), which still retains translational symmetry, the inversion of equation (56) becomes a deconvolution with an non-symmetric kernel function which could, for example, be obtained from a finite element simulation with the non-symmetric tip. As we will see later, equation (60) still allows obtaining surface dipole moments P ⊥ of nanostructures from integration of V * images (i.e. without prior deconvolution and thus without the necessity of knowing the function γ). Hence, P ⊥ is independent of the tip shape, which explains our high reproducibility in the experimentally measured values of this quantity [28].
The most general situation is reached if we drop the assumption of a planar surface. Then, irrespective of the tip shape, neither equation (58) nor equation (60) are valid anymore and γ maintains its full separate dependency on r || and r || . In this case, two situations must be distinguished: 1. A surface which is locally flat, i.e. has flat terraces with a radius which is at least as large as the decay length of γ, can still be considered as completely flat. Due to the rapid decay of γ with |r || − r || | local flatness is already fulfilled if the terrace radius exceeds ≈10 nm. In this case the inversion of equation (56) can still be achieved by deconvolution, one just needs to substitute (example for a symmetric tip) γ pp/axial (|r || − r || |, z) by γ pp/axial (|r || − r || |, z − t d (r || )), where t d (r || ) is the dielectric topographic height of the locally flat region below the QD. Therefore, all that is required is to determine γ axial functions for an entire set of surface-QD distances z − t d (by the methods discussed above) and then select the specific γ axial which has been calculated for the respective z − t d value when deconvolving V * on a locally flat sample region of height t d . For example, this procedure allows for a correct deconvolution of V * for nanostructures on different terraces of a stepped sample surface, where t d varies from terrace to terrace. 2. Under the assumption that t d varies substantially on length scales smaller than ≈10 nm, both axial and translational symmetry are broken and the determination of the Green's function is far from straightforward. In principle, a finite element simulation could solve this problem if the topography is known. In this context, we note that SQDM, in principle, contains information on the dielectric topography in the measurand α rel (see section 4.1). Therefore, we will now discuss the inversion of equation (51) with the goal of estimating both, γ and t d , consistently from experimental α rel data. While t d reveals the dielectric topography, the knowledge of γ potentially improves our recovery of Φ s , via inversion of equation (56), in cases of non-trivial topographies.   us to go beyond the assumption of a flat surface, because for a truly flat surface α rel is unity everywhere and therefore contains no information.
We can again, as in section 4.3.2, separately discuss two situations: Firstly, locally flat surfaces where α rel and thus t d are constant in a sufficiently large area, and secondly, surfaces where this is not the case. Both cases can be rationalized as different limits of treating γ within the framework of highdimensional model representation (HDMR). In the following section we will first introduce HDMR and then discuss various approximations in sections 4.4.3 and 4.4.4.

High-dimensional model representation of γ.
As a consequence of the non-local relation between t d and γ, the latter is a functional Γ of the dielectric topography, i.e.
For simplicity we define t d ≡ t d (r || ). With this, equation (14) becomes The task of recovering the surface topography means to invert α rel (r || , z) for t d . This is a special case of a widespread class of problems, namely finding the relationship between a highdimensional input (t 1 , ... t n ) and the output f (t 1 , ... t n ). In high-dimensional model representation (HDMR), the output f (t 1 , ... t n ) is represented as a hierarchical correlated function expansion in terms of the input variables (t 1 , ... t n ) [34]: The general advantage of HDMR is the high convergence. This is achieved by regrouping the standard Taylor expansion [34].
Within the HDMR framework we can write the functional Γ aŝ Here, we have discretized t d (r || ) on the entire surface to n values t i at the respective positions r ||i , i.e. t i ≡ t d (r ||i ) (figure 9). Note that we do not need to provide the positions r ||i as arguments in equation (64), since the different positions of the discrete t i are encoded in their respective functions f i , f ij etc. Moreover, since a topographic feature cannot screen itself, the cases of r ||i = r || , r ||j = r || , . . . have to be excluded from the sums. This is indicated by a prime at the sums in equation (64).
An intuitive physical interpretation of the terms in equation (64) is possible if we define the (flat) part of the sample surface in a reference area (in which α 0 has been defined) as the origin of the z′ axis and all heights t d accordingly ( figure 9). Then, f 0 is the contribution of the topographic object of height t d at r || on an otherwise completely flat surface with and so on. The second term in equation (64) includes first order screening effects where the contribution of the feature at r || is modified by the presence of each of the non-zero topography elements t i individually. Likewise, the third term describes how two topographic elements t i and t j collectively screen the feature at r || . A quantitative example of first order screening is shown in figure 10.
The aspect of non-locality means that all terms beyond f 0 are, in principle, needed for the correct description of Γ , but at the same time we expect their contributions to decrease with increasing number of parameters. Therefore, we ignore terms beyond f i and proceed in two steps. First, we introduce the zeroth-order approximation in which we only include f 0 . Then, we introduce the first-order approximation to include non-local screening to lowest order into the formalism.

Zeroth-order approximation. The zeroth-order approximation simplifies equation (64) tô
Inserting this into equation (62) one obtains a local description of the gating efficiency To linearize the integrand of equation (68) in t d , we perform a Taylor expansion of f 0 around t d = 0 and obtain (69) If we insert this result into equation (68), we get By assuming a completely planar surface with t d = 0 everywhere, the second integral in equation (70) vanishes. A comparison of equation (70) (also assuming an axially symmetric tip) with equation (14) then yields f 0 (r || , z, r || , 0) = γ axial (|r || − r || |, z).
Since f 0 (r || , z, r || , 0) is independent of t d , equation (71) is also valid for surfaces with t d = 0. Because of (equation (51) and definition of α 0 ) we obtain the zeroth-order approximation for α rel α rel (r || , z) = 1 If equation (73) is to be employed for determining t d (r || ) for a given α rel , the integral kernel needs to be known. We discuss shape and norm of this kernel separately. In the zeroth-order approximation we can write, in analogy to equation (71), We obtain the shape of γ topo from the consideration that a zerodimensional topographic feature is a polarizable object. In the homogeneous gating field above the otherwise flat sample (as assumed in the zeroth-order approximation) it therefore behaves as a local dipole. This dipole is analogous to a pointlike deformation of the surface potential (see section 4.3.1) δΦ s (r || ) = βV b δ(r || − r ||0 ) on a homogeneous background V b . Inserting this into equation (48), we obtain where we have moreover used equation (58) which is valid because we have abstracted the zero-dimensional topographic feature as a pure dipole potential on a flat surface. The integral over the second term in the integrand (equation (76)) collapses at r || = r ||0 such that we obtain If we instead consider the original case of a bias V b homogeneously applied to a surface with a zero-dimensional topographic feature and Φ s = 0 everywhere, we find (equation (48)) Since (equations (69), (71), (74) and (75)) (80) A comparison of equations (77) and (80) reveals that the integral kernel γ topo which determines the influence of a topographic feature on α is proportional to γ axial , i.e. proportional to the kernel that determines the influence of a local deformation of the potential on Φ QD , Thus, in the present approximation also γ topo has axial symmetry. This is fully confirmed by a finite element (FE) analysis with z = 24 Å , d = 7 Å and a cuboid-shaped nanostructure with a lateral diameter of 3.5Å ( figure 10). The norm of γ topo in the zeroth-order approximation quantifies the influence on α rel of a single zero-dimensional feature on a flat surface. Several methods of determining the norm of ) on an otherwise flat surface. Tip and sample (z = 23 Å , d = 7 Å ) are kept at zero potential except for a potential of 1 V which is either applied to the patch marked with a red line or to the cuboid (blue lines) in the FE simulation. (B) Simulated potential Φ QD (r || ) ∝ γ(r || , z, r ||i ) as the QD is scanned across the cuboid-shaped nanostructure in panel (A). For the blue (red) curve, the nanostructure (patch) at r ||i is biased. The solid red curve is simulated in the absence of the blue cuboid protrusion, i.e. disregarding its screening influence, while the dashed red curve is obtained with the cuboid (at zero potential) included in the FE simulation (the slight asymmetry in this curve is too small to be visible). This exemplifies how nanostructures at the surface screen or enhance the Green's function (see figure 3).
γ topo are conceivable. For example, one could perform an FE analysis as in figure 10. However, due to the simplifications inherent in the zeroth-order approximation, using the norm of γ topo that is correct for the case of a single zero-dimensional nanostructure will under practically no circumstances yield correct deconvolution results for t d for realistic surfaces. Therefore, we propose the use of a different norm for γ topo , constructed such that the correct height t d of extended terraces is recovered from the deconvolution of α rel according to equation (73). We note that this choice inevitably leads to errors for isolated nanostructures as discussed above.
To calculate this norm of γ topo , we consider α rel on two extended surface terraces and denote the α rel functions on the two terraces as α upper rel and α lower rel , respectively. The height difference between the two terraces is t d . Then we take advantage of the fact that This equation is valid because the constant offset t d of the complete dielectric topography from the lower to the upper terrace can be compensated by a corresponding offset of the vertical position z of the sensor. We expand where Because α and thus also α rel increases with decreasing tipsample separation, g > 0. Setting z = z 0 + t d in equation (83), we obtain With equation (82) this becomes If we assume without loss of generality that α 0 is defined on the lower terrace, then by definition α lower rel (z 0 ) = α lower (z0) α0(z0) = 1 and equation (86) becomes On the other hand, according to equation (73), the gating efficiency of a flat terrace of empty surface t d above the reference terrace is (88) Comparing equations (87) and (88), and using equations (74) and (75), we find in the zeroth-order approximation In other words, in the zeroth-order approximation the norm of the integral kernel γ topo is equal to gα 0 (z 0 ), where according to its definition in equation (84), g can be measured by a calibration experiment above the empty surface. Note that g varies from tip to tip [31] and thus needs to be determined for each tip and z 0 separately. Typical experimental g values are in the range of 0.03Å −1 . We thus have with equations (81) and (89) It was our goal to derive an expression for γ which we have achieved now in the zeroth-order approximation. With equations (61), (67), (69), (71), (74), (81) and (90) we arrive at (91) With equation (62), this result leads to a straightforward relationship between α rel and t d , namely which allows for the deconvolution of α rel images to obtain t d (r || ) maps. However, as stated above, the way in which the norm of γ topo was determined in this section is not unique and will inevitably lead to errors for isolated nanostructures. In the following section we will therefore develop one possible first-order approximation of Γ which removes this problem. Figure 11. Norm of γ topo in the zeroth-order approximation. The norm of γ topo (and thus also γ itself) for a nanostructure of height t d is the same, no matter whether the nanostructure is isolated on the surface (panel (A)) or part of an extended terrace (panel (B)).

First-order mean field approximation.
For the firstorder approximation we truncate the HDMR expression equation (64) after the second term. Hence, we need to specify the functions f i . While it is not impossible to calculate the f i (r || , z, r || , t d , t i ) and thus γ for a large number of r ||i and t i values via a set of FE simulations, this is cumbersome.
To simplify the f i , we introduce an effective height t (r || ), taking a weighted average over the t i , with weights that decrease with increasing lateral distance |r ||i − r || |. Because t i =t ∀i, this approach resembles a mean field approx imation. The sum over all f i can be expressed as a single function f , and thuŝ (94) With the introduction of a flat effective topography t we have practically eliminated the difference in complexity between f 0 and f . The term f 0 describes the case t = 0 (f (t = 0) = 0), while f expresses all changes that occur when t takes on nonzero values. Hence, we can combine both terms into a single function f m (r || , z, r || , t,t) which reduces to f 0 for t = 0: Note that both t d and t are functions of r || . Here, f m (m stands for mean field) describes γ for a local topographic feature of height t d at r || on an infinite plane of height t . We expand equation (95) around t d = 0 and t = 0 Since the first two terms of this expansion are identical to equation (69), the third term must contain the non-local screening effects. With the following considerations we bring this term into a more intuitive form. First, we express the derivatives as differences for an infinitely small , and Then, we realize that the two cases t d = ,t = 0 and t d = 0,t = − describe the same topography with the sole difference that the tip-QD separation has increased to z = z + in the latter case. Hence we can write and f m (r || , z, r || , 0, − ) = f m (r || , z + , r || , , 0)

If we insert equations (99) and (100) into equation (98) we obtain
(101) A comparison with equation (97) reveals that the first term on the right hand side is simply −∂f m /∂t d , such that we get Inserting this into equation (96) turns the Taylor expansion into The third term in equation (103) describes the influence of an extended terrace (for which the second term vanishes because t = t d ) on γ, while the second term describes the effect of isolated protrusions located on a flat terrace of height t . Equation (103) therefore eliminates the problem present in the zeroth-order approximation (equation (92)), where a single term was used to describe both, isolated topographic features and terraces. For situations of several nanostructures which screen each other, equation (103) interpolates between the two limits of a single nanostructure and a terrace. In this respect, the first-order approximation is the simplest possible consistent formulation for γ. The first and the third terms of the sum in equation (103), in combination, are simply γ axial as calculated for the surface-QD distance z −t instead of z, such that we obtain The second term in equation (104) describes a single protrusion (or depression) of height t d −t on a flat surface. Hence, by analogy with equation (74) ∂f m /∂t d has the shape of γ topo (and γ axial , equation (81)), while its amplitude in the first-order mean field approximation remains to be determined (equation (89) has been derived in the zeroth-order approx imation and is not valid here). Therefore, equation (104) simplifies to (105) For clarity, in this equation we have made the dependence of both t d and t on r || explicit. Inserting this into equation (62), we obtain a description of the gating efficiency in the firstorder mean field approximation, While in the second integral the dependence on t (r || ) is explicit, in the first term it is still implicit in the dependence of γ axial on z −t(r || ). For a constant t , the first integral, i.e. the norm of γ axial , is given by (equation (87)) but for an arbitrary function t (r || ) the value of the first integral in equation (106) is unclear. We introduce a function γ * axial the norm of which is always α 0 (z), i.e.
irrespective of the function t (r || ), and write in analogy to equations (107) and (91) Clearly, for constant t , equation (107) follows from equation (109), and therefore equation (109) yields the correct limit. Note that the argument z −t(r || ) in γ * axial is still important (unlike in equation (91)), since it influences the shape of γ * axial . Specifically, γ * axial becomes narrower for larger t values. Using equation (109) in equation (106), we obtain This equation connects the object function t d (r || ) to the image function α rel (r || ) in the first-order mean field approximation. It replaces equation (92) which is valid in the zeroth-order approximation. Note that, not surprisingly, the first term in equation (110), describing the contribution of the flat effective topography, is analogous to equation (92), while the second term quantifies the contribution of local protrusions (or depressions) around the flat effective topography.
It is important to note that the introduction of the first-order approximation breaks the axial and translational symmetry, as was already noted in section 4.3.2. These symmetries are expressed by the fact that γ depends on (|r || − r || |, z) instead of (r || , z, r || , t d ) (equation (58)), which implies Because axial and translational symmetry are broken, this equation is not valid for the γ in equation (105).
(112) Figure 12. Symmetries of γ. (A) For a flat surface, γ depends only on the distance |r || − r || | and has therefore translational and axial symmetry with respect to r || and r || . Its norm is uniquely defined (see text). (B) If a nanostructure at r || breaks the symmetry of the surface, also γ loses its symmetries. Shown are two cuts through γ(r || , r || ) where either r || or r || is held constant (FE simulated curves; scaling different from panel (A)). Encircled in green are the points where the two cuts intersect and which have hence identical γ values. The blue curve illustrates that γ retains its axial symmetry for r || = r || and in the absence of other nanostructures. The red curve illustrates how the local surface potential is screened or enhanced by a nanostructure (see figure 10.) As a consequence of equation (112), while γ topo in the zeroth-order approximation has a uniquely defined norm (see equation (111)), this is not true in the first-order mean field approximation. Therefore, we calculate its integral over r || and refer to it as amplitude A.
Before equation (110) can be inverted to obtain the object function t d from the image function α rel , we need to find expressions for t (r || ) and the amplitude A of γ topo . t expresses the collective screening action of all elements (r ||i , t i ) on γ(r || , z, r || , t d ). We use the sum to determine t . The weighting function θ is given by an analytical formula fitted to a series of FE simulations ( figure 13).
Since equations (93) and (114) are sums of terms each of which contains one t i only, we consider each topography element (r ||i , t i ) in a separate FE simulation and add the terms according to equation (114). To this end, we place a single cuboid-shaped nanostructure (2 × 2 × 1.3Å 3 ) at a series of distances |r ||φ − r ||i | from r ||φ ( figure 13(A)). To obtain the influence of this protrusion on γ(r || , z, r || , t d ), we introduce a minimal deformation Φ s (r || ) = φδ(r || − r ||φ ) of the surface potential at r ||φ while keeping Φ s = 0 everywhere else. Then, the integral in equation (53) collapses to Φ * (r || ) = φγ(r || , z, r ||φ , 0) ∝ γ r ||φ (r || ). Since we have defined t as a sample property, i.e. independent of the position r || of the sensor, the weight θ has to take into account the effect of the protrusion (r ||i , t i ) on γ r ||φ (r || ) at all r || . Therefore, we calculate θ by averaging over r || as where the potential Φ * (r || ) is calculated for the situation shown in figure 13(A) and Φ * 0 (r || ) is calculated in the absence of the protrusion, i.e. describes the situation without screening. The minus sign in equation (115) is needed because ∆Φ * < 0.
The results of the simulation are shown as black dots in figure 13(B). We employ the empirical function (116) to interpolate the FE simulations. Here b is chosen such that the norm of θ is 1; this is required to yield t = t i for a large area where t i = const ∀i. In order to emulate the primed sum in equation (114), which originates from the fact that a topographic feature cannot screen itself and thus the case of r ||i = r || has to be excluded from the sum, the centre of θ is cut out by a Fermi-type damping function. The size of the cut-out region is chosen empirically such that it represents roughly the size of a single atom.
Equation (116) can be used in conjunction with equation (110) and (114) to obtain the object function t d from the image function α rel . Note, however, that a variety of alternative procedures for the determination of t exist and that, in principle, the entire HDMR series in equation (64) can be emulated by an elaborate method to compute t from the t i . This consideration could be the starting point for alternative formulations of the first-order approximation and beyond.
To determine the amplitude A of γ topo in the first-order mean field approximation, we use equation (105) and obtain Integrating equation (117) over r || yields imaging plane We now consider a zero-dimensional topographic feature of height t d at r ||φ on an otherwise flat surface (t(r ||φ ) = 0), see figure 14(A), and evaluate the integrals in equation (118) at r || = r ||φ , such that we get with equations (72) and (111) Here we also have divided by α 0 . Using equation (119), we can calculate the amplitude A =˜γ topo (|r || − r ||0 |, z)d 2 r || (equation (113)). To this end, we calculate the integral on the right hand side of equation (119) from an FE simulation in the presence of the nanostructure for several t d values ( figure 14) and divide by α 0 , which can be calculated as the integral over γ axial (i.e. in the absence of the nanostructure, equation (72)).
In the FE simulation we assume a minimal deformation Φ s (r || ) = φδ(r || − r ||φ ) of the surface potential at r ||φ while keeping Φ s = 0 everywhere else ( figure 14(A)). Since the amplitude of γ topo depends on the size of the nanostructure we use the correct limit of a single atom. The result is plotted in figure 14(B). From a fit to the FE data (red) we obtain its slope c. The amplitude of γ topo is therefore Together with equations (110), (114), and (116) this constitutes the relation between image function α rel and object function t d in the first-order mean field approximation. We note the similarity of equation (120) with equation (89): The amplitude of γ topo in the first-order mean field approximation is proportional to α 0 , as is the norm of γ topo in the zeroth-order approximation. The difference is, however, that the norm in equation (89) has been derived for a nanostructure within an extended terrace (which is equal to a flat surface!), i.e. in the presence of maximal screening, even though the zerothorder approximation disregards screening by construction. This inconsistency creates the problem that was mentioned in section 4.4.3. In the first-order mean field approximation this problem is removed, which is also evidenced by the fact that c g. We obtain c = 1.17 Å −1 , whereas typical g values are in the range of 0.03Å −1 . This comparably high value of c is, however, compensated to a large extent when calculating α in the first-order mean field approximation (figure 15), as becomes obvious in the next section when we discuss the special case of a single isolated nanostructure in the zeroth and first-order mean field approximations.
It is important to note that we perform the linear fit in figure 14(B) over a limited t d −t interval where the FE results are almost linear. The reason for the saturation to −1 at large negative t d −t values is the strong screening of the potential at the bottom of a deep depression in the topography (figure 3). The value of the integral on the right hand side of equation (119) converges to ˜γ (r || , z, r ||0 , t d −→ −∞)d 2 r || −→ 0. When using the first-order mean field approximation, one therefore has to keep in mind that it produces unphysical behavior for large negative t d −t values.

The example of a single isolated nanostructure on a flat
surface. To all orders, the effect of a single isolated nanostructure at r ||0 on the image function α rel can be expressed as (equation (51)) where we have separated the contribution of the nanostructure at r ||0 from the contribution of all other points on the surface. The splitting in equation (121) can be realized in a thought experiment in which we apply the bias voltage V b (equation (9)) either exclusively to the nanostructure (first case, figure 15(B)) or only to rest of the surface (second case, figure 15(C)). The hypothesis, which we are going to prove now, is that the dependence of the two terms in equation (121) on the distance |r || − r ||0 | has the same shape.
If  we assume the tip to be at zero potential (Φ T = 0). The first term in equation (121) thus generates a potential (122) at the quantum dot ( figure 15(B)), while the second term yields (figure 15(C)) Evidently, both Φ QD_1 (r || ) and Φ QD_2 (r || ) exhibit axial symmetry with respect to the point r ||0 ( figure 15(D)), while in equation (123) the γ function for individual r || does not have this symmetry because of the presence of the nanostructure ( figure 15(C)). This is an important difference to the respective functions in the first-order mean field approximation, where because of the mean field t approximation (equation (93), figure 16) each γ exhibits axial symmetry with respect to r || (equation (105)). Next, we apply a bias offset of −V b to the entire tip + surface boundary in the case of equation (123), ( figure 15(E)). This corresponds to a redefinition of the potential without loss of generality. Because the potential of the entire surface excluding the nanostructure is now zero, we obtain Since the offset Φ T has no dependency on r || , this immediately proves our hypothesis, namely that Φ QD_1 and Φ QD_1 and therefore also the two terms in equation (121) do indeed have the same shape as function of the distance |r || − r ||0 |. In other words, the increase in the SQDM image function α rel (r || ) around r || = r ||0 , caused by the increase of γ r ||0 (r || ) due to the nanostructure at r ||0 , has the same shape but opposite sign as the sum of the decreases in α rel (r || ) that originate from the screening action of the nanostructure on γ r || (r || ) at all r || = r ||0 . This is indicated by the red and green arrows in figure 15(A). However, due to the different potential of the tip electrode, the two cases B and C do not cancel precisely to zero, as the FE simulation in figure 15(D) shows. The strong reduction of the red to the black curve in figure 15(D) explains why the large ratio c/g in section 4.4.4 does not lead to a similarly strong enhancement of α rel in the first-order mean field approximation. 4.4.6. Performance of the first-order mean field approximation. We demonstrate the use of the first-order mean field approximation and its deviation from the zeroth-order  (64) are zero in the zeroth-order approximation and consequently the γ r || (r || ) for r || = r ||0 are unaffected by the presence of the nanostructure. In contrast, in the first-order mean field approximation, the sum is nonzero, unless r || = r ||0 whence t i = 0 ∀ r ||i = r ||0 (equation (65), note that r ||i = r ||0 is excluded from the sum). A major difference between the full description of the potential around a nanostructure as captured in FE simulations on the one hand and the first-order approximation in combination with the mean field t approximation on the other hand is the shape of the individual γ r || (r || ) functions in equation (123). In the t approximation axial symmetry around r || is enforced for every individual γ r || (r || ) function. We highlight the consequences of this provision by two exemplary scenarios, (1) a zero-dimensional nanostructure at r ||0 (figure 16(A)) and (2) a surface step ( figure 16(B)). We compare α rel (r || − r ||0 ) as obtained from the zeroth-order and the first-order mean field approx imations with FE simulations for the two scenarios. Figure 16(A) shows that the zeroth-order approximation (green curve) substanti ally underestimates the strength of the influence of the nanostructure on α rel . The obvious reason is that the factor g which determines the norm of γ has been defined for the closed layer with maximal screening. The firstorder approx imation in principle corrects this deficiency as it increases α rel (red curve in figure 16(A)). However, in combination with the mean field t approximation it also results in an unrealistic sombrero shape of α rel (|r || − r ||0 |).
The origin of the sombrero shape can be understood when comparing the individual axially symmetric γ r || (r || ) functions from the first-order mean field approximation with their FE-simulated counterparts (figure 16(C)). To facilitate this comparison, we plot the relative change γ r ||φ (r || )/γ 0 (r || ), where γ 0 (r || ) is the γ r ||φ (r || ) function in the absence of the nanostructure. As before, the potential is applied at r ||φ = 0.
When the nanostructure is placed at zero where also the potential is applied (r ||0 = r ||φ = 0), a scenario with axial symmetry is created as reflected by the red curves in figure 16(C) which are both symmetric around r || = 0. Since t (0) = 0 (no self-screening), γ r ||φ has the same shape as γ 0 but an amplitude that is larger by a factor 1 + ct d . In our example, we have chosen t d = 1 Å which leads to the dashed red curve. The fitting procedure used for the determination of c (figure 14) results in a good correspondence between the amplitudes of the FE simulation (solid red) and of the first-order approximation. Note that the deviations at both ends of the curve are practically irrelevant since the absolute γ values are very small there.
When the nanostructure is placed beside the point r ||φ = 0 where the potential is applied, at r ||0 = 3 Å (blue) and 6Å ) at r ||0 on an otherwise flat surface as calculated with three different methods. The image resolution for the zeroth-and first-order approximation is 1 pixel Å −1 , for the FE simulation it is 3 voxel Å −1 . (B) Profile of α rel perpendicular to a step edge (height 1 Å ). While the zeroth-order approximation comes very close to the FE result, the α rel curve from the first-order mean field approximation is characterized by strong overshooting in the vicinity of the step edge. (C) Plots of the relative enhancement or damping of γ r ||φ (r || ) with r ||φ = 0 as induced by the zero-dimensional nanostructure from panel (A) placed at different positions r ||0 and calculated with the first-order mean field approximation and FE simulations. The reference curve γ 0 (r || ) is calculated in the absence of the nanostructure.
(yellow), it partially screens the potential at r ||φ , thereby damping the respective γ r ||φ . Since the axial symmetry is now broken, the FE simulation curves are asymmetric. The FE simulations reveal that the screening predominantly attenuates γ in the region above and behind the nanostructure (r ||0 r || ), whereas γ is less affected in the opposite direction (r ||0 r || ). However, since for large offsets |r || − r ||φ | both, γ r ||φ (r || ) and γ 0 (r || ), approach zero very quickly, the maximal absolute damping of γ r ||φ occurs in the immediate vicinity of the nanostructure (r || ≈ r ||0 ). This is also where the absolute enhancement of γ r ||0 (r || ) due to the nanostructure is maximal (for r ||φ = r ||0 ). Since both effects are strongest at the same position, namely at r || = r ||0 , both contributions cancel to a large extent and the shape of α rel upon introduction of the nanostructure is barely changed which leads to the results in figures 15(D) and 16(A). This is in marked contrast to the first-order approximation which leads to the sombrero shape.
In the first-order approximation the axial symmetry of the γ r || (r || ) curves is retained due to the t mean-field approximation. Since t > 0 for the dashed blue and yellow curves, the shapes of γ r || (r || ) and γ 0 (r || ) are slightly different which explains the (weak) symmetric r || dependency of the two curves in figure 16(C). The relative amplitude γ r ||φ (r || )/γ 0 (r || ) is now determined by 1 + gt + c(t d −t) (equations (110) and (120)) which is smaller than 1 since t (r ||φ ) > 0 whereas t d (r ||φ ) = 0, and g c. The precise value of the relative amplitude depends on the weighting function θ used to calculate t around the nanostructure ( figure 13(B)) and rapidly approaches 1 as |r ||0 − r ||φ | increases (dashed yellow curve).
Since the screening factor is applied uniformly to γ r ||φ (independent of r || ) in the first-order mean field approximation, too little screening occurs at the position r ||0 of the nanostructure, whereas too much screening is applied in the region around the nanostructure. This creates the regions of α rel < 1 in the red curve in figure 16(A) and, in general, its sombrero shape. It is interesting to note that, for a single nanostructure, the second term in equation (110) which depends on t d -t is zero if integrated over the entire r || imaging plane, since θ (equation (116)) is normalized to 1.
As a second example to assess the properties of the zeroth and first-order approximation we consider a step edge. Figure 16(B) shows that in this case the t approximation and the resulting enforced axial symmetry of γ r || (r || ) has much stronger consequences, namely an over-damping and corresponding overshooting at the lower and higher terrace side of the step respectively. The zeroth-order approximation on the other hand reproduces the FE simulated α rel curve almost perfectly.
In conclusion, we find that of the two approximations evaluated here, the zeroth-order approximation is the better choice for an actual SQDM image deconvolution since it is devoid of overshooting artifacts. Its major shortcoming, a general underestimation of α rel for isolated objects, can be overcome by the first-order approximation, but the additional mean field t approximation as employed here leads to artifacts.
As an outlook, the first-order approximation could, however, become the better choice if the t approximation and the associated axial symmetry of γ is replaced by a more sophisticated approach. Such an approach could, for example, be based on a set of FE simulations as the ones in figure 16(C).

Surface potential and dipole density
Quite generally, the surface potential Φ s can also be understood as a local surface dipole density Π ⊥ ≡ P ⊥ /A relative to the empty surface. This can be seen as follows: If a charge density σ c (r ) is applied at height z c to the metallic sample surface S , this creates an image charge density σ i (r ) below the surface. Together, they form a capacitor. Because the lower 'plate' with charge density σ i (r ) is in the metal, it is grounded, i.e. Φ i (r ) ≡ 0. For the purpose of the present argument, we fix the zero of the coordinate z normal to the surface at this plate. Furthermore, we assume that locally we can approximate the capacitor at r as a parallel plate capacitor. Then, according to Gauss's law the electric potential at z due to the image charge plate is while the potential of the upper plate is given by we obtain the Helmholtz equation In other words, the surface potential is directly proportional to a perpendicular dipole density. For simplicity, we may set = 0 . Equation (130) accounts properly for the orientation of the dipole density. If σ c > 0 (σ c < 0), then also Π ⊥ > 0 (Π ⊥ < 0), and the dipole density points outward (inward), always towards the positive charge, as it must. The work function change ∆W associated with the surface potential Φ s is simply ∆W = eΦ s , where e is the charge of the electron (with sign).

Dipole moments of nanostructures from 2D integration
While the dipole moment density Π ⊥ is a meaningful quantity for extended uniform surface nanostructures on the sample surface, it is of little meaning at Angstrom-sized inhomogeneous structures, for which the surface dipole moment P ⊥ itself gives a more meaningful description. The latter can be obtained from the former by integration over the area N of the nanostructure where Π ⊥ = 0, Therefore, surface dipoles P ⊥ of individual nanostructures can be obtained by SQDM if the respective structure is located on the empty surface such that Φ s = 0 outside the integration area N . While we stick to the Φ s = 0 case in the following, in fact, relative dipole moments can also be obtained if Φ s provides a constant, non-zero background over the entire integration area. An example for the latter case would be the dipole of a vacancy inside an otherwise closed adsorbate layer. We note that under certain assumptions we find Here, we have considered a situation where the condition Φ s = 0 holds even in an entire area A * around the nanostructure which is defined by the range of γ * (white area in figure  17). In this case we can writë Here we have made use of Fubini's theorem to reverse the order of integration. Since Φ s (r || ) does not depend on r || , we find The inner integral of γ * over A yields 1 for any given r || by the definition of γ * . Note that this is only valid if there is a translational symmetry on the surface. In general the integrals ˜γ (r || , z, r || )d 2 r || =˜γ(r || , z, r || )d 2 r || (equation (112)). This means that the normalization criterion of equation (51) does not hold if there is no translational invariance in γ. Thus we finally arrive at¨A This relation proves that dipole moments of nanostructures can potentially be inferred directly from the measured V * image without the need for the deconvolution procedure.

Summary
SQDM is a powerful imaging technique which can reveal the rich electrostatic landscape at the nanoscale with high resolution. The working principle of SQDM is based on the gating of a QD which is mechanically strongly but electronically weakly coupled to the conductive tip of a non-contact atomic force microscope. The potential difference between the QD and the tip (which determines the charge state of the QD) is influenced by the surface potential and the shape of the surfaces of tip and sample. SQDM can hence be interpreted in the framework of a boundary-value problem of electrostatics. In this paper, we have derived the respective general formalism for metallic as well as dielectric boundaries. In SQDM we measure the constantly changing bias voltage which needs to be applied to the sample during scanning to continually compensate the influence of the locally varying static surface potential onto the QD. These are the primary measurands V ± . The secondary measurands V * and α rel are calculated from V ± . The relative gating efficiency α rel quantifies the contribution of the surface topography on the QD potential Φ QD , while V * quantifies the influence of Φ s on Φ QD for a given topography. Consequently, we always have to define a topography before attempting to recover the surface potential from V * .
Since changes in α rel during scanning cannot be uniquely attributed to either metallic or dielectric topographic features, we have introduced the concept of a dielectric topography t d which combines variations in topography and in dielectric properties of the sample surface into a single quantity. The di electric topography of the sample is the hypothetical metallic topography which would yield the same relative gating efficiency α rel as the actual surface topography. With this definition, we have two unknowns Φ s and t d in the object plane of the sample surface which need to be recovered from the two secondary measurands V * and α rel measured in the imaging plane (at the height z of the QD). As the central part of this paper we have outlined a strategy to achieve this recovery. This is an inverse problem in which the point spread function γ plays a central role. γ is related to the gradient of the Green's function on the sample surface (equation (47)) which, in turn, encodes t d . However, this encoding is non-local (section 4.1) such that γ becomes a functional of t d (equation (61)). There exists no analytical expression for γ, even for the Figure 17. Integration areas required to obtain surface dipole moments. To be able to obtain P ⊥ from an integration of V * the following conditions have to hold: Φ s = 0 in N and Φ s = 0 in A * .
This has the consequence that V * = 0 in A ⊂ A * and V * = 0 on the border of A . The sizes of A and A * are given by the range of γ * . simplest case in which tip and surface are approximated as parallel planes. Hence, we have presented a series of approximative solutions and discussed their implications.
Initially, we limited our considerations to a flat sample surface and discussed the cases of a flat tip (section 4.3.1), an axially symmetric tip (section 4.3.2), and an arbitrarily shaped tip (section 4.3.2). All these approximations are convenient since they yield a single γ function for all points in an image (i.e. γ is translationally invariant), which is important for the recovery of surface dipole moments directly from V * images (section 5.2). However, since a flat surface implies that t d = 0 and α rel = const, these approximations cannot be used to interpret α rel images in terms of t d . To go beyond flat surfaces, we have approximated the functional which relates γ and t d (equation (61)) using the mathematical approach of high-dimensional model representation (section 4.4.2) and truncating the resulting sum (equation (64)) after the first or second terms. In the zeroth-order approximation (section 4.4.3), which only includes the first term of HDMR, the nonlocal nature of the relation between t d and γ is completely neglected (equation (67)). In the course of this derivation we have defined γ topo as a measure for the local difference in γ between a flat surface (γ = γ axial ) and a surface with t d = 0 (equation (79)). We found that γ topo is proportional to γ axial such that in the zeroth-order approximation γ differs from its flat-surface pendant γ axial only insofar as its amplitude now increases or decreases with the local value of t d (equation (91)). This yields a straightforward relation between γ and α rel (equation (92)) which can be used to recover the object function t d via deconvolution.
In the first-order approximation (section 4.4.4) we have included non-local effects into our description of the relation between γ and t d . To reduce the complexity of the problem, this was done via a mean field approximation in which the entire topography around each location r || on the surface is averaged into a single number t via a kernel summation (equation (114)). Similar to the zeroth-order approximation, also in the first-order approx imation γ topo is proportional to γ axial . However, with respect to the zeroth-order approximation there are two notable differences: First, γ axial is calculated for a tip-QD distance z −t instead of z and second, γ topo scales with t d −t instead of t d (equation (105)).
In the first-order approximation a protrusion at r || causes an increase in γ at r || (where t d −t is positive) but simultaneously a decrease in the region around r || (where t d −t is negative). These two effects compensate each other to a large extent ( figure 15). In the zeroth-order approximation, on the other hand, there is no influence of the protrusion on its surrounding and vice versa. Thus it is not surprising that the amplitude of γ topo varies much stronger with the height of the protrusion t d in the first-order approximation compared to the zeroth-order approximation. We derive a ratio of c/g ≈ 40 for the strength of this scaling from finite element simulations ( figure 14) and experiments. With this information, we can use equation (110) which gives the relation between α rel and t d in the first-order approximation to recover the object function t d via deconvolution.
In the last part of the paper, we compare the performance of the zeroth and first-order approximation to finite-element simulations for two prototypical surface structures: An atomsized protrusion and a step edge ( figure 16). The benchmark is made with the help of simulated α rel profiles for a given t d (r || ) function. We find that the zeroth-order approximation reproduces the shape of α rel very well in both cases. This can be explained because an isolated nanostructure is a case for which the integral over all γ(r || , r || ) with respect to r || has the same shape as γ axial (r || ) ( figure 15). However, the amplitude predicted by the zeroth-order approximation for the isolated protrusion is too small because the factor c which determines the scaling of γ with t d is chosen to yield correct results for an extended layer which is characterized by maximal screening. This also explains the generally good correspondence in its prediction of α rel for the step edge.
The first-order approximation accounts for the scaling of γ with t d in a better way, which is directly reflected in the higher amplitude of α rel predicted for the isolated protrusion. However, the mean field approximation which we have used in conjunction with the first-order approximation causes a sombrero-like prediction for α rel , since it is not capable of locating the source of screening at a certain topographic element, but instead produces a smeared-out isotropic screening effect ( figure 16(C)). Consequently, the predicted screening at the position of the protrusion is too weak, whereas it is too strong in the area around the protrusion. The same reason explains the strong overshooting visible in the α rel profile across the step edge. These results suggest that at the current level of theoretical description of γ as reported in this paper, the zeroth-order approximation should be preferred over the first-order approximation for the analysis of SQDM images and the recovery of t d and Φ s images. The convincing results of such a recovery which are reported in [28] confirm this conclusion.

Outlook
Our work establishes the theoretical description of SQDM on the firm basis of the boundary value formalism and connects the Green's function formalism to practically applicable recipes for the recovery of dielectric topography and surface potential images from secondary measurands of SQDM. Our theoretical analysis and benchmarking also clearly indicates where the analysis can be improved further. Most importantly, we have linked the poor overall performance of the first-order approximation to the inability of the mean field approx imation to describe local screening effects. Here we see the biggest chance for a future improvement of the methodology: A more flexible approximation fitted to a series of finite element simulations could overcome this limitation and further improve our ability to recover dielectric topographies and surface potentials from SQDM images. Other data-driven approaches based on simulations seem also feasible. While it seems implausible at first to simulate all possible topographic structures of a surface, one could, for example, calculate a representative subset and use it in a machine learning approach. This approach is severely simplified by the short-sightedness of SQDM. Because γ drops exponentially with |r || − r || |, any topographic structure outside a region with a radius of about 10 nm becomes irrelevant. A second aspect of SQDM image analysis which is worth further research efforts is the determination of γ functions from experimental data on model systems. Here, high-quality ground truth data for Φ s (and possible t d ) is required which can be obtained by state-of-the-art ab-initio methods.