Momentum‐Space Imaging of Ultra‐Thin Electron Liquids in δ‐Doped Silicon

Abstract Two‐dimensional dopant layers (δ‐layers) in semiconductors provide the high‐mobility electron liquids (2DELs) needed for nanoscale quantum‐electronic devices. Key parameters such as carrier densities, effective masses, and confinement thicknesses for 2DELs have traditionally been extracted from quantum magnetotransport. In principle, the parameters are immediately readable from the one‐electron spectral function that can be measured by angle‐resolved photoemission spectroscopy (ARPES). Here, buried 2DEL δ‐layers in silicon are measured with soft X‐ray (SX) ARPES to obtain detailed information about their filled conduction bands and extract device‐relevant properties. This study takes advantage of the larger probing depth and photon energy range of SX‐ARPES relative to vacuum ultraviolet (VUV) ARPES to accurately measure the δ‐layer electronic confinement. The measurements are made on ambient‐exposed samples and yield extremely thin (< 1 nm) and dense (≈1014 cm−2) 2DELs. Critically, this method is used to show that δ‐layers of arsenic exhibit better electronic confinement than δ‐layers of phosphorus fabricated under identical conditions.


(A) SIMS depth profiling of δ-layer samples
Time-of-flight secondary ion mass spectrometry (ToF-SIMS) measurements were made using an IONTOF ToF.SIMS5-Qtac100 system at the Surface Analysis Facility at Imperial College London.The base pressure was below 10 −10 mbar, with a mass resolution /Δ = 10,000.The ToF-SIMS system was operated in static mode, with a 25 keV Bi + primary ion beam in high current bunch mode (HCBM).For depth profiling, a 250 eV, 20 nA Cs + sputter ion beam was used, with an impact angle at 45° to the sample surface, with an estimated sputter rate of (0.03 ± 0.02) nm/s.These settings were found to maximise the depth resolution and relative sensitivity factor, whilst minimising the atomic mixing and straggle depth [S1], [S2].Depth profiles were made over a 300 μm × 300 μm sputter crater, where the analytical region was gated to the central 50 µm × 50 µm to minimise the artefacts of the non-uniform edges of the sputter crater.The depth of the SIMS craters was measured using a Zygo NewView 200 3D optical interferometer; this determines a uniform sputter rate for each SIMS measurement.
To calibrate the secondary ion intensity to a volumetric density, a standard δ-layer sample with a known δ-layer density was used as reference.
Since the secondary ion yield for a given dopant species depends on the matrix it is embedded in, a normalisation procedure is required for a complete interpretation of the SIMS data.In the literature, accurate quantification has been studied for ultra-shallow dopant layers buried in a SiO2 and Si stack, where either the SiO2 or Si secondary ion signals are used for normalisation [S3].Here, we normalise the dopant signal to the sum of the SiO2 and Si signals, since the δ-layers lie within a convoluted SiO2 and Si matrix.
Figure S1 shows the SIMS depth profiles of identically prepared phosphorous (red) and arsenic (blue) δ-layers buried 2 and 3 nm below the Si(001) surface.In comparing Figure S1(b,c) and Figure S1(d,e), we find that an arsenic δ-layer yields a better dopant confinement in the direction perpendicular to the δ-layer at comparative depths: (i) For the 2 nm deep arsenic and phosphorous δ-layer, a full width at half maximum (FWHM) of (2.0 ± 0.2) nm and (2.2 ± 0.2) nm is measured respectively; (ii) For the 3 nm deep arsenic and phosphorous δ-layer, a FWHM of (2.2 ± 0.2) nm and (2.9 ± 0.2) nm is measured, respectively.We attribute this to the smaller diffusion constant of arsenic in silicon [S4], [S5], which allows arsenic to remain more confined, relative to phosphorous, for equivalent thermal budgets of sample growth.It is important to add that the SIMS measure of the FWHM is an upper-bound estimate at best, as the depth profiles are inherently broadened dues to instrumental parameters, ion beam interactions and element-specific detections [S6], [S7].Thus, we expect the actual δlayers to be more confined than the results in Figure S1.
By integrating the area beneath the dopant signal of the SIMS depth profiles, the δ-layer sheet density can be determined.Here, the density is defined in terms of a monolayer (ML), which conveniently translates to the fractional coverage of dopants on the initial Si(001) surface, where 1 ML = 6.78 × 10 14 cm 2 .For phosphorous and arsenic, we find a δ-layer density of (0.36 ± 0.04) ML and (0.27 ± 0.03) ML respectively; or (2.5 ± 0.3) × 10 14 cm −2 and (1.8 ± 0.2) × 10 14 cm −2 , respectively.This difference in the saturation density is expected due to the differences in the dissociation chemistry of the dopant precursors [S8], [S9].

