Spatial charge inhomogeneity and defect states in topological Dirac semimetal thin films of Na3Bi

See allHide authors and affiliations

Science Advances  22 Dec 2017:
Vol. 3, no. 12, eaao6661
DOI: 10.1126/sciadv.aao6661


Topological Dirac semimetals (TDSs) are three-dimensional analogs of graphene, with carriers behaving like massless Dirac fermions in three dimensions. In graphene, substrate disorder drives fluctuations in Fermi energy, necessitating construction of heterostructures of graphene and hexagonal boron nitride (h-BN) to minimize the fluctuations. Three-dimensional TDSs obviate the substrate and should show reduced EF fluctuations due to better metallic screening and higher dielectric constants. We map the potential fluctuations in TDS Na3Bi using a scanning tunneling microscope. The rms potential fluctuations are significantly smaller than the thermal energy room temperature (ΔEF,rms = 4 to 6 meV = 40 to 70 K) and comparable to the highest-quality graphene on h-BN. Surface Na vacancies produce a novel resonance close to the Dirac point with surprisingly large spatial extent and provide a unique way to tune the surface density of states in a TDS thin-film material. Sparse defect clusters show bound states whose occupation may be changed by applying a bias to the scanning tunneling microscope tip, offering an opportunity to study a quantum dot connected to a TDS reservoir.


Three-dimensional (3D) topological Dirac semimetals (TDSs), such as Na3Bi and Cd3As2 (14), express the pseudorelativistic physics of 2D Dirac material graphene but extended to three dimensions. TDSs can yield ultrahigh mobilities (5) and new physics such as the chiral anomaly (6). The close approach of the Fermi energy EF of a Dirac semimetal to the Dirac point ED is limited by fluctuations in the Fermi energy due to charge puddling. Thus, low Fermi-energy fluctuations are critical for the observation of new Dirac point physics such as velocity renormalization (79) and the Dirac plasma (10, 11) at |EFED| < kBT, where kBT is the thermal energy. Potential fluctuations arise in Dirac semimetals when the Fermi energy EF measured relative to the Dirac point energy ED approaches zero, meaning that the carrier density tends to zero and metallic screening disappears. Poorly screened disorder induces spatial fluctuations in the local ED and, hence, fluctuations ΔEF, corresponding to local electron and hole “puddles.” The carriers in puddles, in turn, screen the disorder, with ΔEF determined self-consistently by the disorder and the screening properties of the Dirac materials (12, 13). Puddles have been visualized in graphene using scanning single-electron transistor microscopy (14) and scanning tunneling spectroscopy (STS) (1517), where the fluctuations are largely governed by the underlying substrate, and have also been measured in the Dirac surface state of a topological insulator (18).


Here, we use 20-nm Na3Bi thin films grown via molecular beam epitaxy (MBE) on both semiconducting [Si(111)] and insulating [α-Al2O3(0001)] substrates in ultrahigh vacuum (UHV) to probe ΔEF fluctuations using scanning tunneling microscopy (STM) and STS. Figure 1A shows a large-area (400 nm × 380 nm) topographic STM image of a thin 20-nm film of Na3Bi on Si(111), with several atomically flat terraces >100 nm in size. The inset of Fig. 1A reveals the (1 × 1) Na-terminated atomic lattice (a = 5.45 Å) with an individual Na vacancy. Figure 1 [B (45 nm × 45 nm taken immediately after growth) and C (30 nm × 30 nm taken 7 days after growth)] shows the topography of two atomically flat regions away from step edges or the screw dislocations seen in Fig. 1A. Figure 1D shows an atomically flat region (30 nm × 30 nm) of Na3Bi grown on sapphire [α-Al2O3(0001)]. Whereas atomically flat regions of Na3Bi up to 100 nm × 100 nm can be obtained on Si(111), sparse defect cluster sites (few per 100 nm × 100 nm) give rise to tip-induced ionization ring features (1922) that will be discussed further below. This necessitated focusing on smaller areas free of ionization rings to unambiguously determine the variation in Dirac point.

Fig. 1 Large-area morphology and spectroscopy of 20-nm Na3Bi on Si(111).

