Spatial Charge Inhomogeneity and Defect States in Topological Dirac Semimetal Thin Films
Mark T. Edmonds, James L. Collins, Jack Hellerstedt, Indra Yudhistira, Lídia C. Gomes, João N. B. Rodrigues, Shaffique Adam, Michael S. Fuhrer
11 Spatial Charge Inhomogeneity and Defect States in Topological Dirac Semimetal Thin Films
Mark T. Edmonds , James L. Collins , Jack Hellerstedt , Indra Yudhistira , Lídia C. Gomes , João N. B. Rodrigues , Shaffique Adam , Michael S. Fuhrer School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia 2.
Monash Centre for Atomically Thin Materials, Monash University, Clayton VIC 3800 Australia 3.
Department of Physics and Centre for Advanced 2D Materials, National University of Singapore, 117551, Singapore 4.
Yale-NUS College, 6 College Avenue East, 138614, Singapore † These authors contributed equally to this work * Corresponding author: [email protected] and [email protected]
The close approach of the Fermi energy E F of a Dirac semimetal to the Dirac point E D uncovers new physics such as velocity renormalization, and the Dirac plasma at | E F - E D | < k B T , where k B T is the thermal energy. In graphene, substrate disorder drives fluctuations in E F . Three-dimensional topological Dirac semimetals (TDS) obviate the substrate, and should show reduced E F fluctuations due to better metallic screening and higher dielectric constants. Here we map the potential fluctuations in TDS Na Bi using a scanning tunneling microscope. The rms potential fluctuations are significantly smaller than room temperature (Δ E F ,rms = 4-6 meV = 40-70 K) and comparable to the highest quality graphene on h-BN; far smaller than graphene on SiO , or the Dirac surface state of a topological insulator. Surface Na vacancies produce a novel resonance close to the Dirac point with surprisingly large spatial extent and provides a unique way to tune the surface density of states in a TDS thin-film material.
Three-dimensional topological Dirac semimetals (TDS) such as Na Bi and Cd As express the pseudorelativistic physics of two-dimensional Dirac material graphene, but extended to three dimensions. TDS can yield ultra-high mobilities as well as new physics such as the chiral anomaly. As the Fermi energy E F measured relative to the Dirac point energy E D approaches zero, the carrier density tends to zero and metallic screening disappears. Poorly screened disorder induces spatial fluctuations in the local E D and hence fluctuations Δ E F , corresponding to local electron and hole “puddles”. The carriers in puddles in turn screen the disorder, with Δ E F determined self-consistently by the disorder and the screening properties of the Dirac materials. Puddles have been visualized in graphene using scanning single electron transistor microscopy and scanning tunnelling spectroscopy, where the fluctuations are largely governed by the underlying substrate. Here we use scanning tunneling microscopy and spectroscopy (STM/STS) to probe Δ E F fluctuations in 20 nm Na Bi films on both semiconducting [Si(111)] and insulating [ -Al O (0001)] substrates. Figure 1(a) shows a large area (400 nm x 380 nm) topographic STM image of a thin 20nm film of Na Bi on Si(111), with several atomically flat terraces >100 nm in size. The inset to Figure 1(a) reveals the (1 x 1) Na-terminated atomic lattice (a = 5.45 Å) with an individual Na vacancy. Figures 1(b) (45 nm x 45 nm taken immediately after growth) and 1(c) (30nm x 30nm taken seven days after growth) show the topography of two atomically flat regions away from step edges or the screw dislocations seen in Fig. 1(a). Figure 1(d) shows an atomically flat region (30 nm x 30 nm) of Na Bi grown on sapphire ( -Al O (0001). Whilst atomically flat regions of Na Bi up to 100 nm x 100 nm can be obtained on Si(111), sparse defect cluster sites (few per 100 nm x 100 nm) give rise to tip-induced ionization ring features that will be discussed further below. This necessitated focusing on smaller areas free of ionization rings in order to unambiguously determine the variation in Dirac point. Figure 1(e) shows area-averaged scanning tunneling spectra (STS) of the Na Bi film. STS measures the differential conductance d I /d V as a function of sample bias V , proportional to the local density of states (DOS) at energy eV relative to E F . STS was performed at points on a grid encompassing the entire regions shown in Figures 1(b & c), along with two other regions not shown. Figure 1(e) shows the averaged spectra of these four areas. Two key features are prominent in all spectra: a distinct minimum in the DOS corresponding to the Dirac-point energy ( E D ), and a resonant feature ~35 meV below E D , labelled D. As-grown, the Dirac point is located ~20 meV above the Fermi level indicating p -type doping consistent with angle-resolved photoelectron spectroscopy (ARPES) measurements. 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 Na Bi due to the adsorption of atomic species present in UHV. This adsorption results in the formation of intermittent impurity clusters on the surface, as such all topography and spectroscopic measurements were deliberately performed away from such sites. The resonance feature D is unambiguously tied to E D , and not to the E F , as the relative energy shift of D with respect to E D 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, E D , found by tracking the position of the minimum differential conductance in STS (alternatively, the shift in the defect resonance D in each spectrum gives similar results). Figures 2(a) and 2(b) show the spatial variation of E D for regions ‘A’ and ‘B’ corresponding to Figures 1(b) and 1(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. This modulation in E D represents the puddling of charge density at the surface. The variation in E D is found to be correlated with the individual Na vacancy sites which possess a more positive Dirac point energy, identifying the Na vacancies as acceptors and a significant source of charge disorder (see below and Supplementary Information). In addition, Figure 2(c) shows the local E D of 20 nm Na Bi on -Al O (0001) (labelled Region C) which possesses a larger n -type doping. The upper, middle and lower panels of Fig. 2(d) show histograms of E D relative to E F for the scans in Figs. 2(a)-(c) respectively. From spatial autocorrelation analysis of the puddle maps, we determine spatial coherence lengths ξ of 13.4 ± 5.2 nm, 9.3 ± 2.4 nm and 5.1 ± 1.9 nm for Regions A, B and C respectively (see SI for calculation details). The observed histograms of E D have mean values 19.7 ± 1.7 meV (Region A), -15.4 ± 1.3 meV (Region B) and -47.2 ± 0.6 meV (Region C), and standard deviations Δ E F of 5.6 meV, 4.2 meV and 3.5 meV respectively. These values are comparable to the 5.4 meV observed for graphene on h-BN. However, undersampling results in an underestimate in the magnitude of Δ E F by a factor ( L /ξ) /[( L /ξ) -1] (also likely in Ref. 8), hence the corrected Δ E F values are 6.1 meV and 4.6 meV for Regions A and B respectively, whilst Region C remains essentially unchanged due to the small coherence length. Figure 2(e) plots the measured ξ as a function of E F . The shaded region shows the coherence length estimated within the Thomas-Fermi approximation 𝜉 = √ 𝜋2𝛼𝑔 ( ℏ𝑣 𝐹 𝐸 𝐹 ) where the spin/valley degeneracy g = is the fine structure constant and v F the Fermi velocity. The agreement is good within uncertainty but closer to the lower bound, implying a larger (i.e. more strongly interacting system) consistent with previous studies. Assuming = 0.174, we can calculate an impurity density and mobility for Regions A, B and C using the Thomas-Fermi (TF) and Random Phase (RPA) approximations (see Section S6 in Supplementary Information for full calculation and results). For Region C we infer n imp of 3.6x10 cm -3 from RPA, in good agreement with the doping measured by Hall effect (4.35 x 10 cm -3 ), demonstrating that the dopants are charged impurities. However, the experimental mobility (3540 cm V -1 s -1 ) is significantly lower than expected (19,000 cm V -1 s -1 ) for charged impurity scattering alone, indicating that other sources of disorder (e.g. point defects, grain boundaries) are important and suggests that ultra-high mobility thin-film Na Bi could be achieved provided these other sources of disorder can be eliminated. We now turn to the resonance feature (labelled D) observed in the STS spectra in Fig. 1(d). Figure 3(a) 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. 1(b)]. The high degree of correlation indicates that the defect sites are the origin of the resonance feature (see Fig. S3 in SI for further analysis). These defects correspond to Na surface vacancies (inset of Fig. 1(a)) in position Na(2) of the crystal structure of Na Bi shown in Fig. 3(b). To better understand the origin of the resonance feature, DFT calculations including spin-orbit coupling were performed (see Methods). Figure 3(c) compares an experimental STS for Na Bi on Si(111) (green curve) with the calculated density of states D ( E ) for the following Na Bi structures with an Na(2) vacancy: Bulk Na Bi with one Na(2) vacancy per 8 unit cells (black curve), a Na Bi slab (2x2x2 unit cells) containing an Na(2) vacancy in the interior of the slab (blue curve), and a Na Bi slab (2x2x3 unit cells) containing an Na(2) vacancy at the surface (red curve). The associated bandstructures for each of the structures are shown in Fig. 3(d). In all cases the Na(2) vacancy accepts one electron. Additionally we find that indeed 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. 3(d)). 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 additional discussion in Supplementary Information). Since the enhanced D ( E ) at the surface results from a bandstructure effect, we expect it is not localised at defect sites but rather varies on a length scale set by the Fermi wavelength, in agreement with observations. The Na(2) surface vacancy concentration thus provides a unique knob to tune the surface density of states of this Dirac material, which could be used e.g. to tune electron-electron interactions, or tune the coupling to magnetic impurities. Qualitatively different behavior was observed for occasional defect clusters. Figure 4(a) where the STM topography of a 60 nm x 60 nm region of Na Bi shows a large concentration of singular defects sites, as well as one defect which appears in topography to consist of a multi-vacancy cluster. The spectroscopic signature of such defects is very different from the quasi-bound resonant state of the singular defects [see Figure 3], as shown in the fixed-bias d I /d V maps (b)-(d) taken at -196 meV, -216 meV and -236 meV respectively. Due to 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. This phenomenon has also been observed for defects in graphene , and dopants in semiconductor systems 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. We have demonstrated using scanning tunneling microscopy and spectroscopy the existence of charge puddles in the topological Dirac semimetal, Na Bi. The ultra-low Dirac-point energy fluctuations, which occur over length scales of approximately 5-15 nm, are of the order of 4-6 meV = 40-70 K well below room temperature and comparable to the highest quality graphene on h-BN. The ultra-low potential fluctuations in this 3D Dirac system will enable the exploration of novel physics associated with the Dirac point.
In addition, we observed for the first time defect-associated quasi-bound and bound states in a 3D TDS, which open the possibility to tune electron-electron interactions, or tune the coupling to magnetic impurities by varying the sodium surface vacancy concentration.
Experimental Methods:
The 20 nm Na Bi thin film was grown in a ultra-high vacuum (UHV) (10 -10
Torr) molecular beam epitaxy (MBE) chamber and then transferred immediately after the growth to the interconnected Createc LT-STM operating in UHV (10 -11
Torr) for STM/STS measurements at 5 K. For Na Bi 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 not less than 1:10, calibrated by quartz microbalance. The Bi rate used was ~0.03 Å/s, and Na was ~0.7 Å/s. The pressure during growth was less than 3x10 -9 Torr.
Growth on Si(111) -
To prepare an atomically flat substrate, a Si(111) wafer was flash annealed in order to achieve 7 x 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 a 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 pure oxygen atmosphere at 1050 ⁰ C to achieve an atomically flat surface. Ti/Au (5/50nm) contacts were deposited on the corners of the substrate, and wirebonded to a contact busbar on the sample plate to allow for in-situ transport measurements in a 1T perpendicular magnetic field at 5K. The sapphire was annealed in UHV to 400 ⁰ C for 1 hour to remove atmospheric species. Na Bi films were then grown using a two-step growth method as reported previously in [Ref. 19]. The Na Bi film used in this study had a final growth temperature was 330 ⁰ C. A PtIr STM tip was prepared and calibrated using an Au(111) single crystal and the Shockley surface state before all measurements. STM differential conductance (d I /d V ) was measured using a 5 mV AC excitation voltage (673 Hz) that was added to the tunnelling bias. Differential conductance measurements were made under open feedback conditions with the tip in a fixed position above the surface. Data was prepared and analysed using MatLab and WSxM software. Density functional theory methods:
First-principles calculations based on density-functional theory are used to obtain electronic bands and density of states of bulk, bilayer and trilayer Na Bi, with and without Na(2) vacancies. The first-principles approach is based on Kohn-Sham density functional theory (KS-DFT), as implemented in the Quantum ESPRESSO code. The exchange correlation energy is described by the generalized gradient approximation (GGA) using the PBE functional. Interactions between valence and core electrons are described by Troullier-Martins pseudopotentials. The Kohn-Sham orbitals were expanded in a plane-wave basis with a cutoff energy of 50 Ry, and for the charge density, a cutoff of 200 Ry was used. The Brillouin-zone (BZ) was sampled using -centered grids, following the scheme proposed by Monkhorst-Pack. For convergence of charge density, an 8x8x4 grid was used for bulk while a 6x6x1 grid was used for bilayer and trilayer. For the density of states, a finer grid of 32x32x16 (26x26x1) is employed for bulk (bilayer and trilayer). The calculation of the spin-orbit splitting was performed using non-colinear 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, guarantee that this size is enough to avoid spurious interactions between neighbouring images.
References Elias, D. C., Gorbachev, R. V., Mayorov, A. S., Morozov, S. V,. Zhukov, A. A., Blake, P., Ponomarenko, L. A., Grigorieva, I. V., Novoselov, K. S., Guinea, F., and Geim, A. K. Dirac cones reshaped by interaction effects in suspended graphene, Nature Phys. 7, 701 (2011) Luican, A., Li, G., and Andrei, E. Y. Quantized Landau level spectrum and its density dependence in graphene, Phys. Rev. B 83, 041405(R) (2011). Chae, J., Jung, S., Young, A. F., Dean, C. R., Wang, L., Gao, Y., Watanabe, K., Taniguchi, T., Hone, J., Shepard, K. L., Kim, P., Zhitenev, N. B., and Stroscio, J. A. Renormalization of the graphene dispersion velocity determined from scanning tunneling spectroscopy, Phys. Rev. Lett. 109, 116802 (2012). Crossno, J., Shi, J. K., Wang, K., Liu, X., Harzheim, A., Lucas, A., Sachdev, S., Kim, P., Tanihuchi, T., Watanabe, K., Ohki, T. A., Fong, K. C. Observation of the Dirac fluid and the breakdown of the Wiedemann-Franz law in graphene, Science 343,1058 (2016) Bandurin, D. A., Torre, I., Krishna Kumar, R., Ben Shalom, M., Tomadin, A., Principi, A., Auton, G. H., Khestanova, E., Novoselov, K. S., Grigorieva, I. V., Ponomarenko, L. A., Geim, A. K., Polini, M. Negative local resistance caused by viscous electron backflow in graphene, Science 351, 1055 (2016) Wang, Z., Sun, Y., Chen, X.-Q., Franchini, C., Xu, G., Weng, H., Dai, X., and Fang, Z. Dirac semimetal and topological phase transitions in A Bi (A=Na, K, Rb), Phys. Rev. B 85, 195320 (2012) Liu, Z. K., Zhou, B., Zhang, Y., Wang, Z. J., Weng, H. M., Prabhakaran, D., Mo, S.-K., Shen, Z. X., Fang, Z., Dai, X., Hussain, Z., Chen, Y. L. Discovery of a three-dimensional topological Dirac semimetal, Na Bi, Science 343, 864 (2014) Xue, J., Sanchez-Yamagishi, J., Bulmarsh, D., Jacquod, P., Deshpande, A., Watanabe, K., Taniguchi, T., Jarillo-Herroro, P., and LeRoy, B. Scanning tunnelling microscopy and spectroscopy of ultra-flat graphene on hexagonal boron nitride, Nature Mater. 10, 282 (2011) Zhang, Y., Brar, V. W., Girit, C., Zettl, A., Crommie, M. F. Origin of spatial charge inhomogeneity in graphene, Nature Phys. 5, 722 (2009) Martin, J., Akerman, N., Ulbricht, G., Lohmann, T., Smet, J. H., von Klitzing, K., and Yacoby, A. Observation of electron-hole puddles in graphene using a scanning single-electron transistor, Nature Phys. 4, 144 (2008) Beidenkopf, H., Roushan, P., Seo, J., Gorman, L., Drozdov, I., Hor, Y. S., Cava, R. J., and Yazdani, A. Spatial fluctuations of helical Dirac fermions on the surface of topological insulators, Nature Phys. 7, 939 (2011) Liu, Z. K., Jiang, J., Zhou, B., Wang, Z. J., Zhang, Y., Weng, H. M., Prabhakaran, D., Mo, S.-K., Peng, H., Dudin, P., Kim, T., Hoesch, M., Fang, Z., Dai, X., Shen, Z. X., Feng, D. L., Hussain, Z., and Chen, Y. L. A stable three-dimensional topological Dirac semimetal Cd As , Nature Materials 13, 677 (2014) Borisenko, S., Gibson, Q., Evtushinky, D., Zabolotnyy, V., Büchner, B., and Cava, R. J. Experimental realization of a three-dimensional Dirac semimetal, Phys. Rev. Lett. 113, 027603 (2014) Liang, T., Gibson, Q., Ali, M. N., Liu, M., Cava, R. J., and Ong, N. P. Ultrahigh mobility and giant magnetoresistance in the Dirac semimetal Cd As , Nature Materials 14, 280 (2015) Xiong, J., Kushwaha, S. K., Lian, T., Krizan, J. W., Hirschberger, M., Wang, W., Cava, R., Ong, N. P. Evidence for the chiral anomaly in the Dirac semimetal Na Bi, Science 350, 413 (2015) Adam, S., Hwang, E. W., Galitski, V. W., and Das Sarma, S. A self-consistent theory for graphene transport, Proceedings of the National Academy of Science 104, 18392 (2007) Ramakrishnan, N., Milletarì, M., and Adam, S. Transport and magnetotransport in three-dimensional Weyl semimetals, Phys. Rev. B
92, 245120 (2015) Samaddar, S., Yudhistira, I., Adam, S., Courtois, H., Winkelmann, C. B. Charge puddles in graphene near the Dirac point, Phys. Rev. Lett. 116, 126804 (2016) Hellerstedt, J., Edmonds, M. T., Ramakrishnan, N., Liu, C., Weber, B., Tadich, A., O’Donnell, K. M., Adam, S., and Fuhrer, M. S. Electronic properties of high-quality epitaxial topological Dirac semimetal thin films, Nano Letters 16, 3210 (2016) Edmonds, M. T., Hellerstedt, J., O’Donnell, K. M., Tadich, A., and Fuhrer, M. S. Molecular doping the topological Dirac semimetal Na Bi across the charge neutrality point with F4-TCNQ, ACS Applied Materials and Interfaces 8, 16412 (2016) Brar, V. W., Decker, R., Solowan, H.-M., Wang, Y., Maserati, L., Chan, K. T., Lee, H., Girit, C. O., Zettl, A., Louie, S. G., Cohen, M. L., and Crommie, M. F. Gate-controlled ionization and screening of cobalt adatoms on a graphene surface, Nature Phys. 7, 43 (2011) Marczinowski, F., Wiebe, J., Meier, F., Hashimoto, K., and Wiesendanger, R. Effect of charge manipulation on scanning tunneling spectra of single Mn acceptors in InAs, Phys. Rev. B 77, 115318 (2008). Zhang, Y., Liu, Z., Zhou, B., Kim, Y., Hussain, Z., Shen, Z.X., Chen, Y., and Mo, S.-K. Molecular beam epitaxial growth of a three-dimensional topological Dirac semimetal Na Bi, Applied Physics Letters 105, 031901 (2014) Horcas, I., Fernandez, R., Gomez-Rodriguez, J. M., Colchero, J., Gomez-Herrero, J., Baro, A. M. WSXM: A software for scanning probe microscopy and a tool for nanotechnology, Rev. Sci. Instrum. 78, 013705 (2007) Kohn, W., Sham, L. J. Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965). Giannozzi, P., et al., Quantum Espresso: a modular and open-source software project for quantum simulations of materials, J. Phys.-Cond. Matter. 21, 395502 (2009). Perdew, J. P., Burke, K., Ernzerhof, M. Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996). Troullier, N., Martins, J. L. Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B. 43, 1993 (1991). Monkhorst, H. J., and Pack, J. D. Special points for Brillouin-zone integrations, Phys. Rev. B. 13, 5188 (1976).
Acknowledgements
J. L. C., M. T. E., J. H and M. S. F. are supported by M. S. F.’s ARC Laureate Fellowship (FL120100038). M. T. E. is supported by ARC DECRA fellowship DE160101157. S. A., I. Y. and J. N.
B. R. are supported under the National Research Foundation of Singapore’s fellowship program (NRF-NRFF2012-01). L. C. G. acknowledges A. H. Castro Neto and is supported by 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 high-performance computing facilities.
Author Contributions
M. T. E. and J. L. C. contributed equally to this work. M. T. E., J. H., and M. S. F. devised the experiments. M. T. E. and J. L. C. 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 composed the manuscript.
Additional Information
The authors declare no competing financial interests. Supplementary information can be found.
Figure 1. Large area morphology and spectroscopy of 20 nm Na Bi on Si(111). (A) Large area (400nm x 380 nm) topographic STM image (bias voltage V = -3 V and tunnel current I = 50 pA). Inset: Atomic resolution STM image with lattice constant 5.45 Å (taken on separate 20 nm Na Bi film) showing an individual Na vacancy at the surface. (B) STM topography ( V = 300 mV and I =
250 pA) on a 45 nm x 45 nm region of Na Bi. (C) STM topography (V = -
550 mV and I = 200 pA) on a 30 nm x 30 nm region of Na Bi. (D) STM topography ( V = -50 mV and I =
100 pA) on a 30 nm x 30 nm region of Na Bi on sapphire ( -Al O (0001) . (E) Area-averaged STS 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, whilst the blue spectra corresponding to the region in (C) and green spectra (topography not shown) were taken one week after growth. Feature E D reflects the Dirac point, which corresponds to the minimum in the LDOS, and D represents the resonance feature (discussed in text). Figure 2. Charge puddling profiles of p -type and n -type Na Bi. (A) Dirac point energy map of 45 nm x 45 nm (90 pixels x 90 pixels) region of Na Bi ( V = -250 mV and I = 250 pA), corresponding to Region A represented in Figure 1(B). Scale bar is 15nm. (B) Dirac-point energy map of 30 nm x 30 nm (60 pixels x 60 pixels) region of Na Bi corresponding to Region B of Figure 1(C) ( V = -150 mV and I = 200 pA). Scale bar is 10 nm. (C) Dirac-point energy map of 30 nm x 30 nm (60 pixels x 60 pixels) region of Na Bi grown on -Al O (0001) (labelled Region C) ( V = -150 mV and I = 200 pA). Scale bar is 10 nm. (D) Upper, middle and lower panels representing histograms of the Dirac point energy maps in (A)-(C) respectively. The histograms are colour-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 v F = ms -1 and = 0.069, whilst the lower bound (dashed line) uses v F = ms -1 and = 0.174[19]. Figure 3. Determining the bound state defect resonance. (A) Map of the d I /d V magnitude at the defect resonance energy, where defects are shown as red circles. (B) Crystal structure of Na Bi, with the surface terminated Na labelled Na(2) (gold), with the remaining Na atoms in blue with the Na bonded to Bi in the hexagonal lattice labelled Na(1) and the Bi atoms in grey. (C) Comparison between DFT calculations of the DOS for Na Bi with a Na(2) vacancy at the surface of a 2x2x3 cell (red curve), a Na(2) vacancy inside a 2x2x2 cell (blue) and bulk Na(2) vacancy (black curve) and the experimental STS curve for a 20 nm Na Bi 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. (d) Accompanying electronic bandstructures for bulk Na Bi with a Na(2) vacancy (left panel), Na(2) vacancy inside a 2x2x2 cell (middle panel) and Na(2) vacancy at the surface of a 2x2x3 cell (right panel). Inset: Brillouin zone of bulk Na Bi and the projected (001) surface. Figure 4. Tip-induced ionization rings around large defects. (A) STM topography (V = -250 mV and I = 250pA) on a 60 nm x 60 nm region of Na Bi. Fixed bias d I /d V maps taken at (B) -196 mV, (C) -216 mV and (D) -236 mV over the same region as (A) showing a ring like feature centred around a large vacancy site highlighted by the dashed circle in (A). Supplementary Information Spatial Charge Inhomogeneity and Defect States in Topological Dirac Semimetal Thin Films
Mark T. Edmonds , James L. Collins , Jack Hellerstedt , Indra Yudhistira , Lídia C. Gomes , João N. B. Rodrigues , Shaffique Adam , Michael S. Fuhrer School of Physics and Astronomy, Monash University, Clayton VIC 3800, Australia 2.
Monash Centre for Atomically Thin Materials, Monash University, Clayton VIC 3800 Australia 3.
Department of Physics and Centre for Advanced 2D Materials, National University of Singapore, 117551, Singapore 4.
Yale-NUS College, 6 College Avenue East, 138614, Singapore † These authors contributed equally to this work * Corresponding author: [email protected] and [email protected]
TABLE OF CONTENTS: S1. Determining the Dirac point position from dI/dV spectra S2. Demonstrating that charge puddling is tied to Na(2) vacancies S3. Correlation between resonance feature in dI/dV spectra with lattice defects S4. Calculation of puddle coherence length from autocorrelation analysis S5. Theory discussion on correlation length, impurity density and mobility S6. DFT calculations of Na vacancies in the Na Bi lattice S1. Determining the Dirac point position from the dI/dV spectra
For a Dirac semimetal such as graphene, or in the present case of 3D Dirac semimetal Na Bi, the density of states (DOS) should be at a minimum at the Dirac-point Energy (E D ). This manifests in the scanning tunneling spectra as a minima in the differential conductance (dI/dV) as a function of sample bias, which is proportional to the local density of states. The simplest approach to determining the minima in the dI/dV (i.e. Dirac point) is fitting the spectra to a low order polynomial. However, in the case of Na Bi the presence of a defect-associated bound energy resonance ~33meV below E D means standard polynomial fitting does not always yield accurate determination of the Dirac point as shown in Fig. S1. To counter this problem a non-parametric local regression procedure in MATLAB was employed to produce fits of the dI/dV spectra to then identify and track critical points of the DOS (i.e. E D and also the defect resonance) with reduced systematic error. This procedure then allowed us to construct puddling maps from either (1) the change in the minima directly, or (2) from the change in the defect resonance position D (which also reflect changes in E D ) and then correcting for the constant energy separation of E D – D = 35 meV. Figure S1. Determining the Dirac point position from the dI/dV spectra.
A single dI/dV spectra (black curve) taken at an arbitrary point on an atomically flat region of Na Bi on Region B (Fig. 1(c) main manuscript). This spectra is fit to a low-order polynomial (red curve) and a non-parametric local regression (blue curve). The large discrepancy between the minima position in the dI/dV spectra and the minima obtained from the low-order polynomial highlights the need for a detailed fitting procedure to ensure accurate determination of Dirac point position. S2. Demonstrating that charge puddling is tied to Na(2) vacancies
In order to demonstrate that the defects attributed to Na vacancies of the lattice in Na Bi are correlated with the fluctuations in the energy potential we compare the STM topography and charge puddle map of Region A in Fig. S2(a) and (b) respectively. Black circles are used to highlight the defect sites in the charge puddles map in Fig. S2(b). In Fig. S2(c) histograms representing the entire Region A (orange stairs) and only the black marker locations (blue bars) are plotted. A shift to higher energy of ~4 meV is observed. This is consistent with our DFT calculations which indicate that the Na(2) vacancy accepts one electron (see Fig. S5). We conclude that Na(2) vacancies acts as negatively charged impurities, increasing the local hole doping. In order to determine the expected doping from the sodium vacancies we compare the number of defects observed in the topography with the predicted impurity density calculated using the random phase approximation. In Fig. S2(a) ~133 sodium vacancies are located within the 45 nm x 45 nm scan region, equating to 6.6x10 sodium vacancies per cm , giving 3.3x10 cm -3 for 20 nm of Na Bi (a similar defect concentration is found for Region B). This is similar to the predicted impurity density of 2.9-4.5 x10 cm -3 calculated for Region A in Table S1. This demonstrates that sodium vacancies are a major source of the impurity density and disorder in pristine Na Bi films grown on Si(111). However, it should be noted for the n -type doped Na Bi there are clearly other sources of doping than sodium vacancies which are acceptors. At present we cannot determine the exact source of the large n -type doping, or whether those dopants lie at the substrate interface or distributed throughout the bulk of the Na Bi film.
Figure S2 . STM topography and charge puddling map of p -type Na Bi on Region A. (A) STM topography corresponding to Region A (Fig. 1(B)) of the main manuscript. (B) Charge puddling map corresponding to Fig. 2(A) of the main manuscript, with black circles representing the defect sites. (C) Two histograms of subsets of the Dirac-point energy in Region A, which have been normalized for clarity. The orange stairs represent the Dirac-point energy over all of Region A (with the mean of 20 meV, marked as a dashed orange line) and the blue bars represent the Dirac-point energy measured only at the black marker locations (with the mean of 24 meV marked as a blue dashed line). S3. Correlation between resonance state in dI/dV spectra with lattice defects
To further demonstrate the relationship between the magnitude of the bound resonant state and lattice defects we examine the defect resonance and lattice defects of Region A in Fig. S3. In Fig. S3(a) the dI/dV resonance intensity is plotted for Region A, along with markers highlighting the defect positions. In Fig. S3(b) we plot the dI/dV resonance histograms of all of Region A (orange stairs) and only the defect sites marked by green circles (blue bars). We omit the defects highlighted in red due to those regions showing faint signatures of charging ring behaviour, which affects the resonance intensity. The resonance differential conductance magnitude over all of Region A has a mean value at 4.9pA and standard deviation of 0.90pA, whilst at only the green defect locations there is a mean value of 5.4pA with standard deviation 0.66pA. This shift to higher dI/dV intensity at the defect sites highlights the direct correlation between the resonance state and defect location.
Figure S3 . Spatial dependence of defect resonance. (A) Map of the dI/dV magnitude at the defect resonance energy corresponding to Fig. 3(A) in the main manuscript. Markers indicate the location of lattice defects, with the defects assigned red and green marker, based on whether their location shows (red) and does not show (green) faint signatures of charging-ring behaviour. (B) Two histograms of subsets of the resonance magnitude in Region A, which have been normalized for clarity. The orange stairs represent the resonance differential conductance magnitude over all of Region A (with the mean of 4.9pA marked as a dashed orange line) and the blue bars represent the resonance differential conductance measured only at the green marker locations (with the mean of 5.4pA marked as a blue dashed line). S4. Calculation of correlation length, impurity density and mobility
The coefficients of an autocorrelation matrix
𝐶(𝑥 , 𝑥 ) for charge puddle fluctuations in E D , is predicted to exhibit exponential spatial dependence in two-dimensions. With the assumption that there should be no necessary anisotropy of the measurement, it is convenient to use a radially-averaged C ( R ). The average coherence length of C ( R ) can then be taken as the puddle length scale 𝜉 ; this is the distance at which C ( R ) = 𝐶(0)( − 1) , where 𝑒 is Euler’s number. As the experimental length scale L is not much larger than the observed fluctuation length scale, it is suspected that undersampling is responsible for producing anti-correlation coefficients in C ( R ) at length scales greater than 𝐿 /π . We attempt to account for this in two ways: (1) We perform an exponential regression over the range L, which does not account for undersampling-related uncertainty in C ( R ). (2) We perform linear regression of log 𝐶(𝑅) , weighted by an estimate for error 𝛿𝐶(𝑅) . On a length scale R, there are 𝐿 − 𝑅 pairs of measurements which are used for C(R), however the number of statistically different measurements is of order − 𝑅 )/𝜉 . Taking 𝛿𝐶(𝑅) then as the inverse standard deviation of measurements, 𝛿𝐶(𝑅) ≅ (1 + 𝐿 −𝑅 𝜉 ) −2 . By performing least-squares fitting of log 𝐶(𝑅) over the range of positive 𝐶(𝑅) , we achieve an estimate of 𝜉 from goodness of fit analysis after 𝛿𝐶(𝑅) for a given 𝜉 propagates. By taking a lower-bound ( 𝜒 = 1 best fit) and upper-bound (1 𝜎 overlap of fit to model with error), the mean of these two values is then our estimate for 𝜉 . Additionally, as charge puddles are expected to follow a Gaussian energy distribution, we place an estimate of 𝛿|𝐸 𝐹 | as the standard error, where mean 𝜎 𝑚𝑒𝑎𝑛 = 𝜎/√𝑁 , and further assume that the number of measurements is 𝑁 = 𝐿 /𝜉 . S5. Theory discussion on correlation length, impurity density and mobility
To theoretically compute the correlation functions in Na Bi, we make two assumptions. First, the impurity potential comes from a random three-dimensional (3D) distribution of charged impurities with density . Second, we assume that it is possible to find a global screening function that can describe the effects of these local screening variations. With these assumptions, we can write the autocorrelation function of the screened Coulomb potential as (1) (2) where is impurity density, is the spatial distance between two points, and . Here we have used long-range Coulomb potential (3) where is the electron charge, is the effective dielectric function of the material, and is the momentum transfer between incoming and outgoing momenta of the scattered electron. The dielectric function is (4) where is polarizability function, is the density of states of Na Bi, and . Introducing effective fine structure constant , we get the inverse screening length . Here is degeneracy factor, is Fermi velocity and is effective Fermi wavevector. imp n ),( eff nq riq enq qVqdnrC ),( )(=)( qrqrnqq qdqen )(sin),(8= qrqrqq qdqen )(sin8= rgIekn ,24= ,2 )(2sin)(2~=),( effeff2eff22 20 rk rkkzdrzI imp n r )/(2= eff kq ,4=)( qeqV e q )()(1=),( eff qqVnq )(~)()(1= eff qnqV ,)(1= qnq )( q F vkgn )()/()(~ eff nqq )/(= F ve effeffs )/(~2=)( kqgnq g F v eff k Thomas-Fermi polarizability is given by or , hence . In this case can be evaluated analytically as (5) Hence the autocorrelation function can be calculated as (6) We would like to remark that the autocorrelation function has an exponential spatial profile, in contrast with autocorrelation function in monolayer graphene, that has a gaussian profile.
We can then define Thomas-Fermi correlation length as nothing but Thomas-Fermi effective screening length Since in our experiment the Fermi energy is larger than the fluctuation in energy , effective wavevector can very well be approximated by Fermi wavevector , leading to (7) Here, we have used effective low energy dispersion relation of Na Bi, in the last line. We plotted the upper bound and lower bound of the Thomas-Fermi correlation length plotted in figure 2(d) of the main paper using the formula above with and , respectively. We have used to take spin and valley degeneracy into account. The ratio of the random phase approximation (RPA) polarization function and the density of states is given by the sum of two components: a vacuum part and a finite density part )(=)( effTF nq TF q effTF /2= kgq ),( rzI TF rk rkzdrzI TF effeff222 20 .4= eff2 zrk ez rkqIeknrC ,24=)( effTFTF22effimpTF .)(2= )eff(TF22effTF imp rnq eenq n )(1/=)( effTFeffTF nqn .12= eff kg F E rms E F kk eff F kg TF .||2= FF Evg FFF kvnsgnE )(= eff5 F v eff5 F v g )]/(2=[~ RPA F kqx )(~ V x )(~ M x x xxxxxx (8) Here , where is ultraviolet cutoff. Analytical evaluation of RPA autocorrelation function is not possible. However, the numerics of RPA autocorrelation function is found to also have exponential decay with spatial distance just like Thomas– Fermi autocorrelation function. Hence, we can obtain the RPA correlation length by fitting the numerics of the ratio of RPA autocorrelation function at finite and at to exponential function (9) We have found that RPA correlation length doesn’t differ much from Thomas– Fermi correlation length.
2 Impurity density
The fluctuation in the screened, charged impurity-induced potential is given by (10) (11) Since in our experiment the Fermi energy is larger than the fluctuation in energy , effective wavevector can very well be approximated by Fermi wavevector , leading to (12) where we have used effective low energy dispersion relation of Na Bi, in the second line. Inverting this relation, we obtain charged impurity density (13) .log32=~ xx )/(2= F k r RPA r r .exp0)=( )( RPARPARPA rrC rC rms E = VE rC ),( )()(2= nq qVqdn qq qdqen eff imp2eff2 )(,24= k nvkgJ F .)(2~=),(
222 20 kzdkzJ F E rms E F kk eff FFF k nvkgJE imp222rms )(,24= ,)(,24= imp32 FFF
E nvkgJ FFF kvnsgnE )(= .)(,24 1= F FF v EEkgJn Thomas– Fermi dimensionless polarizability is given by , hence (14) Therefore, charged impurity density for rms » EE F is given by (15) Within RPA, the charged impurity density for rms » EE F is given by (16) (17) where need to be evaluated numerically. Both Thomas– Fermi and RPA upper and lower bound for charged impurity density for region A, B, and C are computed in table S1. Table S1: Charged impurity density for rms » EE F calculated using Thomas– Fermi and RPA for both region A, B, and C. Region A Region B Region C Here, we have used to take into account valley and spin degeneracy and ultraviolet cutoff for RPA.
3 Mobility
The conductivity of Na Bi is given by (18) TF q .4==),(
222 20TF zzdkzJ .)( ||)2(= FF v EEgn )(,24 1= F FF v EEkgJn ,)(2~=),( kzdkzJ k RPA J cm10 n cm10 n g F k ),(3 )(= FFF
EEve where is density of states at Fermi level and is transport scattering time. The Boltzmann transport scattering time for Coulomb impurities within Born approximation is given by (19) (20) where is charged impurity density, is the momentum transfer between incoming plane-waves with wavevector and outgoing wavefunctions with wavevector and . Hence, the conductivity is given by .2112= imp4222 gHnkghe F (21) Therefore, mobility is (22) Substituting equation of for rms » EE F (eq. 13), we get 𝜇 ≈ 𝑒ℏ𝜋 ( 𝑣 𝐹 𝐸 𝑟𝑚𝑠 ) 𝑔𝛼2𝜋 ,𝑘 𝐹 )𝐻(√ 𝑔𝛼2𝜋 ) (23) FFF vkgE )(2cos1)( )()(2 '2=)( '2233imp k k EEqqVdnE FF FFFF kkVdvkn
222 20222imp FFF kgdekvn gHkvn FF .)(2~ )(1=)(
222 2310 F kzdzH imp n |'=| kk q k ' k )/(2= F kq ne = .2121= imp2 gHnkhe F imp n Within Thomas-Fermi approximation, we can obtain analytical expression for (24) Hence, the mobility for rms » EE F is given by (25) RPA mobility for rms » EE F is given by similar expression 𝜇 ≈ 𝑒ℏ𝜋 ( 𝑣 𝐹 𝐸 𝑟𝑚𝑠 ) 𝑅𝑃𝐴 (√ 𝑔𝛼2𝜋 ,𝑘 𝐹 )𝐻 𝑅𝑃𝐴 (√ 𝑔𝛼2𝜋 ) (26) (27) where both J RPA need to be evaluated numerically. Both Thomas– Fermi and RPA upper and lower bound for mobility for region A, B, and C are computed in table S2. TABLE S2: Mobility for rms » EE F calculated using Thomas– Fermi and RPA for both region A, B, and C. Region A Region B Region C Here, we have used to take into account valley and spin degeneracy and ultraviolet cutoff for RPA. We can see that both Thomas– Fermi approximation and RPA yield the same mobility. )( zH
222 2310TF )(1=)( zdzH zz TF2rms gHEvge F ,)(2~ )(1=)( F kzdzH RPA H sVm sVm g F k S6. DFT calculations of Na vacancies in the Na Bi lattice
We have used Density Functional Theory (DFT) calculations to investigate the electronic properties of Na Bi slabs with Na-vacancies at the surface. There are two inequivalent Na atoms in the unit cell: Na(1) atoms [represented by dark blue spheres in Fig. 3(b)]; and Na(2) atoms [represented by the light blue and yellow spheres in Fig. 3(b)]. Total energy calculations indicate that the Na(2) surface vacancies [i.e. absence of the yellow atom in Fig. 3(b)] are energetically more favorable than Na(1) surface vacancies. As can be seen in the leftmost bandstructure of Fig. 3(d), the introduction of a Na(2) vacancy in the bulk Na Bi (2x2x2 supercell) preserves the gapless Dirac semimetal structure of pristine Na Bi. Nevertheless, as can be seen in panels (a) and (c) of Fig. S4, the position of the Dirac point is slightly moved away from the Gamma-point, while the lowest energy bands become flatter in the vicinity of the Dirac point. The removal of a Na(2) atom shifts the Fermi level to the valence band, as expected for the removal of this low electronegativity atom. Integrating the DOS in Fig. S5 for the thin-slabs from the Fermi level to the band gap gives one hole per Na(2) vacancy. Equivalently, we note that the removal of a single Na(2) atom removes 9 valence electrons from the system, but only 8 spin non-degenerate bands below E F . This indicates that the Na(2) vacancy can be seen as effectively accepting one electron. Furthermore, there is no defect-induced flat band or localized state in the bandstructure, nor a resonance in D(E) as shown by the full black curve in Fig. S5.
Confinement of the electrons along the z-direction by considering a finite slab [along the vertical direction in Fig. 3(b)] a bandgap is opened (with its size dependent on the thickness of the Na Bi as expected). This gap is significantly larger in the slab with a Na(2) vacancy at the surface, when compared with the gap of a slab without Na-vacancies – see panels (b) and (e) of Fig. S4 corresponding to bilayers of Na Bi, respectively, without and with Na(2) vacancy at the surface. Nevertheless, as expected, increasing the thickness of the Na Bi slab with a Na(2) vacancy at the surface decreases the magnitude of the gap – compare panels (e) and (f) of Fig. S4 where we show the bandstructure of a Na Bi bilayer slab and of a Na Bi trilayer slab, both with a Na(2) vacancy (2x2 supercell). The slabs with Na(2) vacancies at the surface show pronounced peaks in the
D(E) near the gap edges – see green and red curves in Fig. S5. Such resonances are not present in the
D(E) of either the pristine slab (i.e., slab without vacancy) or the slab with an interior Na(2) vacancy, as shown respectively by the dashed-black and the full-blue curves of Fig. S5. Notably, the trilayer slab with a Na(2) surface vacancy, shows not only a narrower bandgap, but also a more pronounced valence band resonance than the conduction band resonance. This result is consistent with the experiment, where no
D(E) resonance is seen at the conduction band. Furthermore, the ab initio
D(E) peaks appear at approximately the same energy position (~50meV below the band edge for the trilayer) of the resonant feature experimentally observed (~35meV below E D ). As can be seen in Fig. S4(e,f), the peaks in the D(E) result from the formation of a "Mexican hat" shaped valence and conduction band edges in the bandstructure. The fact that by increasing the slab thickness (from bilayer to trilayer) the
D(E) peaks magnitude changes only slightly, while the confinement gap decreases sharply (from 1.39 eV to 0.70 eV), suggests that the peaks do not scale with the size of the confinement gap, but rather originate from the presence of a surface. This is further confirmed from the absence of a resonant peak (as well as absence of “Mexican hat” in the bandstructure) when the Na(2) vacancy is in the middle of the 2x2x2 slab – see Fig. S4(d) and the full-blue line in Fig. S5. From the DFT calculations we can then conclude that the observed resonance in the D(E) in the experiment is a result of the modification of the bandstructure due to the presence of the surface Na(2) defect, coupled with the formation of a confined state.
Figure S4. Electronic structure of pristine and defective Na Bi. (A) Pristine bulk. (B) Pristine bilayer simulated with a 1x1 supercell. (C) Bulk with a Na(2) vacancy simulated with a 2x2x2 supercell. (D) Bilayer with an interior Na(2) vacancy simulated with a 2x2 supercell. (E) Na(2) vacancy at the surface of a bilayer (simulated with a 2x2 supercell). (F) Na(2) vacancy at the surface of a trilayer (simulated with a 2x2 supercell). All electronic dispersions were calculated including SOC effects.
Figure S5. Total density of states (DOS) of pristine and defective Na Bi.
The grey shadow stands for the DOS of the pristine bulk Na Bi, while the dashed-black curve stands for the DOS of the pristine bilayer slab. The black-full curve stands for the DOS of the bulk Na Bi with a Na(2) vacancy (2x2x2 supercell), the blue curve identifies the DOS of a bilayer slab with an interior Na(2) vacancy (2x2 supercell), while the green and red curves stand for the DOS of, respectively, a bilayer and a trilayer slab with a surface Na(2) vacancy (2x2 supercell). References Adam, S., Jung, S., Klimov, Nikolai, N. N., Zhitenev, N. B., Stroscio, J. A., and Stiles, M. D., Mechanism for puddle formation in graphene,
Phys. Rev. B , 84,235421 (2011). Samaddar, S., Yudhistira, I., Adam, S., Courtois, H., and Winkelmann, C. B.,. Charge puddles in graphene near the Dirac point,
Phys. Rev. Lett. , 116,126804, (2016). Hellerstedt, J., Edmonds, M. T., Ramakrishnan, N., Liu, C., Weber, B., Tadich, A., O’Donnell, K. M., Adam, S., and Fuhrer, M. S. Electronic properties of high-quality epitaxial topological Dirac semimetal thin films,
Nano Letters
16, 3210 (2016) LV, M., and Zhang, S.-C.,
Dielectric function, friedel oscillations, and plasmons in Weyl semimetals,
International Journal of Modern Physics B , 27(25), 1350177 (2013). Skinner, B. Coulomb disorder in three-dimensional Dirac systems,
Phys. Rev. B , 90,060202 (2014). Das Sarma, S., Hwang, E. H., and Min, H., Carrier screening, transport, and relaxation in three-dimensional Dirac semimetals,
Phys. Rev. B , 91,035201 (2015). Ramakrishnan, N., Milletari, M., and Adam, S. Transport and magnetotransport in three-dimensional Weyl semimetals,