(B) Hall measurements of δ-layer samples
Eight-terminal Hall bars were fabricated from the phosphorous and arsenic δ-layer samples after the SX-ARPES experiments.Figure S2(a) is a schematic of the final, etched Hallbar mesa.After cleanroom processing, the samples were then electrically characterised at ≈ 10 K with Hall measurements, whose results are shown in Figure S2(b).By fitting to these data, we can quantify the carrier concentration and mobility for each sample.This is shown in Table S1.Note, we only show the data for the 3 nm deep arsenic δ-layer, as the 2 nm arsenic δ-layer sample did not display any transport due to open-circuit contacts; we attribute this to the potential failure of the etching process during cleanroom processing.
We find that the carrier concentration determined from the Hall experiments in Figure S2 is smaller than the δ-layer densities extracted from SIMS.This is expected, as not all incorporated donors provide electrons and there are several mechanisms for dose loss; deactivated donors in interstitials, vacancies or due to surface oxidation.Additionally, it is important to state that the samples were HF etched prior to the deposition of the Al contacts shown in Figure S2(a), which is required to achieve ohmic contact to the δ-layer.This inevitably brings the δ-layer closer to the surface and leads to another regrowth of the surface oxide; hence, a further dose loss of the δ-layers is expected relative to the initial SX-ARPES measurements.Nevertheless, the carrier concentrations remain well above the metal-insulator transition.In particular, since arsenic is known to suffer from clustering at high concentrations [S10], [S11], this further results in incomplete activation, and we attribute this to the smaller carrier concentration of the arsenic δ-layers compared to phosphorous in Table S1.

(C) ARPES of samples with and without a δ-layer along XW
By aligning the ARPES measurement plane to coincide with the ΓX symmetry plane and tuning the photon energy to 380 eV, the   (  ) band dispersion along the XW symmetry line is mapped.Figure S3 shows a comparison of the ARPES data for a Si(001) sample without a δ-layer (Figure S3(a)), with a 2 and 3 nm deep arsenic (Figure S3(b)) and phosphorous (Figure S3(c)) δ-layer.The data immediately show that the sample without a δ-layer has no occupation of the conduction band.This is expected as the lack of any doping implies that the Fermi-level resides within the bandgap.However, both types of δ-layer samples show a bright conduction band pocket at the bulk X-point, which is a clear signature of the metallic nature of these δ-layers.A crucial point is that the photon energy is identical for all the ARPES data shown in Figure S3, so the inelastic mean free path (IMFP) of the emitted photoelectrons is constant.Thus, for a δ-layer buried deeper beneath the surface, we would expect it to yield a weaker spectral intensity, due to the sub-surface nature of the δ-layer.This is observed in Figure S3(b) and Figure S3(c), where the 3 nm deep δ-layers both show a diminished intensity relative to the 2 nm deep δ-layer.This provides a clear indication that the occupied conduction band state originates from the sub-surface δ-layer.