(A) Large-area (400 nm × 380 nm) topographic STM image (bias voltage V = −3 V and tunnel current I = 50 pA). Inset: Atomic-resolution STM image with a lattice constant of 5.45 Å (taken on a separate 20-nm Na3Bi film) showing an individual Na vacancy at the surface. (B) STM topography (V = 300 mV and I = 250 pA) on a 45 nm × 45 nm region of Na3Bi. (C) STM topography (V = −550 mV and I = 200 pA) on a 30 nm × 30 nm region of Na3Bi. (D) STM topography (V = −50 mV and I = 100 pA) on a 30 nm × 30 nm region of Na3Bi on sapphire [α-Al2O3(0001)]. (E) Area-averaged scanning tunneling spectra (vertically offset for clarity) corresponding to four different regions of the sample. Black spectra corresponding to the region in (B) and red spectra (topography not shown) were taken immediately after growth, whereas the blue spectra corresponding to the region in (C) and green spectra (topography not shown) were taken 1 week after growth. Feature ED reflects the Dirac point, which corresponds to the minimum in the LDOS, and D represents the resonance feature (discussed in the main text).

Figure 1E shows area-averaged scanning tunneling spectra of the Na3Bi film. STS measures the differential conductance dI/dV as a function of sample bias V, proportional to the local density of states (LDOS) at energy eV relative to EF. STS was performed at points on a grid encompassing the entire regions shown in Fig. 1 (B and C), along with two other regions not shown. Figure 1E shows the averaged spectra of these four areas. Two key features are prominent in all spectra: a distinct minimum in the LDOS corresponding to the Dirac point energy (ED) and a resonant feature ~30 meV below ED, labeled D. As-grown, the Dirac point is located ~20 meV above the Fermi level, indicating p-type doping. Similar doping (~25 meV) has been reported on similar thickness Na3Bi films on Si(111) measured with angle-resolved photoelectron spectroscopy, validating our assumption that the minimum in the LDOS reflects the Dirac point (23). Seven days after growth, the Dirac point has shifted to approximately 15 meV below the Fermi level, reflecting a gradual global n-type doping of the Na3Bi due to the adsorption of atomic species present in UHV. Photoelectron spectroscopy (not shown) shows the emergence of carbon-related species (such as CO and CHx) after several days in UHV. This adsorption results in the formation of intermittent impurity clusters on the surface; hence, all topography and spectroscopic measurements were deliberately performed away from these sites. The resonance feature D is unambiguously tied to ED, and not to EF, because the relative energy shift of D with respect to ED remains unchanged within experimental accuracy during the transition from p-type to n-type doping.

We first turn to the spatial variation of the Dirac point energy, ED, found by tracking the position of the minimum differential conductance in STS (regions B and C) or, alternatively, the shift in the defect resonance D (region A; see discussion provided in section S2). Figure 2 (A and B) shows the spatial variation of ED for regions A and B corresponding to Fig. 1 (B and C, respectively). In both Dirac point energy maps, it can be seen that a clear, continuously connected local potential modulation emerges, correlated on a scale much larger than the crystal lattice or point-spectroscopy grid pixel size (0.5 nm). This modulation in ED represents the puddling of charge density at the surface. The individual Na vacancy sites are shown as black dots in Fig. 2 (A to C). The local variation in ED is found to be positively correlated with the individual Na vacancy sites that have a more positive Dirac point energy (we have verified this statistically; see the Supplementary Materials); this indicates that the Na vacancies act as acceptors and a significant source of charge disorder. In addition, Fig. 2C shows the local ED of 20-nm Na3Bi on α-Al2O3(0001) (labeled region C), which has a larger n-type doping. The upper, middle, and lower panels of Fig. 2D show histograms of ED relative to EF for the scans in Fig. 2 (A to C, respectively). From spatial autocorrelation analysis of the puddle maps, we determine spatial coherence lengths ξ of 13.4 ± 5.2, 9.3 ± 2.4, and 5.1 ± 1.9 nm for regions A, B, and C, respectively (see the Supplementary Materials for calculation details). The observed histograms of ED have mean values of 19.7 ± 1.7 meV (region A), −15.4 ± 1.3 meV (region B), and −47.2 ± 0.6 meV (region C) and SD EF of 5.6, 4.2, and 3.5 meV, respectively. These values are comparable to the 5.4 meV observed for graphene on hexagonal boron nitride (h-BN) (15). However, undersampling results in an underestimate in the magnitude of ΔEF by a factor (L/ξ)2/[(L/ξ)2 − 1], where L is the scan size in each region [also likely in the study of Xue et al. (15)]; hence, the corrected ΔEF values are 6.1 and 4.6 meV for regions A and B, respectively, whereas region C remains essentially unchanged because of the small coherence length.

