Observation of surface superstructure induced by systematic vacancies in the topological Dirac semimetal Cd
3As
2Christopher J. Butler,1,*Yi Tseng,1Cheng-Rong Hsing,2Yu-Mi Wu,1Raman Sankar,3,4Mei-Fang Wang,1Ching-Ming Wei,2,†
Fang-Cheng Chou,4,5,6and Minn-Tsong Lin1,2,7,‡
1Department of Physics, National Taiwan University, Taipei 10617, Taiwan
2Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan
3Institute of Physics, Academia Sinica, Taipei 11529, Taiwan
4Center for Condensed Matter Sciences, National Taiwan University, Taipei 10617, Taiwan
5National Synchrotron Radiation Research Center, Hsinchu 30076, Taiwan
6Taiwan Consortium of Emergent Crystalline Materials (TCECM), Ministry of Science and Technology, Taipei 10622, Taiwan
7Research Center for Applied Sciences, Academia Sinica, Taipei 11529, Taiwan
(Received 12 October 2016; revised manuscript received 28 December 2016; published 28 February 2017) The Dirac semimetal phase found in Cd3As2is protected by a C4rotational symmetry derived from a corkscrew arrangement of systematic Cd vacancies in its complicated crystal structure. It is therefore surprising that no microscopic observation, direct or indirect, of these systematic vacancies has so far been described. To this end, we revisit the cleaved (112) surface of Cd3As2using a combined approach of scanning tunneling microscopy and ab initio calculations. We determine the exact position of the (112) plane at which Cd3As2naturally cleaves, and describe in detail a structural periodicity found at the reconstructed surface, consistent with that expected to arise from the systematic Cd vacancies. This reconciles the current state of microscopic surface observations with those of crystallographic and theoretical models, and demonstrates that this vacancy superstructure, central to the preservation of the Dirac semimetal phase, survives the cleavage process and retains order at the surface.
DOI:10.1103/PhysRevB.95.081410
Dirac and Weyl semimetals host low-energy excitations described by the three-dimensional (3D) Dirac/Weyl Hamilto- nian [1]. As opposed to the two-dimensional (2D) linear dis- persion found in graphene or at topological insulator surfaces [2,3], they host analogous “cones” of three-dimensional linear dispersion, emanating from an even number of paired Weyl nodes in the bulk Brillouin zone (BZ). Where time-reversal or inversion symmetries are broken, pairs of isolated Weyl nodes of nonzero chiral charge can exist, giving rise to a Weyl semimetal phase with associated exotic physical properties [4–7]. Otherwise, Weyl node pairs are forced to coexist as chiral-charge-neutral Dirac nodes, but protected from gap opening perturbations by additional crystal symmetries. This Dirac semimetal (DSM) phase was first proposed in BiO2and A3Bi (A= Na, K, Rb) [1,8], and was realized in Na3Bi [9]
and then in the more stable Cd3As2[10–12].
Numerous transport observations have revealed unusual properties associated with the DSM phase in Cd3As2[13–19], including those involving surfaces, such as chirality transfer between Weyl nodes through Fermi arc surface states [20], and even possible unconventional superconductivity induced at quantum point contact interfaces [21,22].
Dirac nodes are protected by crystal symmetries such as Cn rotational, or nonsymmorphic glide reflection symme- tries. In the case of Cd3As2, they are protected by a C4 rotational symmetry around the [001] axis [23], associated with a peculiar arrangement of systematic Cd vacancies in [001]-oriented corkscrews, in an otherwise antifluoritelike structure. This arrangement has been described in detail as
part of a thorough bulk structural characterization reported by Ali et al. [24]. An apparent discrepancy in the position of the Dirac points (DPs) within the Brillouin zone revolves around these systematic vacancies. The DPs have been located along the (001) direction in angle-resolved photoemission spectroscopy (ARPES) investigations reported by Wang et al.
[23], among others [12,24,25], but along the (112) direction in work reported by Liu et al. [10]. This discrepancy may be attributable to a possible randomization of the Cd vacancy distribution near the cleaved surface [25]. However, it was recently demonstrated that if these Cd vacancies were to be significantly disordered or filled, the symmetry-protected band crossing points, and therefore the DSM phase, would be destroyed [16,26]. Even with recent scanning tunneling microscopy (STM) observations [24,25], it remains that no direct or indirect microscopic observation of these systematic (or random) Cd vacancies has been reported in Cd3As2, despite their central role in the formation of the DSM phase, and the position of the Dirac nodes within the BZ. Thus a gap potentially exists between the prevailing crystallographic understanding of Cd3As2 (with its theoretical implications) and existing microscopic results.
In this Rapid Communication we revisit the most easily attainable surface of Cd3As2, the cleaved (112) surface, using a combined approach of atomically resolved STM measure- ments and density functional theory (DFT) calculations. We construct an atomistic model of the (112) surface expected according to the bulk structure including systematic Cd vacancies, and show that it is predictive of our microscopic observations. This brings into agreement the crystallographic description and microscopic observations of the systematic Cd vacancy arrangement underpinning the DSM phase in Cd3As2. Figure1(a)shows the conventional tetragonal Cd3As2unit cell reported recently [27], with a portion of the atoms cut
FIG. 1. Bulk crystal and electronic structure of Cd3As2. (a) A view of the conventional unit cell, with a cutaway plane showing a (112)-oriented surface. The pink diamond indicates a (112) plane.
(b) The projection of the structure along the (112) axis, showing the pseudohexagonal nature of the surface lattice. (c) A photograph of a typical Cd3As2 platelet as used in this work. (For scale, 1 small square= 1 mm2). The primary facet has a hexagonal shape (as indicated by the yellow dashed line) and is identified as a (112) surface. (d) A bulk band structure diagram, and (e) a zoom-in view of the region along -Z near the Fermi level, showing the location of the DPs.
away to illustrate the (112) surface. The pink diamond shows one of the (112) planar cuts through the unit cell, parallel to which are the natural cleavage planes. The projection of this structure looking along the normal to the cleavage planes is displayed in Fig.1(b), showing a pseudohexagonal lattice. A photograph of a single-crystal Cd3As2platelet as used here is shown in Fig.1(c). Platelets used in this work were grown as described previously [27], and have large shiny primary facets with hexagonal shapes, corresponding to the (112) plane of the crystal. Figure1(d)shows a bulk electronic band structure diagram calculated for this structure, with high-symmetry axes taken from the structure’s smaller primitive cell [28].
Figure1(e)gives a zoom-in view showing the DP located near the Fermi level between the and Z points of the BZ. All band structure and charge density calculations presented in this work were performed in the framework of DFT, using the projector- augmented wave (PAW) method [29,30] as implemented in theVASPpackage [31–33], with the Perdew-Burke-Ernzerhof exchange-correlation functional [34] and fully accounting for spin-orbit coupling. Structural optimization was performed using a 2× 2 × 1 Monkhorst-Pack k-point mesh [35], and a kinetic energy cutoff of 300 eV. The force on each atom was less than 0.01 eV ˚A−1.
In order to perform STM measurements on the (112) surface, bulk single crystals of Cd3As2 were cleaved at room temperature in UHV conditions, before quickly being transferred to an Omicron LT-STM held at 4.5 K. An elec- trochemically etched tungsten tip was used, and the obtained STM images were processed using theWSXMsoftware package
FIG. 2. Interpretation of topographic step heights. (a) A large- scale STM topograph (V = 2.0 V, I = 0.1 nA), with the topographic height profile, taken along the white dashed line, shown in (b).
Small “substeps” can be seen in some areas, marked by the white arrows in (a) and visible in the height profile. (c) A side view of the structural model [viewing along the [111] direction, and with the (112) surface normal pointing upwards] shows the division of the structure into trilayerlike slabs stacked in alternating orientations.
Planar projections of the two resulting inequivalent surface lattices, labeled I and II, are shown in (d) and (e), respectively.
[36]. Figure 2(a)shows the typical step-terrace morphology observed in large-scale STM images taken at the (112) surface of Cd3As2. Previously, the precise position of the cleavage planes through the crystal had yet to be determined, though it has been suggested by Jeon et al., based on STM surface lattice images, that the cleaved surface terminates with an As atomic layer [25].
We determine the precise position of the cleavage planes along the (112) normal by inspecting the step heights of the step-terrace morphology along the dashed line in Fig.2(a).
The height profile along this line is shown in Fig.2(b). The predominant step height, around 7.2 ˚A, corresponds to twice the As-As interlayer spacing parallel to the (112) normal.
Figure2(c)shows a side view along the [111] direction of the crystal stacking. It can be seen that the basic stacking units are formed by layers of As atoms sandwiched on each side by Cd layers. Each of the observed topographic steps corresponds to a pair of these “trilayer sandwich” stacking units. The pairing is due to the inequivalence of the planes separating the stacking units, which we label as cleavage planes I and II in Fig.2(c). [And since for each type of cleavage plane, the
resulting pair of cleaved surfaces on either side are structurally equivalent with one another, these labels carry over from the cleavage planes to the resulting surface: Cleavage through plane I (or II) always results in a type I (or II) surface.] As shown in Figs.2(d)and2(e), the main structural difference between the two surfaces is that three quarters of the systematic Cd vacancies of the Cd3As2 structure are located in the Cd atomic layer which terminates surface II. This results in linear channels of Cd vacancies arranged along the [111] direction of the conventional unit cell, indicated by red circles in Fig.2(c).
In fact, cleavage at both planes I and II is observed to occur.
In Fig.2(a), small substeps (marked with white arrows) can be seen protruding from the major steps, and have a height corresponding to half the major step height, or one trilayer sandwich, as is shown in the topographic line profile shown in Fig.2(b). This indicates that both cleavage planes I and II are observable, but that one is strongly energetically preferred over the other. Here the question remains: Which cleaved surface, I or II, is the predominant one? This can be answered by detailed inspection of the observed surface atomic lattice, to which we will return below.
At this stage, we conclude that the step heights of the observed step-terrace morphology are consistent with recent descriptions of the Cd3As2 structure, and that the cleavage planes are located between adjacent Cd atomic layers, which would suggest a Cd-terminated surface rather than the As surface previously suggested [25]. A Cd-terminated surface is in fact more in line with previously reported x-ray photo- electron spectroscopy measurements, which have shown Cd 4d doublets dominating in intensity over As 3d doublets [10,27]. (The Cd 4d and As 3d levels have approximately the same photoionization cross sections in the low photon energy range [37].)
Figure 3 shows high-resolution STM topographs of the cleaved (112) surface. The expected pseudohexagonal atomic surface lattice and the previously reported 2× 2 reconstruction are observed. The most outstanding feature is a periodic, uniaxial zigzag pattern, highlighted with dashed white lines in topographs shown in Figs. 3(a) and 3(b), which are acquired in the same position but using V = −0.5 and −0.1 V, respectively. The overall periodicity of the surface superlattice can be characterized with a rectangular unit cell (solid white lines), with lattice parameters a= 17.9 ˚A, b = 15.6 ˚A.
Large-scale topographs show the ubiquity of this pattern over the entire measured surface [28]. These features were not observed in STM measurements reported by Jeon et al. [25], and possible reasons for this discrepancy are discussed below.
Figure3(c)shows a scheme for reorienting the Cd3As2unit cell so that one of its faces corresponds to the rectangular surface lattice detailed in Figs. 3(a) and 3(b). The relation between the reoriented and original lattice vectors is expressed by the unitary transformation [a,b,c]= [a,b,c]U, where
U =
⎡
⎣ 1 0.5 0.5
−1 0.5 0.5
0 −0.5 0.5
⎤
⎦.
This alternative bulk unit cell now forms the basis of our slab model for ab initio calculations, and unlike the unit cell shown in Fig.1(a), allows the construction of a suitable slab
FIG. 3. Interpretation of atomically resolved topography. Zoom- in STM topographs taken using V = −0.5 and −0.1 V, with I = 1 nA, are shown in (a) and (b), respectively. Both images show a zigzaglike superstructure, which can be characterized by a 2D unit cell shown as a white rectangle in each image. A scheme for setting up a reoriented unit cell suitable for slab-model calculations is shown in (c). The rectangular upper face of the unit cell represents the 2D unit cell identified in (a) and (b). Side views of the resulting structurally optimized slab models for surfaces I and II are shown in (d) and (e), respectively. (f) A zoom-in view of the measured STM topography, and the simulated STM images for surfaces I and II, shown in (g) and (h), respectively, are given for comparison.
with a choice of (112) surface terminations facing the vacuum.
Slab models terminated by surfaces I and II were constructed with a thickness of eight “trilayer sandwich” stacking units [28]. Structural optimization resulted in a small deviation from bulklike terminations (i.e., terminations without any reconstruction) as shown in Figs. 3(d) and 3(e). A simple bulklike termination for either surface I or II would result in either of the Cd surface lattices as shown in Figs.2(d)and2(e).
Instead, in both cases surface relaxation results in a subsidence of the surface Cd lattice down to around the same level as the second (As) layer. Despite this intermingling of the surface Cd and As lattices, the number density of observed atoms in STM images is not increased compared to that expected for a bulklike termination. This may be because the density of states (DOS) is strongly dominated either by the Cd 5s states (conduction band) or by the As 4p states (valence band) [23], so that one lattice, either the Cd or As lattice (see below), provides the overwhelming contribution to the tunnel current and topography image. These effects greatly complicate the
interpretation of STM images, since it is no longer guaranteed that surfaces I and II can be identified simply by their patterns of surface lattice vacancies.
Having obtained optimized structures for the two candidate surfaces, we calculate the partial charge density distribution in order to obtain simulated STM topographs. In order for the simulated and experimental images to be comparable, these partial charge density plots should be obtained by integrating the local density of states (LDOS) between the Fermi level and the energy corresponding to the tip voltage. However, this is not as straightforward as usual: As has been consistently reported in ARPES investigations [11,12,27] and scanning tunneling spectroscopy measurements [25,27], the DPs are found∼200 meV below the Fermi level [28]. This is at odds with the expectation, based on the charge neutrality condition [23], that the Dirac nodes reside at or near the Fermi level, as is also predicted using ab initio calculations [11,23,24].
This n-type shift has been speculated to result from donor impurities, with a prime candidate being As vacancies [38]. For the STM image shown in Fig.3(b), taken at−0.1 eV, the bias is at (EDP+ 100) meV. The LDOS integral to obtain the closest corresponding partial charge density map should then be made between EF = (EDP+ 200) meV, and (EDP+ 100) meV.
Simulated STM maps for candidate surfaces I and II were obtained from calculated partial charge density distributions by considering the Tersoff-Hamann model [39] and using a tip-sample distance of∼4 ˚A [28]. These maps are shown in Figs. 3(g)and 3(h), respectively, alongside a zoom-in STM topography image in Fig. 3(f). In each simulated map, the 2D projected atom positions for the uppermost layer of both Cd and As are illustrated. Surprisingly, though the cleavage plane was located between two adjacent Cd layers, and the topmost Cd and As layers reside at about the same height, the dominant contribution to the observable surface lattice comes from As, possibly due to a dominance of As 4p orbitals in the DOS.
Both of the simulated STM images successfully reproduce the expected rectangular periodicity of the observed super- structure, and an additional periodicity which may correspond to the 2× 2 reconstruction observed by Jeon et al. Despite the challenges detailed above, we attempt to match some key features of the simulated STM images with features in the measured STM topograph. Most obviously, the simulated STM for surface I exhibits a pronounced zigzag pattern, which matches the zigzag appearing in the measured topographs shown in Figs. 3(a)and 3(b). The center line of the zigzag passes through a line of protrusions, several of which are marked with yellow dots in Fig.3(g). We find corresponding atoms in the measured topograph [also marked with yellow dots in Fig.3(f)], which also lie along the center line of the zigzag. In both the measured image and the image simulated
for surface I, this line of atoms possesses the maximum topographic height. No obviously comparable zigzag, nor any corresponding set of high protrusions, can be found in the simulated STM map for surface II [Fig.3(h)]. Based on these arguments, we identify the surface observed in measurements as surface I.
Our initial results, showing that the (112) cleavage plane is located between two adjacent Cd atomic layers, but that the terminating layer is a relaxed mixture of Cd and As atoms, clarifies upon the interpretation given by Jeon et al. It was previously thought that the observation of an apparently complete atomic surface lattice (without obvious systematic vacancies) indicated the full hcp planar lattice expected for an As-terminated (112) surface [25]. While the corrugations observed in STM topographs likely correspond to the uppermost As layer, they do not directly reflect the actual surface structure.
A more concerning difference is the absence in previous work of the surface superstructure described here. This might be explained in a number of ways: (i) a difference in samples, with the samples used by Jeon et al. possessing a structural phase different from that described in current models [24,27], and (ii) a difference in the STM tip condition. In any given STM experiment, the tip apex has its own unknown orbital character, which determines the details of the tunneling transition matrix between the tip and sample, introducing unknown selection effects into topography maps. It is not impossible that different tip conditions may explain the absence in previous reports of the observed superstructure.
In our main result, we provide a microscopic description of a periodic surface superstructure caused by a pattern of systematic Cd vacancies, the symmetry of which underpins the DSM phase in Cd3As2. This observation reconciles the status of real-space microscopic observations with the established structure determined using x-ray diffraction (XRD), and with observations of the broad range of physical phenomena thought to accompany it. It also shows that the systematic ordering of vacancies at the surface survives the cleavage process. This surface characterization should help in the un- derstanding of surface and interface phenomena which might be sensitively termination dependent, such as surface Fermi arcs and other possible trivial or nontrivial surface states, or possible unconventional superconductivity at heterostructure interfaces.
The authors thank W.-L. Lee for stimulating discussions.
This work was supported in part by the Ministry of Science and Technology and National Science Council of Taiwan through the following Grants No. MOST 104-2119-M-002-029 and No. NSC 101-2119-M-002-007.
C.J.B. and Y.T. contributed equally to this work.
[1] S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe,Phys. Rev. Lett. 108, 140405 (2012).
[2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim,Rev. Mod. Phys. 81,109(2009).
[3] M. Z. Hasan and C. L. Kane,Rev. Mod. Phys. 82,3045(2010).
[4] T. Ojanen,Phys. Rev. B 87,245112(2013).
[5] P. Hosur and X. Qi,C. R. Phys. 14,857(2013).
[6] A. C. Potter, I. Kimchi, and A. Vishwanath,Nat. Commun. 5, 5161(2014).
[7] Y. Baum, E. Berg, S. A. Parameswaran, and A. Stern,Phys. Rev.
X 5,041046(2015).
[8] Z. Wang, Y. Sun, X. Q. Chen, C. Franchini, G. Xu, H. Weng, X.
Dai, and Z. Fang,Phys. Rev. B 85,195320(2012).
[9] Z. K. Liu, B. B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D.
Prabhakaran, S. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen,Science 343,864(2014).
[10] Z. K. Liu, J. Jiang, B. B. Zhou, Z. J. Wang, Y. Zhang, H. M.
Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen,Nat. Mater. 13,3990(2014).
[11] M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, A. Bansil, F.-C.
Chou, and M. Z. Hasan,Nat. Commun. 5,3786(2014).
[12] S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy, B.
B¨uchner, and R. J. Cava,Phys. Rev. Lett. 113,027603(2014).
[13] L. P. He, X. C. Hong, J. K. Dong, J. Pan, Z. Zhang, J. Zhang, and S. Y. Li,Phys. Rev. Lett. 113,246402(2014).
[14] J. Cao, S. Liang, C. Zhang, Y. Liu, J. Huang, Z. Jin, Z.-G. Chen, Z. Wang, Q. Wang, J. Zhao, S. Li, X. Dai, J. Zou, Z. Xia, L. Li, and F. Xiu,Nat. Commun. 6,7779(2015).
[15] H. Pan, K. Zhang, Z. Wei, B. Zhao, J. Wang, M. Gao, L. Pi, M.
Han, F. Song, X. Wang, B. Wang, and R. Zhang,Appl. Phys.
Lett. 108,183103(2016).
[16] S.-T. Guo, R. Sankar, Y.-Y. Chien, T.-R. Chang, H.-T. Jeng, G.-Y.
Guo, F.-C. Chou, and W.-L. Lee,Sci. Rep. 6,27487(2016).
[17] T. Liang, Q. Gibson, M. N. Ali, M. Liu, R. J. Cava, and N. P.
Ong,Nat. Mater. 14,4143(2014).
[18] J. Feng, Y. Pang, D. Wu, Z. Wang, H. Weng, J. Li, X. Dai, Z.
Fang, Y. Shi, and L. Lu,Phys. Rev. B 92,081306(2015).
[19] C. Li, L.-X. Wang, H. Liu, J. Wang, Z. Liao, and D. Yu,Nat.
Commun. 6,10137(2015).
[20] P. J. W. Moll, N. L. Nair, T. Helm, A. C. Potter, I. Kimchi, A. Vishwanath, and J. G. Analytis,Nature (London) 535,266 (2016).
[21] L. Aggarwal, A. Gaurav, G. S. Thakur, Z. Haque, A. K. Ganguli, and G. Sheet,Nat. Mater. 15,4455(2015).
[22] H. Wang, H. Wang, H. Liu, H. Lu, W. Yang, S. Jia, X.-J. Liu, X. C. Xie, J. Wei, and J. Wang,Nat. Mater. 15,38(2016).
[23] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang,Phys. Rev. B 88,125427(2013).
[24] M. N. Ali, Q. Gibson, S. Jeon, B. B. Zhou, A. Yazdani, and R. J. Cava,Inorg. Chem. 53,4062(2014).
[25] S. Jeon, B. B. Zhou, A. Gyenis, B. E. Feldman, I. Kimchi, A. C.
Potter, Q. D. Gibson, R. J. Cava, A. Vishwanath, and A. Yazdani, Nat. Mater. 13,851(2014).
[26] X. Yuan, P. Cheng, L. Zhang, C. Zhang, J. Wang, Y. Liu, Q. Sun, P. Zhou, D. W. Zhang, Z. Hu, X. Wan, H. Yan, Z. Li, and F. Xiu, arXiv:1601.02316.
[27] R. Sankar, M. Neupane, S.-Y. Xu, C. J. Butler, I. Zeljkovic, I. Panneer Muthuselvam, F.-T. Huang, S.-T. Guo, S. K. Karna, M.-W. Chu, W.-L. Lee, M.-T. Lin, R. Jayavel, V. Madhavan, M. Z. Hasan, and F.-C. Chou,Sci. Rep. 5,12966(2015).
[28] See Supplemental Material athttp://link.aps.org/supplemental/
10.1103/PhysRevB.95.081410for an illustration of the primitive Cd3As2 unit cell and the corresponding BZ, a typical large- scale topography map, and a scanning tunneling spectroscopy curve obtained at the (112) surface, structures in CIF format for the reoriented bulk and both slab models, examples of simulated STM for various tip-sample distances, and a shortPYTHONcode used to obtain simulated STM images fromPARCHGproduced byVASP.
[29] P. E. Bl¨ochl,Phys. Rev. B 50,17953(1994).
[30] G. Kresse and D. Joubert,Phys. Rev. B 59,1758(1999).
[31] G. Kresse and J. Hafner,Phys. Rev. B 48,13115(1993).
[32] G. Kresse and J. Furthm¨uller, Comput. Mater. Sci. 6, 15 (1996).
[33] G. Kresse and J. Furthm¨uller,Phys. Rev. B 54,11169(1996).
[34] J. P. Perdew, K. Burke, and M. Ernzerhof,Phys. Rev. Lett. 77, 3865(1996).
[35] J. D. Pack and H. J. Monkhorst,Phys. Rev. B 16,1748(1977).
[36] I. Horcas, R. Fern´andez, J. M. G´omez-Rodr´ıguez, J. Colchero, J. G´omez-Herrero, and A. M. Baro, Rev. Sci. Instrum. 78, 013705(2007).
[37] J.-J. Yeh, Atomic Calculation of Photoionization Cross Sections and Asymmetry Parameters (Gordon and Breach, Langhorne, PA, 1993).
[38] D. P. Spitzer, G. A. Castellion, and G. Haacke,J. Appl. Phys.
37,3795(1966).
[39] J. Tersoff and D. R. Hamann,Phys. Rev. Lett. 50,1998(1983).