(D) Quantifying the δ-layer carrier density
We discuss here the procedure used to quantify the carrier concentration from the ARPES Fermi surface.All the relevant ARPES data are shown in Figure S4 for both the arsenic and phosphorous δ-layer samples, where   -  slices are taken through the conduction valleys.Note, in the bottom panel of Figure S4(d), the +  valley does not appear as expected and we attribute this to a slight rotational misalignment during the acquisition of the ARPES data.
The total area enclosed by each valley is determined consistently by thresholding the ARPES data at two iso-contour limits; values at which the maximum intensity drops to 25% (upper bound estimate of the enclosed area) and 50% (lower bound estimate of the enclosed area).These two limits are shown in Figure S4(a-d) as shaded regions around each valley; coloured purple for the +  valley (middle panel) and green for the in-plane (±  ,±  ) valleys (bottom panel).An estimate of the carrier concentration,   , is then determined by applying Luttinger's Theorem [S12], which states that the total area enclosed by the Fermi surface is directly proportional to the number density of electrons.A simple expression for this can be derived from the free-electron (Sommerfeld) model, where: (5) Here,   is the -dimensional volume of -space, where in the case for the twodimensional δ-layer system,  = 2.The result of applying Equation ( 5) is shown in the upperpanel of Figure S4(a-d), which shows a breakdown of the carrier concentration for each one of the six valleys.We find that within the experimental uncertainty, there is an equal occupancy amongst each valley, where the sum yields the total carrier concentration of each δ-layer sample.Since the  = 2 state of the Γ-valley cannot be clearly resolved within the Fermisurfaces shown in the middle panel of Figure S4(a-d), the contribution of the 2Γ state is extracted from the fits to Figure S6(a-d).A free-electron like (or circular) Fermi-surface is assumed.In doing so, the carrier concentration estimate slightly increases by 10% and the total carrier concentration for all δ-layer samples is found to lie within (0.88 ± 0.10) × 10 14 cm -2 .
Alternatively, it is known that a more accurate method to determine the total area enclosed by the Fermi surface is to extract the maximum gradient of the energy integrated photoemission intensity from the experimental ARPES data [S13].This so-called gradient method is based on the theoretical many-body definition of the Fermi surface and this analysis was implemented for the 2 nm arsenic δ-layer sample, using the data in Figure S4(a).In doing so, the total carrier concentration was found to be approximately the same; ~15% higher.

(E) Quantifying the δ-layer thickness by fitting the kz response
Here, we compare the two methods of quantifying the δ-layer confinement () discussed in the main manuscript and show that they are in mutual agreement (Figure S5(g)).
The first method is to apply Equation (1) from the main manuscript, where  is the reciprocal of the difference in the longitudinal span of the Γand Δ-valleys (Figure 3).The second method relies on fitting the   response of the data to the convolution of a pair of mean free path (MFP) broadened Lorentzian's, to another Lorentzian whose width is proportional to  −1 (Figure S5(a)).A Lorentzian curve shape is used since the Fourier Transform of an exponentially decaying wavefunction is Lorentzian, and a pair of curves due to the ±  momenta of the initial states that are in superposition and form a standing wave perpendicular to the δ-layer.The TPP-2M formalism for silicon [S14] is used to estimate the photoelectron MFP,   , via the incident photon energy, ℎ; we find   spans 1.2 -2.0 nm for ℎ between 380 -820 eV.Furthermore, we neglect the contribution from the spectrometer energy resolution, which translates into a momentum broadening of ∆ ≈ 0.02Å −1 about the same ℎ range.The initial position of the MFP-broadened peaks is estimated by finding the corresponding Fermi wave-vector,   , of the Δ-valleys in   -  .A least-squares approach is employed to find the optimal value of  required to fit the   response of both the Γand Δ-valleys.
Horizontal (  ) and vertical (  ) cuts of the Fermi-surface are shown in Figure S5(a), with the 1Γ and 2Γ states being identified.We observe that the first quantum well state (1Γ) brackets the second quantum well state (2Γ), and that the intensity of the   line-profile across the Γ-valley is dominated by the 1Γ contribution.Thus, when fitting the   -response of the data, we found near identical results for   intervals that are either localised around 1Γ, or integrated across the whole Γ-valley.Additionally, by taking the Fourier Transform of the probability density function solutions from the Schrödinger-Poisson model of δ-layers derived in Figure S7 (discussed further in Section G), we find that they also provide a good agreement to the   ARPES spectral response.This is identically shown for the Δ-valley in Figure S5(b), where the 1Δ singlet state can be identified.Furthermore, we again find that the Fourier Transform coincides with theoretical expectations, which also appears more localised along   ; the sharper   -response is associated with the 1Δ wave function being broader in real-space (see Figure S7).
Figure S5(c-f) shows the best estimate of  using the model fit approach outlined in Figure S5(a-b) and discussed above.We find that arsenic δ-layers yield a better confinement relative to phosphorous and that the fits to both the Γand Δ-valleys yield consistent values of .Furthermore, Figure S5(g) confirms that either of the approaches discussed above can be used to estimate , as the values are in mutual agreement.This works because the convolution of Lorentzian's yields a Lorentzian whose width is the arithmetic sum of the component widths.