Fig. 2 Charge puddling profiles of p-type and n-type Na3Bi.

(A) Dirac point energy map of the 45 nm × 45 nm (90 pixels × 90 pixels) region of p-type Na3Bi on Si(111) (V = −250 mV and I = 250 pA), corresponding to region A represented in Fig. 1B. Scale bar, 15 nm. (B) Dirac point energy map of the 30 nm × 30 nm (60 pixels × 60 pixels) region of n-type Na3Bi Si(111) corresponding to region B of Fig. 1C (V = −150 mV and I = 200 pA). Scale bar, 10 nm. (C) Dirac point energy map of the 30 nm × 30 nm (60 pixels × 60 pixels) region of n-type Na3Bi Si(111) grown on α-Al2O3(0001) (labeled region C) (V = −150 mV and I = 200 pA). Scale bar, 10 nm. (D) Upper, middle, and lower panels representing histograms of the Dirac point energy maps in (A) to (C), respectively. The histograms are color-coded to reflect the intensity scale in the corresponding Dirac point energy map. (E) Plot of spatial coherence length as a function of Fermi energy, comparing the experimental data for region A (red triangle), region B (blue triangle), and region C (purple diamond) to theoretical predictions (shaded region) where the upper bound (solid line) is defined using vF = 2.4 × 105 ms−1 and α = 0.069, whereas the lower bound (dashed line) uses vF = 1.4 × 105 ms−1 and α = 0.174 (24).

Figure 2E plots the measured ξ as a function of EF. The shaded region shows the coherence length estimated within the Thomas-Fermi (TF) approximation Embedded Image, where the spin/valley degeneracy g = 4, α is the fine structure constant, and vF the average Fermi velocity of the kx, ky, and kz directions. Both α and vF are expected to be universal material parameters; however, experimental measurements have not yet tightly constrained them, with α predicted to be between 0.069 and 0.174 and vF to be between 1.4 × 105 and 2.4 × 105 ms−1 (2, 24). Given this uncertainty, Fig. 2E plots the theoretical upper bound as a solid line (using vF = 2.4 × 105 ms−1 and α = 0.069) and the lower bound as a dashed line (using vF = 1.4 × 105 ms−1 and α = 0.174). The agreement is good within uncertainty but closer to the lower bound, implying a larger α (that is, a more strongly interacting system), consistent with previous studies (24, 25). Assuming α = 0.174, we can calculate an impurity density and mobility for regions A, B, and C using the TF approximation and random-phase approximation (RPA) (see section S6 for the full calculations and results). For region C, we infer an nimp of 3.6 × 1018 cm−3 from RPA, in good agreement with the doping measured by the Hall effect (4.35 × 1018 cm−3), demonstrating that the dopants are charged impurities. However, the experimental mobility (3540 cm2 V−1 s−1) is significantly lower than expected (19,000 cm2 V−1 s−1) for charged impurity scattering alone, indicating that other sources of disorder (for example, point defects and grain boundaries) are important and suggests that ultrahigh-mobility thin-film Na3Bi could be achieved, provided that these other sources of disorder can be eliminated.

