Nanoscale imaging of mobile carriers and trapped charges in delta doped silicon p–n junctions

Integrated circuits and certain silicon-based quantum devices require the precise positioning of dopant nanostructures, and hydrogen resist lithography can be used to fabricate such structures at the atomic-scale limit. However, there is no single technique capable of measuring the three-dimensional location and electrical characteristics of these dopant nanostructures, as well as the charge dynamics of carriers and trapped charges in their vicinity. Here, we show that broadband electrostatic force microscopy can be used for non-destructive carrier profiling of atomically thin n-type (phosphorus) and p-type (boron) dopant layers in silicon, and their resulting p–n junctions. The probe has a lateral resolution of 10 nm and a vertical resolution of 0.5 nm, and detects the capacitive signature of subsurface charges in a broad 1 kHz to 10 GHz frequency range. This allows the bias-dependent charge dynamics of free electrons in conducting channels and trapped charges in oxide–silicon interfaces to be investigated. Broadband electrostatic force microscopy can be used to non-destructively image n-type and p-type dopant layers in silicon devices with a lateral resolution of 10 nm and a vertical resolution of 0.5 nm.

T he patterned doping of silicon at the nanometre scale, for both classical and quantum applications, brings new diagnostic challenges. For example, the visualization and electrical characterization of dopant structures within a single transistor is difficult due to the reduced transistor size in current complementary metal-oxide-semiconductor (CMOS) technology. This has serious implications for quality control and security in integrated circuit (IC) manufacturing, where tampering of the local dopant profile has been shown to compromise IC integrity 1 . Thus, the development of techniques that can overcome current diagnostic limitations is essential [2][3][4] . In particular, a technique capable of non-destructive, nanoscale measurements, including quantitative electrical characterization, would be valuable for both IC failure analysis 5 and the creation of quantum architectures 6 based on quantum confined dopant nanostructures.
Using scanning tunnelling microscopy (STM) hydrogen resist lithography, both acceptors 7 and donors 8 can be patterned with atomic resolution into dopant nanostructures 6 such as sheets 9 , wires 10 and dots 11 . This atomic-resolution doping has also been implemented in three-dimensional (3D) structures 12 and can be integrated into a CMOS platform 13 . However, STM is relatively insensitive to the subsurface 14 when used to measure such structures. Scanning microwave microscopy (SMM) can be used non-destructively to determine the 3D location of structures once they are buried within the silicon wafer 15 .
In this Article, we show that donor-acceptor interface nanostructures can be fabricated by STM and then quantitatively characterized using broadband electrostatic force microscopy (bb-EFM). This technique can sense the impedance gradient under different local biasing conditions across a broad frequency range from 1 kHz to 10 GHz, and offers a factor of two improvement in lateral resolution and a factor of five improvement in sensitivity compared with SMM. Quantitative finite-element simulations, using a semiconductor drift-diffusion model, allow us to extract the vertical distribution and polarity of the carriers with high confidence and without the need for a special calibration sample. Furthermore, bb-EFM measurements made at varying frequencies allow us to estimate the density of oxide-interface trap states, which determine the performance of the final device 16 .
Unlike destructive imaging techniques-such as transmission electron microscopy (TEM) 17 , focused ion beam scanning electron microscopy 18 , secondary ion mass spectrometry (SIMS) 19 , atom probe tomography (APT) or scanning spreading resistance microscopy (SSRM) 20 -which require milled, bevelled or cross-sectioned samples, EFM is non-destructive. The technique can also be performed in an ambient environment using a standard atomic force microscope (AFM) with a lock-in amplifier and a signal generator. In addition, measurements are performed in tapping mode, which provides an important advantage over SMM where the best resolution requires the more invasive contact mode.