(F) 𝚪𝚪and 𝚫𝚫-band fitting procedures
The procedures used to fit the ARPES data shown in Figure 4 of the main manuscript are discussed here.A summary of the Γ-band fits are shown in Figure S6(a-d) and the Δ-band fits are shown in Figure S6(e-h), for all δ-layer samples.Each panel shows the ARPES data, overlaid with the best fit parabolic bands, where the number in brackets identifies the type of cut that is taken through the ARPES data.There are two cuts used in the fitting of the ARPES data; (i) momentum distribution curves (MDCs) or horizontal cuts (shown above the ARPES data) and (ii) energy distribution curves (EDCs) or vertical cuts (shown to the right of the ARPES data).We find that the Γ-band fits can be deconvolved into two discrete states, denoted as 1Γ and 2Γ, whereas the Δ-band fits are deconvolved into a single state, 1Δ.
In highly doped δ-layers, a strong electrostatic potential is induced that confines electrons to the plane.The high concentration of donor electrons leads to a substantial overlap of their wavefunction and screening of the attractive cores.As a result, the local disorder is smoothed out.It is therefore convenient to ignore the placement of individual atoms and describe the layer by a constant average charge density.The system now exhibits planar symmetries, where the electrostatic potential only varies in the plane-perpendicular direction.The infinite δ-layer is thus described by one-dimensional Schrödinger-like equations: Here,  ∥ ≈ 0.98  corresponds to the longitudinal anisotropic electron mass in silicon and  ⊥ ≈ 0.19  the transverse effective mass given in terms of the bare electron mass   .Note that the equations are not given in terms of the full electronic wavefunction () = ()(, ), but instead of envelope functions   .The complete solutions of the wave equations include rapid Bloch oscillations that can be omitted in the effective mass approximation.The Schrödinger-like equations are re-written as systems of first order linear equations.The eigenvalues   and envelope functions are computed using a root-finding Schrödinger solver and the second order Runge-Kutta method.Note that the model neglects inter-valley coupling.This is appropriate for the current samples, as the valley splitting of the 1Γ band becomes vanishingly small for δ-layers wider than 0.2 nm [S17].
Assuming the electron wavefunctions are separable, the plane-parallel solutions to the wave equation correspond to a two-dimensional electron gas (2DEG).Each of the quantised bands from the plane-perpendicular envelope functions can thus be occupied by electrons obeying the state filling rules of a 2DEG.Taking the two-and four-fold degeneracy of the bands into account, the following set of equations are obtained: where  Γ () =  ⊥ /ℏ 2 and  Δ () = �  ∥  ⊥ /ℏ 2 are the density of states per unit area of the corresponding bands.The Fermi energy is found by the charge neutrality condition, i.e. the total number of electrons   = ∑    equals the total number of donors , and by finding a self-consistent solution for: where   are the number of bands below the Fermi energy.The electrons are distributed to the different bands according to Equation (9).The potential () =   [()] +   [()] is obtained from the resulting electron carrier distribution.The exchange-correlation potential   is computed with the XC functional by Ceperley and Adler and by employing scaled atomic units in the parametrization by Perdew and Zunger [S18].The Hartree potential   is obtained from Poisson's equation by applying the finite difference method.The method maps the problem onto a grid, such that the second derivative of the potential becomes: where  are the indices of the grid.By adopting the matrix representation and imposing the boundary conditions lim →±∞ () = 0, the Poisson equation takes the form: where   is the number of discrete points centred around  = 0. Multiplication with the inverse matrix  −1 yields the Hartree potential.
Unlike the effective mass approximation in Ref. [S16], the donor atoms are not confined to one doping layer but are subject to a sharp Gaussian segregation profile, in line with Ref. [S17].Adding up the opposing electrostatic potentials from dopants and electrons and the exchange-correlation contribution gives the total potential of the system  � ().The calculation is repeated with () =  � () until self-consistency is reached.It occurs when ∫ [ +1 () −   ()] ≈ 0 is satisfied, where  and  + 1 refer to successive iterations.
Figure S7 summarises the results obtained when using the self-consistent Schrödinger-Poisson model of δ-layers.The wave-functions of the occupied Γ-and Δ-band states are shown in Figure S7(a-b), along with their corresponding carrier concentrations.We find an equal occupancy amongst the Γ-and Δ-bands for metallic δ-layers, whose density is of the order 10 14 cm −2 .Figure S7(c) also shows the evolution of the energy minima of each sub-band as a function of the FWHM of the δ-layer; as the δ-layer becomes broader, the effect of energy quantisation becomes smaller.