We now turn to the resonance feature (labeled D) observed in the scanning tunneling spectra in Fig. 1E. Figure 3A plots the peak differential conductance of the resonant feature D, with the red markers indicating the position of defect sites observed in the topography of region A (Fig. 1B). The high degree of correlation indicates that the defect sites are the origin of the resonance feature (see fig. S6 for further analysis). Figure 3B plots dI/dV point spectra taken at a defect site and up to 7 nm away. It is clear the resonance is never completely suppressed, demonstrating that the resonance is not atomically localized at defect sites and has an extended state. These defects correspond to Na surface vacancies (inset of Fig. 1A) in position Na(2) of the crystal structure of Na3Bi shown in Fig. 3C. To better understand the origin of the resonance feature, density functional theory (DFT) calculations including spin-orbit coupling were performed (see Materials and Methods and section S7 for details). Figure 3D compares an experimental STS for Na3Bi on Si(111), where ED is shifted to 0 eV (green curve) with the calculated DOS D(E) for the following Na3Bi structures with an Na(2) vacancy: bulk Na3Bi with one Na(2) vacancy per 8 unit cells (black curve), an Na3Bi slab (2 × 2 × 2 unit cells) containing an Na(2) vacancy in the interior of the slab (blue curve), and an Na3Bi slab (2 × 2 × 3 unit cells) containing an Na(2) vacancy at the surface (red curve). The associated band structures for each of the structures are shown in Fig. 3E. In all cases, the Na(2) vacancy accepts one electron. In addition, we find that Na(2) vacancies at the surface produce a peak in D(E) at an energy close to the experiment (red curve), due to the formation of a “Mexican hat”–shaped valence band edge (see Fig. 3E). This structure and associated D(E) peak are not present for Na(2) vacancies in the bulk or in the interior of confined slabs (see the additional discussion in the Supplementary Materials). Because the enhanced D(E) at the surface results from a band structure effect, we expect that it is not atomically localized at defect sites but rather varies on a length scale set by the Fermi wavelength, in agreement with observations. The extended nature of the resonant state is confirmed from DFT calculations by projecting the DOS on different surface atoms and also on an increasing number of layers to reveal that the resonance spreads across the surface layer and also somewhat penetrates the bulk (calculations are shown in section S7 and figs. S9 and S10). The Na(2) surface vacancy concentration thus provides a unique knob to tune the surface DOS of this Dirac material, which could be used, for example, to tune electron-electron interactions or the coupling to magnetic impurities.

Fig. 3 Determining the bound state defect resonance.

(A) Map of the dI/dV magnitude at the defect resonance energy, where defects are shown as red circles. (B) dI/dV point spectra (taken on region B) on/off the defect site, at locations corresponding to the defect site (black), then 1 nm (red), 3 nm (purple), 5 nm (green), and 7 nm (brown) away from the defect. (C) Crystal structure of Na3Bi, with the surface-terminated Na labeled Na(2) (gold), with the remaining Na atoms in blue with the Na bonded to Bi in the hexagonal lattice labeled Na(1) and the Bi atoms in gray. (D) Comparison between DFT calculations of the DOS for Na3Bi with an Na(2) vacancy at the surface of a 2 × 2 × 3 cell (red curve), an Na(2) vacancy inside a 2 × 2 × 2 cell (blue curve), and a bulk Na(2) vacancy (black curve) and the experimental STS curve for a 20-nm Na3Bi film (green curve). Energy scales of all spectra have been corrected so that 0 eV reflects the Dirac point. A vertical offset has been applied for clarity. (E) Accompanying electronic band structures for bulk Na3Bi with an Na(2) vacancy (left), an Na(2) vacancy inside a 2 × 2 × 2 cell (middle), and an Na(2) vacancy at the surface of a 2 × 2 × 3 cell (right).

Qualitatively different behavior was observed for occasional defect clusters. Figure 4A, where the STM topography of a 60 nm × 60 nm region of Na3Bi, shows a large concentration of singular-defect sites, as well as one defect that appears in topography to consist of a multivacancy cluster. The spectroscopic signature of these defects is very different from the quasi-bound resonant state of the singular defects (see Fig. 3), as shown in the fixed-bias dI/dV maps (Fig. 3, B to D) taken at −196, −216, and −236 meV, respectively. Because of the long-ranged electrostatic interaction of the defect with the scanning tip, a ring-like structure emerges that increases in spatial extent as the tip-sample potential becomes more negative. dI/dV maps were also performed at fixed bias while varying the tunneling current and also showed this ring-like structure that varied in width with changing tunneling current, excluding the possibility of a static LDOS feature. This phenomenon has also been observed for defects in graphene (19), dopants in semiconductors (20, 21), and topological insulator systems (22) as a tip-induced “ionization charging ring,” where the sudden increase in charge of a bound defect state at a particular tip potential affects the screening cloud in the substrate at the tip position. In this case, the defect state is truly localized, in contrast to the quasi-bound state observed for single Na vacancies, and offers an opportunity to study a quantum dot connected to a TDS reservoir, an area of significant theoretical interest (2628).