bb-EFM
The bb-EFM set-up is shown in Fig. 1 and its operation is fully explained in the Methods. In brief, a signal generator is connected to a conducting AFM probe. The electric potential V between the probe and the sample deflects the cantilever under an electrostatic force F es following the relation The complex capacitance gradient dC*/dz, which we denote C′, gives information on the local impedance below the tip and thus on the presence of dopants (changes of C′ relative to an offset C′ 0 are denoted ΔC′ = C′ − C′ 0 ). Equation (1)  Integrated circuits and certain silicon-based quantum devices require the precise positioning of dopant nanostructures, and hydrogen resist lithography can be used to fabricate such structures at the atomic-scale limit. However, there is no single technique capable of measuring the three-dimensional location and electrical characteristics of these dopant nanostructures, as well as the charge dynamics of carriers and trapped charges in their vicinity. Here, we show that broadband electrostatic force microscopy can be used for non-destructive carrier profiling of atomically thin n-type (phosphorus) and p-type (boron) dopant layers in silicon, and their resulting p-n junctions. The probe has a lateral resolution of 10 nm and a vertical resolution of 0.5 nm, and detects the capacitive signature of subsurface charges in a broad 1 kHz to 10 GHz frequency range. This allows the bias-dependent charge dynamics of free electrons in conducting channels and trapped charges in oxide-silicon interfaces to be investigated.
by a distance z, C = Aε/z. For a fixed value of A and a dielectric permittivity ε that is independent of z, C′ = −Aε/z 2 . Thus, for given values of A and ε we can determine the distance z between the surface of a medium and a conducting plane produced, for example, by a delta doping layer below the surface. When applying an amplitude-modulated electric field V ac = V ac,0 cos(2πf mod t) cos(2πf carrier t) with a carrier frequency f carrier and amplitude V ac,0 , we detect, with the lock-in amplifier, a cantilever oscillation at a frequency f mod . The electrostatic force producing this deflection is a measure of the local tip-sample capacitance gradient 21 F es,ωmod can then be studied as a function of V bias and f carrier . Note that f carrier is set independently from f mod and is not influenced by cantilever mechanics 21 . As a result, the frequency range available to bb-EFM (1 kHz-10 GHz) is orders of magnitude larger than those typically used in current-sensing techniques like impedance microscopy 22,23 (<10 MHz), scanning capacitance microscopy 24 (1 MHz), scanning nonlinear dielectric microscopy 25 (1 GHz) or SMM 15,26 (1-20 GHz), all of which rely on a very narrowband matching of the electronic detection circuitry to achieve the required electrical sensitivity. Additionally, we note that the capacitance gradient C′ signal offers more spatially localized information 27,28 , and an improved electrical sensitivity, limited only by the thermal noise of the cantilever 29 .