Table S1 :Figure S2 :
Figure S2: Summary of the magnetotransport measurements of δ-layer samples.(a) Schematic of the Hall-bar mesa created after cleanroom processing.A cross-section of the δ-layer device stack is shown to the left, contacted with aluminium vias.A zoom in of the Hall-bar geometry is shown to the right, where the orientation of the magnetic field is shown in blue, and the measurement terminals for the current and voltages   ,   and   are indicated.(b) Hall measurements with linear fits, allowing the carrier concentration to be determined.All measurements were made at ≈10 K.

Figure S3 :
Figure S3: ARPES measurements along the XW symmetry line for samples with and without a -layer.(a) ARPES data for a Si(001) sample without a δ-layer, showing that there are no occupied conduction band states.The binding energy origin has been set to the valence band maximum.(b)-(c) ARPES data for a 2 nm deep (b) arsenic and (c) phosphorous δ-layer buried in Si(001).The binding energy origin is set to the Fermi-edge, where a zoom-in of the occupied conduction band is shown at two different depths.The colour scale is consistent across the figure and the data is take using a photon energy ℎ = 380 eV.

Figure S4 :
Figure S4: Summary of the SX-ARPES determined carrier concentration for each -layer sample.The dopant species and depth for each sample is labelled and ordered such that each column represents a summary of the data obtained; (a-b) arsenic and (c-d) phosphorous δ-layers buried (a,c) 2 nm and (b,d) 3 nm deep.(top) Breakdown of the carrier concentration for each of the six valleys, where the solid black line shows the average occupancy of each valley.(middle) ARPES iso-energetic slice extracted through   -  for the +  valley at ℎ = 380 eV (bottom) ARPES iso-energetic slice extracted through   -  for the in-plane (±  ,±  ) valleys.The iso-slices are integrated over the binding energy,   , from -250 meV to   .The isocontours are also shown, which show the best estimate for the total enclosed area for each valley.

Figure S5 :
Figure S5: Extracting the δ-layer confinement from the ARPES   -response.(a) Example of the model fit for one Γvalley dataset, whose chi-squared ( 2 ) value is stated: (middle)   -  Fermi-surface; (upper) the   -integrated line-profile which shows the best fit (blue) to the data (red) by convolving the MFP-broadened (orange) and confinement-broadened (yellow) Lorentzian's; (lower)   and (right)   line profile cuts taken at three different points and colour-coded by the arrows.The white lines indicate the -integration windows for the horizontal and vertical line-cuts.In the lower   line profiles, the Fourier Transform (FT) of the probability density function solutions from a Schrödinger-Poisson model of δ-layers (Figure S7) is in good agreement to the data.(b) The same as (a), but for a Δ-valley dataset.The ARPES   -response appears sharper since the 1Δ wave function is broader in real-space (see Figure S7).(c-f) Plot of the δ-layer confinement vs photon energy for the model fits defined in (a) and (b) for all samples.The circle and triangle data-points are from the Γand Δ-fits respectively.The black crosses indicate outlier fits; the   line profiles were noisy or had very poor fits to the data here.The shaded areas represent the uncertainty from the fitted range of values.(g) Comparison of the δ-layer confinement found using Equation 1 in the main manuscript (green) and the model fits (black) shown in (c-f).Both methods provide estimates to the δ-layer confinement that are in mutual agreement.

Figure S6 :
Figure S6: Summary of the ARPES curve fits to the (upper) and (lower) -bands of the 2 and 3 nm arsenic and phosphorous -layer samples.The dopant species and depth for each sample is labelled and ordered such that each column represents a summary of the data and fits for each δ-layer sample.(a-d) Summary of the fits to the Γ-bands along the bulk XW symmetry line at a photon energy ℎ = 380 eV.The purple and blue parabolas are the fits to the ARPES data, showing that the 1Γ and 2Γ states can be deconvolved.The line profiles above the ARPES image show the momentum distribution curves (MDCs) taken at three different binding energies, labelled (1)-(3); the line profile to the right, labelled (4) shows the corresponding fit to the energy distribution curve (EDC), taken through the centre of mass of the Γ-band.(e-h) Summary of the fits to the Δ-bands along the bulk ΓX symmetry line at a photon energy ℎ = 720 eV.The green parabola shows the best fit to the 1Δ-band.The line profile above the ARPES image shows the fitted MDC, labelled (1) and the line profile to the right, labelled (2), shows a similarly fitted EDC, both of which are taken through the centre of mass of the Δ-band.