Fig. 4 Tip-induced ionization rings around large defects.

(A) STM topography (V = −250 mV and I = 250pA) on a 60 nm × 60 nm region of Na3Bi. Fixed-bias dI/dV maps taken at (B) −196 mV, (C) −216 mV, and (D) −236 mV over the same region as (A) showing a ring-like feature centered around a large vacancy site highlighted by the dashed circle in (A).


We have demonstrated, using STM and STS, the existence of charge puddles in the TDS, Na3Bi. The ultralow Dirac point energy fluctuations, which occur over length scales of approximately 5 to 15 nm, are of the order of 4 to 6 meV = 40 to 70 K, well below room temperature and comparable to the highest-quality graphene on h-BN. The ultralow potential fluctuations in this 3D Dirac system will enable the exploration of novel physics associated with the Dirac point (711). In addition, we observed defect-associated quasi-bound and bound states in a 3D TDS, which opens the possibility to tune electron-electron interactions or tune the coupling to magnetic impurities by varying the sodium surface vacancy concentration.


Growth and measurements

The 20-nm Na3Bi thin films were grown in a UHV (10−10 torr) MBE chamber and then transferred immediately after the growth to the interconnected CreaTec low-temperature scanning tunneling microscope operating in UHV (10−11 torr) for STM/STS measurements at 5 K. For Na3Bi film growth, effusion cells were used to simultaneously evaporate elemental Bi (99.999%; Alfa Aesar) in an overflux of Na (99.95%; Sigma-Aldrich) with a Bi/Na flux ratio of not less than 1:10, calibrated by quartz microbalance. The Bi rate used was ~0.03 Å/s, and the Na rate was ~0.7 Å/s. The pressure during growth was less than 3 × 10−9 torr.

Growth on Si(111). To prepare an atomically flat substrate, a Si(111) wafer was flash-annealed to achieve 7 × 7 surface reconstruction, confirmed using STM. During the growth, the substrate temperature was 330°C for successful crystallization. At the end of the growth, the sample was left at 330°C for 10 min in an Na overflux to improve the film quality before cooling to room temperature.

Growth on sapphire. The sapphire substrate was annealed in atmosphere at 1350°C and then in pure oxygen atmosphere at 1050°C to achieve an atomically flat surface. Ti/Au (5/50 nm) contacts were deposited on the corners of the substrate and wire-bonded to a contact busbar on the sample plate to allow for in situ transport measurements in a 1-T perpendicular magnetic field at 5 K. The sapphire was annealed in UHV at 400°C for 1 hour to remove atmospheric species. Na3Bi films were then grown using a two-step growth method, as reported previously by Hellerstedt et al. (24) and Edmonds et al. (25). The Na3Bi film used in this study had a final growth temperature of 330°C.

A PtIr STM tip was prepared and calibrated using an Au(111) single crystal and the Shockley surface state at ~−0.5V and flat LDOS near the Fermi level before all measurements. STM differential conductance (dI/dV) was measured using a 5-mV ac excitation voltage (673 Hz) that was added to the tunneling bias. Differential conductance measurements were made under open feedback conditions with the tip in a fixed position above the surface. Data were prepared and analyzed using MATLAB and WSxM software (29).

DFT methods