Subsurface imaging and lateral resolution
Structures consisting of atomic layers ('δ-layers') of B and P were patterned onto a Si(100) surface and subsequently encapsulated by a 15-nm-thick epitaxial Si layer, as detailed in the Methods. Figure  1 shows the resulting surface topography. In Fig. 2a, the same AFM tapping mode topography measurement is overlaid with simultaneously acquired EFM data. In Fig. 1, the surface topography is extremely flat (average roughness of <1 nm) and gives virtually no information about the dopant structure location. The encapsulation was therefore successful, and the applied electric field has no measurable impact on the topographic image. However, the superimposed electrostatic ΔC′ image in Fig. 2a displays an excellent contrast between the spatially averaged values of ΔC′ P = 0.40 aF nm −1 and ΔC′ B = 0.25 aF nm −1 , measured above the nanostructured P and B bars, and the surrounding Si substrate background. This corresponds to a signal-to-noise ratio (SNR) of ΔC′ P /C′ Noise = 12 and ΔC′ B /C′ Noise = 7 (see lower profile line in Fig. 2c), six times better than SMM where we obtained ΔC′ P /C′ Noise = 2 with an AFM probe of comparable sharpness 15 . One curious observation from Fig. 2a, explained in detail in the following section, is that the ΔC′ signal intensity from the B δ-layers is lower than for P δ-layers.
The lateral resolution of the EFM was estimated using patterned B bars of decreasing widths, separated by constant gaps of ~60 nm, as displayed in Fig. 2b, and Extended Data Fig. 1 showing STM measurements prior to encapsulation. The first two sets of two bars The electrostatic force (change of capacitance derivative ΔC′(z)) on a conductive AFM probe is measured to obtain the 2D footprint of a buried dopant layer with nanometre lateral resolution. The doping profile in the vertical dimension is obtained by probing ΔC′ as a function of sample bias V bias , probe-sample distance z and measurement frequency f. The finite element model (FEM) shown at the bottom right is used to extract quantitative information. For frequency-dependent measurements, a signal generator delivers an amplitude-modulated signal to the AFM probe with a fixed kHz modulation frequency and a carrier frequency f carrier between 100 kHz and 10 GHz. The induced electrostatic force is recovered through the lock-in amplifier (LIA) by direct amplitude modulation. Electrostatic force images are acquired together with topography in tapping mode.
with widths of 32 nm and 16 nm are easily resolved and a constant level is reached in the middle of each bar. However, the following bars with widths of 2 × 8 nm and 5 × 4 nm are only resolved as peaks, while thinner bars (that is, 4 nm and 2 nm) are not resolved at all. The line patterns on the right half of Fig. 2b, with 30-and 15-nm gaps, are nearly unresolved. Fitting a sum of two logistic functions I to the edges of the 32-nm and 16-nm bars, we calculate a lateral resolution of 2δ = 10 ± 5 nm. This value is surprisingly good for a pattern buried 15 nm below the surface and is approximately five times better than the lateral resolution achieved by SMM (cf. 2δ SMM = 55 nm) 15 using similar samples and tip-apex radii (R apex,EFM = 26.4 ± 0.2 nm and R apex,SMM = 20 ± 1 nm 15 ; for tip calibration see Extended Data Fig. 2c and Methods). A direct comparison of bb-EFM and SMM on the same P patterned sample, presented in Extended Data Figs. 2 and 3, confirms these EFM results (2δ = 13 ± 2 nm) and demonstrates that bb-EFM can clearly distinguish bars separated by 30 nm, while this is impossible with SMM 15 . We expect the various measurement parameters to affect lateral resolution differently in the force-sensing (C′) and current-sensing (C) techniques. Modelling in Extended Data Fig. 4 shows that, while the tip radius has little impact here, the advantage of the C′ measurement increases at larger tip-sample distance. In both models and experiments, the C′ measurement exhibits maxima in both lateral resolution and dopant contrast at a measurement frequency of ~10-100 MHz (Extended Data Figs. 4 and 5).

Quantitative dopant profiling
In addition to determining the lateral position of a patterned δ-layer, it is also important to measure its vertical position and electrical characteristics. In particular, the carrier sign and concentration are key determinants of the operational behaviour of a device. Figure 3c describes the bb-EFM contrast mechanism with a simplified equivalent circuit. For Si regions without a δ-layer, an effective capacitance (air + oxide + depletion layer) in series with the bulk resistance is formed below the tip. With the δ-layer present, this buried, highly conductive sheet acts as a 2D electron gas (2DEG) capacitively coupled to the tip and bulk substrate in series. The AFM probe apex diameter R and distance z from the surface define the local field penetration depth below the surface and the probed sample volume. Application of a d.c. bias between sample and tip attracts or repels carriers depending on their sign, and modulates the depletion capacitance under the tip depending on their concentration 30 . Additionally, in the δ-layer region, we assume that the bias modulates the occupation of δ-layer bands and therefore the 2DEG sheet resistance. Measuring C′ as a function of tip bias V bias and the tip-sample distance z should therefore allow us to determine the carrier type and vertical carrier concentration profile.
We show this quantification for a test sample containing six P δ-layer bars at two different depths and three different, independently verified, dopant sheet densities (ρ sh,1 = (2.87 ± 0.03) × 10 14 cm −2 , ρ sh,2 = (6.23 ± 0.03) × 10 13 cm −2 and ρ sh,3 = (1.86 ± 0.03) × 10 13 cm −2 ) 15 within an As-doped (3 × 10 14 cm −3 ) Si wafer. Assuming a Gaussian δ-layer distribution with σ = 1.5 nm, this corresponds to a dopant concentration peak value of ρ = ρ sh σ −1 (2π) −1/2 , giving ρ 1 = (7.63 ± 0.08) × 10 20 cm −3 , ρ 2 = (1.66 ± 0.08) × 10 20 cm −3 and ρ 3 = (0.49 ± 0.08) × 10 20 cm −3 . The P δ-layers are structured as sketched in Fig. 3a, with their concentration maximum located at depths of h nom = (3.8 ± 0.1) nm and h nom = (8.9 ± 0.1) nm, as confirmed by SIMS measurements 15 . Although the topography of the test sample is flat ( Supplementary Fig. 2), the EFM image in Fig. 3a clearly visualizes the presence of the patterned bars. To extract the dopant profile we acquired C′(V bias ) and C′(z) curves at different positions on the sample. The Si C′(V bias ) curve in Fig. 3b, acquired over an area with no δ-layer, resembles a classical high-frequency capacitance versus voltage curve for n-type dopants, where we find electron accumulation with high capacitance for high positive biases, depletion with reduced capacitance between 0 and 2 V, and a further decrease of the capacitance in inversion. Acquiring the same curve on the P stripes, we find the C′ difference between accumulation and depletion/inversion is much lower compared to Si. Interestingly, the P bars patterned with a low dopant density exhibit two decays, while the highly doped bars, at both 4 nm and 9 nm below the surface, show only one. This is again in agreement with the picture of a highly doped δ-layer with low R sh acting as an electrically floating metallic plane in which the measured C′ is dominated by the geometrical capacitance between the probe and δ-layer. An applied d.c. bias modulates only the small epitaxial Si region between the probe apex and δ-layer, which has little impact on the overall signal, and results in a single decay in the C′(V) curve (Extended Data Fig. 6b). For a low-doped δ-layer, the sheet resistance is high and the geometrical capacitance plays a minor role, and a simpler, 1D picture of a metal-insulator-semiconductor structure applies. Here, the depletion capacitance C si (V) dominates while the increasingly negative applied bias pushes the depletion layer deeper into the substrate. The resulting decays in the C′(V) curve represent variations in the vertical dopant level (Extended Data Fig. 6).
The experimental C′(z) approach curves in Fig. 3d indicate that the δ-layers can only be sensed when the AFM probe is located very close to the surface (z < 30 nm). For distances above 50 nm, most of the applied potential drops off in the air gap between tip and sample, and C′(z) curves acquired with and without δ-layers practically overlap. Consequently, it is critical to carry out the electrical characterization of the δ-layers with low tapping amplitudes. With our microscope, using PtSi-FM probes, amplitudes between 5 nm and 20 nm were optimal.
To quantitatively extract the active dopant profile from our measurements, we developed an FEM that captures the specific geometry and boundary conditions of the experimental arrangement and solves for the potential (shown in Fig. 3c) and carrier distribution (for details see Methods). Using this model, we calculate C′(z,V,f,h,ρ) for a large range of values of the h and ρ parameters. Best-fit h and ρ values for the measured δ-layer are then extracted through comparison of the calculated and experimentally measured C′ values.
This model is significantly more realistic than previous electrostatics-based FEM models 15 , and allows us to interpret bias-dependent measurements. We note that the simulation is fully based on a drift-diffusion model and therefore cannot be used to explain quantum confinement effects that might be especially relevant at low temperatures where medium-and low-doped silicon is not conductive. However, for measurements at ambient conditions, the doped silicon semiconductor model used here offers a good approximation. This can be seen from the calculated conduction band profile of a P δ-layer shown in the inset of Fig. 3c. Here, we indeed observe that only for high dopant densities (3 × 10 21 cm −3 , solid lines) does the effective conduction/donor band energy reside below the Fermi level (E CB − E F < 0) and thus act like a floating metallic plane. A bias potential applied either to the tip or the substrate only modifies the energy profile in the silicon above the δ-layer. By contrast, for lower dopant concentrations (3 × 10 18 cm −3 , dashed lines) the applied bias also modulates the energy profile below the δ-layer.
Finally, for determination of the dopant density profile we calculated C′(V) and C′(z) curves from the FEM model and fit them to the experimental data. For the bare silicon data, the calculated  C′(V) curve fits very well, assuming a tip work function 31 of 5.7 eV and a tip radius of R = 53 ± 2 nm, only slightly higher than the manufacturer-specified nominal tip radius. This agreement with simulation for the bare silicon confirms the conceptual viability of the model. Using these fixed values of the tip radius and work function parameters, the dopant density ρ and layer depth h are then extracted from the P δ-layer region measurements. The results are summarized in  Fig. 7), the sensitivity of the bb-EFM technique to buried dopants is expected to be strongest at medium to high doping levels (>5 × 10 16 cm −3 ), while low and medium dopant densities lead to a decrease in lateral resolution.