First-principles calculations based on DFT were used to obtain electronic bands and DOS of bulk, bilayer, and trilayer Na3Bi, with and without Na(2) vacancies. The first-principles approach is based on Kohn-Sham DFT (30), as implemented in the Quantum ESPRESSO code (31). The exchange correlation energy was described by the generalized gradient approximation using the Perdew-Burke-Ernzerhof functional (32). Interactions between valence and core electrons were described by Troullier-Martins pseudopotentials (33). The Kohn-Sham orbitals were expanded in a plane-wave basis with a cutoff energy of 50 rydberg (Ry), and for the charge density, a cutoff of 200 Ry was used. The Brillouin zone was sampled using Γ-centered grids, following the scheme proposed by Monkhorst and Pack (34). For the convergence of charge density, an 8 × 8 × 4 grid was used for the bulk, whereas a 6 × 6 × 1 grid was used for the bilayer and trilayer. For the DOS, a finer grid of 32 × 32 × 16 (26 × 26 × 1) was used for the bulk (bilayer and trilayer). The calculation of the spin-orbit splitting was performed using noncolinear calculations with fully relativistic pseudopotentials. For the bilayer and trilayer models, we used periodic boundary conditions along the three dimensions, with vacuum regions of 10 Å between adjacent images in the direction perpendicular to the layers. Convergence tests with greater vacuum spacing guaranteed that this size was enough to avoid spurious interactions between neighboring images.


Supplementary material for this article is available at

section S1. Determining the Dirac point position from dI/dV spectra

section S2. Demonstrating the spatial and energy correspondence between STS measurements of the Na3Bi Dirac point and defect quasi-bound state

section S3. Demonstrating that charge puddling is correlated with Na(2) vacancies

section S4. Correlation of resonance state in dI/dV spectra with lattice defects

section S5. Calculation of puddle coherence length from autocorrelation analysis

section S6. Theory discussion on correlation length, impurity density, and mobility

section S7. DFT calculations of Na vacancies in the Na3Bi lattice

fig. S1. Determining the Dirac point position from the dI/dV spectra.

fig. S2. STM topography of region A (as in Fig. 1C) showing the 1 × 1 Na3Bi (001), Na(2)-terminated surface, with vacancy point defects and unidentified impurities visible.

fig. S3. Frequency histogram of the measured Dirac point and Na(2) vacancy quasi-bound state (STS peak, ~30 mV below the Dirac point) features extracted from STS of region A.

fig. S4. Radially averaged correlation profiles for spatial profiles of key STS features in region A.

fig. S5. STM topography and charge puddling map of p-type Na3Bi on region A.

fig. S6. Spatial dependence of defect resonance.

fig. S7. Electronic structure of pristine and defective Na3Bi.

fig. S8. Total DOS of pristine and defective Na3Bi.

fig. S9. Comparison of calculated total DOS and projected DOS on surface for Na3Bi with Na(2) vacancy.

fig. S10. Comparison of calculated total DOS and projected DOS on individual atomic layers for Na3Bi with Na(2) vacancy.

table S1. Charged impurity density for EF » Erms calculated using Thomas– Fermi and RPA for both region A, B, and C.

table S2. Mobility for EF » Erms calculated using Thomas– Fermi and RPA for both region A, B, and C.

References (3539)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: J.L.C., M.T.E., J.H., and M.S.F. acknowledge support from M.S.F.’s Australian Research Council (ARC) Laureate Fellowship (FL120100038). M.T.E. acknowledges support from ARC Discovery Early Career Researcher Award fellowship DE160101157. M.T.E., J.L.C., and M.S.F. acknowledge support from CE170100039. S.A., I.Y., and J.N.B.R. acknowledge support from the National Research Foundation of Singapore’s fellowship program (NRF-NRFF2012-01). L.C.G. acknowledges A. H. Castro Neto and support from the National Research Foundation–Prime Minister’s Office Singapore under its Medium-Sized Centre funding. The first-principles calculations were carried out on the CA2DM (Centre for Advanced 2D Materials) high-performance computing facilities. Author contributions: M.T.E., J.H., and M.S.F. devised the experiments. M.T.E., J.L.C., and J.H. performed the MBE growth and STM/STS measurements. I.Y. and S.A. assisted with the theoretical understanding and interpretation of the charge puddling. L.C.G. and J.N.B.R. performed the DFT calculations. J.L.C., M.T.E., and M.S.F. analyzed the data and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article