Electrical characteristics of carrier types and traps
In addition to the local carrier profile, the carrier type can also be inferred from the C′(V) curve. As shown in Fig. 4a and Extended Data Fig. 8, application of a bias to the tip or the substrate mainly modulates the carrier distribution between the surface and the highly doped δ-layer, while the electron concentration in the δ-layer and below changes only slightly. When negative voltages are applied to n-type δ-doping, electrons above the δ-layer are repelled from the surface, while for positive voltages they are attracted towards the surface. In contrast, for the C′(V) curve obtained on the p-type B pattern, C′ rises for negative voltages where the holes are attracted towards the surface and falls for positive bias where the holes are repelled (Fig. 4c). The C′(V) curve is therefore essential to identify the carrier type and complements the C′ image. Note that we measured the C′(V) curves at 1 MHz to remove any capacitive contribution from interface traps, which become visible at lower frequencies, as discussed in the following. Table 2 displays the results of the δ-layer dopant density and vertical peak position quantification using FEM, as described above. We note that while the P-and B-doped regions show clearly different C′(V) curves, the region containing mixed P and B doping is dominated by the P response and is almost indistinguishable from the pure P region. The n-type (counter)doping seems to completely eradicate the p-type contrast. Although the single-species pattern B and P measured depths h are hardly distinguishable (Table 2), this observation in the mixed-species pattern might be explained by an enhanced segregation of the P atoms towards the surface during Si overgrowth 33 compared to the B atoms. In this way, the P atoms could effectively shield contributions to the response from the B dopants, which is confirmed by the simulated C′(V) curves from different overlapping B and P dopant distribution widths displayed in Extended Data Fig. 9. The hypothesis is further supported by the SIMS profiles in ref. 7 , and would also explain the overall lower signal contrast on B compared to P observed in Fig. 2.
In addition to determining the carrier type, frequency-dependent C′ measurements provide information about imperfections such as oxide-interface charges or interface traps near the surface of the Si. Such imperfections are of great practical significance because they can lead to reductions in semiconductor device performance 16 . Although fixed oxide charges Q ox either repel or attract carriers and thus effectively add to the applied potential in the C′(V) curve, interface traps Q it can trap holes as well as electrons and thus modify the C′(V) spectra in a more complex way 30,34 . These interface traps can be effectively modelled as an additional capacitance C it in parallel with the depletion capacitance C dep (Fig. 4a). Because interface traps can only follow low-frequency alternating electric fields 35 , it has been shown by Castagné and Vapaille using standard, on-wafer capacitance voltage spectroscopy 35 that, by comparing high-and low-frequency C(V)s, the density of interface states can be successfully extracted following the equation where C ox is the oxide capacitance, Q is the elementary charge and C lf and C hf are the low-and high-frequency capacitances, respectively. The terms in parenthesis are just ratios and the capacitance C ox in the first term can be normalized by the area to c ox = C ox /A. Thus, we can adapt this equation to our measurements by exchanging the capacitance variables with our measured capacitance gradient data C′ at low and high frequency, respectively: While the oxide capacitance in the parentheses corresponds to the capacitance gradient measured in accumulation, C′ ac , we calculate the value for c ox assuming an oxide thickness of 2 nm and ε r = 3.8.
In this way, we directly obtain quantitative values of D it from our measurements without the need for FEM. We measured C′(V) curves at four frequencies from 1 kHz to 1 MHz and visualize the signal C(V,f) in Fig. 4e in a contour plot. We use equation (4) to calculate D it (V) from data acquired at high-frequency (100 kHz) and low-frequency (1 kHz) with respect to the relaxation rate of the traps. Trap densities obtained for Si in this way agree with values typically reported in the literature for native Si/SiO 2 interfaces 36 . For the P δ-layer we observe an increase in D it at positive gate voltages, with a maximum of up to D it = 1.5 × 10 12 cm −2 eV −1 (Fig. 4d). This higher density of interface traps above the P δ-layer is related to surface segregation of the P dopants, which increase the lattice mismatch and induce the formation of trivalent Si atoms, which act as amphoteric traps 34 . Interestingly, we observe a significantly lower trap density of D it = 0.7 × 10 12 cm −2 eV −1 above the B δ-layer, which can again be explained by the lower surface segregation of B compared to P.
Finally, from C′(f) spectra acquired at 0-V bias at different sample locations (Fig. 4f), we find that the interface traps mainly affect the capacitance data in the frequency range below 1 MHz. For frequencies above this value, our C′(V) model, ignoring interface traps, fits well to the data. At lower frequencies we see a clear deviation between model and experiment, especially for Si, where we observe the aforementioned behaviour of C′(f) at frequencies below 100 kHz. Only when we consider interface trapping, shown with the dashed lines, does our model converge better with the experimental data. This result both confirms the presence of interface traps and demonstrates the strength of bb-EFM in its ability to detect and quantify such traps.

Conclusions
We have shown that bb-EFM can be used to determine both the 3D position and the carrier density of buried dopant nanostructures within a Si wafer. The technique can simultaneously distinguish regions of atomically thin n-type and p-type doping, as well as distinguish trapped charges from mobile carriers induced by dopants and local gating. The technique operates in tapping mode, is non-destructive and allows pinpointing of the lateral and vertical positions of δ-layers with 10-nm and 0.5-nm resolution, respectively. This is a factor of 5 better than the resolution reported previously for SMM 15 and a factor of 10 better than reported for SCM 37 on similar samples. The improved resolution is due to a high electrical sensitivity of 1 pF m −1 and the ability to sense the local capacitance gradients across a wide frequency range from 1 kHz to 10 GHz. This allows interface traps in regions containing different dopant types to be identified.
In the devices tested here, the bb-EFM data reveal that there is probably a significant segregation of P towards the surface, whereas the B remains well confined-a fabrication problem for P that could be solved by growing a locking layer to limit the segregation 38 . When combined with semiconductor FEM for parameter quantification, bb-EFM could be used as a diagnostic imaging tool for the next generation of classical and quantum devices. Furthermore, the ability to detect the insertion of hardware Trojans 1 , particularly when implemented below the gate level, would be invaluable for the security sector. It has been shown that by changing the dopant polarity of existing transistors, the cryptographically secure function of certain Intel processors can be compromised 39 . This Trojan could not be detected using optical reverse engineering because only the dopant masks are modified. Although non-destructive X-ray tomography of integrated circuits 3,4 can provide full 3D images not provided by the bb-EFM (which cannot look underneath metallic screening layers), the X-ray technique is insensitive to the electrons and holes directly responsible for chip functionality. Thus, the X-ray and bb-EFM techniques are complementary and could be used in tandem for non-destructive chip inspection.

Methods
Sample preparation. The 3D phosphorus sample was fabricated in an Omicron VT-STM system with a base pressure of less than 2 × 10 −10 mbar. An n-type (15 Ω cm, 500-μm thickness) As-doped Si(100) wafer with deep etched registration markers was cut into 9 × 2 mm 2 samples. After a pre-ultrahigh-vacuum chemical clean with acetone and isopropanol, the sample was loaded into vacuum, thermally cleaned and passivated with hydrogen in a standard process 40 . Electrons from an STM tip were used to depassivate a bar-shaped area of 3.5 × 0.5 µm followed by a saturation dosage of PH 3 (0.09 Langmuir (L)). This results in molecular adsorption of PH 3 exclusively in the hydrogen evacuated region 41 . Subsequently, second (0.003 L) and third (0.002 L) bars were lithographically defined on the surface following the same procedure. The three bars with locally different dopant densities were encapsulated with 5 nm of Si by molecular beam epitaxy from a Si sublimation source. Si growth included a so-called locking layer technique 38 that was developed to reduce P dopant segregation. Following this process, the first 10 monolayers of the sample were grown at a low sample temperature of ~60 °C. A subsequent short rapid thermal anneal of 15 s at 500 °C incorporated the P atoms into the surface 32 , making them electrically active. For further Si growth, the sample temperature was kept at 250 °C. A constant deposition rate of one monolayer per minute was used throughout. After reaching the final thickness, a 2-min anneal to 450 °C flattened the surface, making it suitable for high-quality hydrogen passivation and subsequent lithography steps. The three bars in the top layer were patterned and dosed with the same parameters as the bottom layer bars, but rotated by 90°.
The combined B and P sample was prepared on a Si(100) wafer (1.5 Ω cm, 500-µm thickness). The Si surface was cleaned, then passivated with hydrogen in ultrahigh vacuum, using a standard process 40 . Patterns of clean Si were written by depassivating the Si(100)-2 × 1:H surface 42 using the electron beam from an STM tip, as for the other sample, above. The patterned surface was first exposed to a background PH 3 pressure of ~1 × 10 −9 mbar for 2 min, which resulted in adsorption of the PH 3 molecules exclusively within the depassivated regions 41 . Subsequently, the STM was used again to depassivate an additional region, which was then exposed to a background B 2 H 6 pressure of ~1 × 10 −9 mbar again for 2 min. The P and B atoms were incorporated into the surface via a 2-min anneal at 350 °C (ref. 32 ) and subsequently the donor nanostructure was encapsulated with 15 nm of silicon, grown by molecular beam epitaxy from a Si sublimation source. A low sample temperature of ~250 °C during Si sublimation provided both low donor surface segregation and relatively smooth Si crystal growth 43 . Further details on sample preparation are provided in ref. 7 .
Experimental set-up. The experimental set-up for 3D dopant profiling is shown in Fig. 1 and consists of a 5600 AFM (Keysight Technologies) operated in tapping mode and interfaced with a signal generator through a coaxial cable and tip holder capable of delivering high-frequency signals up to the tip (SMM nose cone). Additionally, a bias voltage can be applied to the tip or the substrate. We used a MXG5183B analog signal generator to cover a frequency range from 1 kHz to 10 GHz. A sinusoidal voltage signal with an angular frequency f carrier was applied between a conductive probe and the bottom of the sample, and the force in equation (1) was measured.
We found that PtSi-FM cantilevers (Nanosensors) with a nominal spring constant of 2.8 N m −1 and a resonance frequency of 67 kHz offered the best trade-off among mechanical stability, wear resistance and electrical conductivity. We also tested highly doped Si and diamond probes (Nanosensors), standard Pt-coated Si probes (Nanosensors) and full Pt probes (RMN). The first two tips showed too low conductivity at frequencies above a few MHz, and the full Pt tips did not allow for stable tapping-mode imaging at low amplitudes. The standard Pt-coated tips tended to wear off and modify very quickly, which is disadvantageous for calibration. When the carrier frequency was beyond the mechanical resonance of the cantilever, the cantilever continued to bend statically due to the nonlinear dependence of the actuation force on the applied voltage. Amplitude modulation of the carrier at a low frequency, f mod , of ~1 kHz led to an oscillation of the down-modulated force, which was detected with an enhanced SNR by a lock-in amplifier sensing directly at the modulation frequency f mod or by demodulating the phase shift of the mechanical cantilever oscillation at f mod .
Frequency calibration was carried out as in ref. 44 . In short, an image was acquired on the region of interest in lift mode at three different lift heights from the surface (for example, 20 nm, 30 nm and 50 nm). While the x-scan was active, the y-scan was stopped and the measurement frequency was changed at every new line, producing a frequency spectrum for every lift height. Extrapolation of the scanning distance defined the background of the frequency spectrum, which was used to normalize the data. Measurements were carried out in a dry nitrogen atmosphere.
FEM. FEM was carried out with COMSOL Multiphysics 5.4 (2D axisymmetric, electrostatics + semiconductor module, semiconductor equilibrium, frequency domain perturbation). The simulation geometry for quantification resembled the experimental conditions. The model consisted of a 15-µm-high AFM tip modelled as a truncated cone with a cone angle θ and cantilever extending the cone end by 10 µm. The tip had a spherical apex with radius R and was located at a distance z above the sample. The sample comprised a 250-µm-thick Si substrate (we verified that the exact thickness is not relevant for our calculations) with a constant background doping of 3 × 10 14 cm −3 and 3 × 10 15 cm −3 for the two different samples, respectively. The additionally grown 15-nm-thick Si layer on the surface was modelled as being undoped. The P and B δ-layers were modelled as acceptor and donor doping distributions with a Gaussian shape centred at h nm below the surface with σ = 1.5 nm, and a lateral width corresponding to the area of the patterned δ-layer. To account for the segregation of dopants towards the surface, we modelled the region between tip and doping layer with a doping density 10 times smaller than the peak value of the δ-layer. The probe surrounding was air with the same size as the Si region. Fermi-Dirac statistics were used to obtain the correct results for high dopant concentrations. The available mobility models-Arora and Fletcher Models-were implemented. A Shockley-Read-Hall generation recombination model was used. To achieve accurate results and convergence of the solution, the meshing was set to scale with the size of the involved geometries. In this way the mesh size on the boundary of the apex and membrane was always at least 20 times smaller than the tip apex. The model was meshed from the inside (apex + δ-layer region) to the outside (cone and surrounding area). Outside parts were meshed automatically by Comsol. The downward directed Maxwell stress tensor in z direction was calculated from the small-signal solution on the tip boundaries, which had been set to Terminal. From these data the capacitance gradient was calculated. To determine the tip geometry, a combination of experiments and simulations were carried out 15,45 . Electrostatic force approach curves (z from 2 to 500 nm) were simulated for a wide range of tip radii R and cone angle θ. Using a fitting procedure with the simulated approach curves, we extracted the tip radius and cone angle that led to the best agreement with the experimental approach curve.
After estimation of the tip geometry, the capacitance gradient C′ was calculated for a wide range of dopant concentrations ρ and depth h, changing the measurement parameters (bias V, frequency f and tip-sample distance z, respectively) while keeping the geometry fixed. A look-up table/interpolation function was generated from these data to fit the simulated data to the experimental data. The fitting results are ρ and depth h. Note that the calculated forces from the small-signal solution assume small a.c. voltages of <1 V to satisfy the harmonic perturbation approach. To relate experiments obtained at excitation voltages of 2-3 V with the simulations, an additional integration-differentiation step over the experimental voltage range was necessary.

Data availability
All data needed to evaluate the conclusions in this paper are present in the paper and/or the Supplementary Information. Additional data related to this paper can be requested from the authors. The data created during this research are available at https://doi.org/10.5281/zenodo.3899692.