Experimental realization of a topological p-n junction by intrinsic defect-grading
T. Bathon, S.Achilli, P.Sessi, V.A.Golyashov, K.A.Kokh, O.E.Tereshchenko, M.Bode
EExperimental realization of a topological p–n junctionby intrinsic defect-grading
T. Bathon, S. Achilli, P. Sessi, ∗ V. A. Golyashov, K. A. Kokh,
4, 5, 6
O. E. Tereshchenko,
3, 5, 6 and M. Bode
1, 7 Physikalisches Institut, Experimentelle Physik II,Universit¨at W¨urzburg, Am Hubland, D-97074 W¨urzburg, Germany Fisica, Universit`a Cattolica di Brescia, via dei Musei 41, I-25121 Brescia, Italy A.V. Rzanov Institute of Semiconductor Physics, Siberian Branch,Russian Academy of Sciences, 630090 Novosibirsk, Russia V.S. Sobolev Institute of Geology and Mineralogy, Siberian Branch,Russian Academy of Sciences, 630090 Novosibirsk, Russia Novosibirsk State University, 630090 Novosibirsk, Russia Saint-Petersburg State University, 198504 Saint-Petersburg, Russia Wilhelm Conrad R¨ontgen-Center for Complex Material Systems (RCCM),Universit¨at W¨urzburg, Am Hubland, D-97074 W¨urzburg, Germany (Dated: September 17, 2018)
A junction between an n- and p-type semiconductor results in the creation of a deple-tion region whose properties are at the basis of nowadays electronics. If realized usingtopological insulators as constituent materials, p-n junctions are expected to manifestseveral unconventional effects with great potential for applications. Experimentally,all these fascinating properties remained unexplored so far, mainly because prototyp-ical topological PNJs, which can be easily realized and investigated, were not readilyavailable. Here, we report on the creation of topological PNJs which can be as narrowas few tenths of nm showing a built-in potential of 110meV. These junctions are in-trinsically obtained by a thermodynamic control of the defects distribution across thecrystal. Our results make Bi2Te3 a robust and reliable platform to explore the physicsof topological p-n junction.
The recent discovery of topologically protected surface states in the binary chalchogenides Bi Te [1], Bi Se [2] andSb Te [3] has tremendously revitalized the interest in these narrow gap semiconductors that have been studied fordecades because of their excellent thermoelectric properties [4]. Topological insulators (TIs) are materials insulatingin the bulk but metallic on their surface due to the existence of linearly dispersing gapless states which, contraryto the trivial surface states usually found in metals and semiconductors, are protected by time reversal symmetry[5]. The strong spin-orbit coupling characterizing these materials perpendicularly locks the spin to the momentum,resulting in a chiral spin texture which forbids backscattering [6, 7] and makes spin currents intrinsically related tocharge currents [8].These unconventional properties make TIs an ideal platform to realize exotic states of matter, such as Majoranafermions [9] and magnetic monopoles [10]. Furthermore, they appear to be ideal candidates to realize magneto-electricand spintronics devices with low power consumption [11]. Within this framework, the creation of p–n junctions (PNJ),which are the building blocks of several semiconducting devices, such as diodes, sensors, solar cells, or transistors, intopological materials would represent the first step towards the direct application of this fascinating class of materials[12]. In TIs, the massless character of surface carriers and the helical spin texture inversion at the p–n interface arepredicted to result in several unconventional phenomena, e.g. gapless chiral edge mode [12], Klein tunneling [13],Veselago lenses [14], exciton condensation and charge fractionalization [15].Experimentally, all these fascinating properties remained unexplored so far, mainly because prototypical topologicalPNJs, which can be easily realized and investigated, were not readily available. One issue towards the successfulfabrication of topological PNJs is the presence of defects which are intrinsically incorporated in TI materials duringthe growth process. Thereby, the resulting TI crystals are heavily doped at a relatively high bulk carrier concentrationlevel [16, 17]. Since bulk carriers give rise to leakage currents that make capacitive charging impossible their existencerepresent a severe obstacle towards the creation of topological PNJs by applying a gate bias between opposite surfacesof a TI thin film [15].Recent studies showed, however, that—depending on the growth conditions—large deviations from the nominalcomposition can be obtained [18]. These defects, which can be anti-sites or vacancies, are characterized by different a r X i v : . [ c ond - m a t . m e s - h a ll ] D ec qu i n t up l e l a y e r Te(I)Bi(II)Te(III)Bi(IV)Te(V)p-type n-type
Growth direction
Bi-rich Te-rich
FIG. 1: Lower panel: Photographic image of a Bi Te crystal grown by the modified Bridgman technique. As schematicallyillustrated by different colors a doping gradient leads to a position-dependent transition from p- to n-type conductivity withthe chemical potential lying in the upper and lower part of the Dirac cone, respectively. This transition is accompanied by aninversion of the helical spin texture, as indicated by green arrows in the upper left and right panel. The crystal structure isshown in the hatched box. formation energies. Here, we demonstrate that their inevitable presence can be turned into a positive effect bythermodynamically controlling their distribution across the crystal. In particular, we show how appropriate growthconditions result in the creation of bulk crystals that intrinsically contain topological PNJs at their surface, therebyavoiding the complicated fabrication of heterostructures and the problems related to the creation of interfaces betweendifferent materials [19–21]. Atomic scale scanning tunneling microscopy (STM) and spectroscopy (STS) measurementscombined with ab initio calculations evidence that Te segregation [17, 22, 23] results in well-separated Te- and Bi-richregions that display p- and n-transport character, respectively. Spatially resolved Hall and Seebeck measurementsconfirm the transition from p- to n-like transport. At the p–n interface the oppositely charged defects, i.e. donors oracceptors, compensate, resulting in a substantial drop of conductivity by almost two orders of magnitude. Scanningtunneling spectroscopy reveals a built-in potential of about 110 meV, thereby considerably shifting the Dirac pointin between the p- and n-region. Spectroscopic data obtained within the p–n compensation region indicate that thejunction width amounts to several 10 nm only.Figure 1 schematically presents some essential background information related to our approach. Within the hatchedbox of Fig. 1 the crystal structure of a quintuple layer of the binary chalchogenide Bi Te is shown. It consists ofalternating Te and Bi layers. Adjacent quintuple layers are weakly bound by van der Waals forces thereby offeringa natural cleavage plane, a situation highly advantageous for STM experiments. The bottom panel of Fig. 1 showsa photographic image of a Bi Te crystal grown by the modified Bridgman technique (see Experimental Section fordetails). The crystal has a nominal composition with a Te content of 61 mol%. As recently shown, the conductionof Bi Te crystals changes from p- to n-type at approximately this Te concentration [24]. Due to Te segregation,however, the stoichiometry will not be constant but continuously change during the growth process from Bi-rich toTe-rich, as indicated by different color shadings in the bottom panel of Fig. 1. Correspondingly, we expect Diracpoints which are energetically located below and above the Fermi level, respectively, with opposite rotational senseof the spin polarization. As we move from p- to n-doped surface areas we expect to pass through a transition regionwith equal concentrations of Bi- and Te-induced charge carriers where the carrier concentration becomes minimal,thereby realizing an intrinsic TI.In fact, the differently doped areas can clearly be distinguished by STM and STS. Fig. 2(a) displays an overviewtopographic STM image obtained at the (0001) surface of the Bi-rich crystal region. Four different defects labeled(i)-(iv) are visible on the surface. Since the defects symmetry and extension reflects the perturbation introduced bythe bonding structure within the crystal [25] (see bond geometry sketched in the inset of Fig. 1), a detailed analysis ofthe atomically resolved images reported in the top row of Fig. 2(b) allows to identify their location within the quintuplelayer structure of Bi Te , i.e. Te(I)-Bi(II)-Te(III)-Bi(IV)-Te(V), where Te(I,V) and Te(III) represent two inequivalentTe planes. Following this procedure the defects (i)-(iv) are all located within Te layers, and they can be Te vacancies(V Te ) or antisites (Bi Te ). Due to the lack of chemical sensitivity of STM a definite assignment can only be achieved -200 0 2001234 E - E /meV F VB CB E D -400 -200 0 20012345 E - E /meV F VB CB vvi viiviii nm iiiiii iv nm f iii ivv vi vii viii nm 2 nm 4 nm3 nm 3 nm 4 nm E D d I / d U / a r b . un i t s CB d I / d U / a r b . un i t s VB CB vvi viiviii nm iiiiii iv nm ad c iii ivv vi vii viii nm 2 nm 4 nm3 nm 3 nm 4 nm E D ii nm2 nm iiii nm2 nm be FIG. 2: (a) Constant-current image obtained over a Bi-rich crystal region. Four different defects (i)-(iv) can be recognized.They are imaged on the atomic scale by STM in the top row of (b) and compared to ab-initio calculations which, due to theelevated computational cost required by the large number of atoms, have been restricted over the areas identified by the hatchedboxes visible in the upper frames. This comparison between theory and experiment allows for their unambiguous identification(see text for details). (c) STS data evidence a Fermi level in close proximity of the top of the valence band (VB), an observationwhich allows to identify this crystal region as p-doped. E D indicates the Dirac energy of the surface state. (d) Constant-currentimage obtained over a Te-rich crystal region with four characteristic different defects labeled (v)-(viii). (e) Again comparisonof atomic scale images with ab-initio calculations allows to unequivocally position them in a well-defined atomic plane of thequintuple layer structure. (f) STS data identify this crystal region as n-doped due to the proximity of the Fermi level with theconduction band (CB). by a comparison of the experimental data with simulated STM images obtained from ab initio calculations [bottomrow of Fig. 2(b)]. Based on the good general agreement achieved between the experimental and theoretical images wecan make the following defect assignment: defect (i) is a Bi Te(I) antisite; defect (ii) is a Te surface vacancy (V
Te(I) );defect (iii) is a Te vacancy in the third layer (V
Te(III) ); and defect (iv) is a Bi antisite in the fifth layer (Bi
Te(V) ).Since all defects are located in Te layers this region is expected to be Te-poor. This is confirmed also by theoreticalcalculations which show that these defects are characterized by the lowest formation energy in Bi-rich regions [16].The electronic properties of this sample region can be characterized by an investigation of the local density of statesas inferred by the STS spectra reported in Fig. 2(c). Following the energy level positioning scheme adopted in Ref. [26]we obtain a valence band maximum close to the Fermi level, thereby indicating the p-doped character of this sampleregion.Investigation of the opposite side of the crystal reveals four other defects (v)-(viii) [Fig. 2(d)], the STM appearanceof which is quite different from the data presented in Fig. 2(a)-(b) before. Again the symmetry and lateral extensionof the various defects in atomic resolution data were used to estimate the position of the defects. Based on theseestimations ab initio
DFT calculations were performed. The comparison presented in Fig. 2(e) reveals a very reasonableagreement. Our results show that all four defects are located in the two Bi layers of the Bi Te quintuple layer, inagreement with the expected Te-rich stoichiometry. Namely, the following assignments were concluded: defect (v) is aTe Bi(II) antisite; defect (vi) is a Bi vacancy (V
Bi(II) ); defect (vii) is a Te
Bi(IV) antisite; and defect (viii) is a Bi vacancy,V
Bi(IV) . As for the p-doped region, our defects identification is consistent with the formation energies calculatedby Scanlon et al. [16], which show that in Bi Te grown under Te-rich conditions the lowest energy defects are Te Bi antisites and Bi vacancies, with the last ones playing a much less significant role. Comparison of the STS spectrashown in Fig. 2, panels (c) and (f), reveals that the different types of defects present in Bi- and Te- rich areas stronglyinfluence the respective electronic properties, i.e. Te-rich regions show a rigid spectral shift towards negative energieswith respect to Bi-rich areas. This negative shifts amounts to approximately 140 meV between the p- and n-dopedsurface areas of the Te-rich and Bi-rich parts of the crystal, respectively. This observation implies that a topological - r / W c m Temperature T /K B CA R e s i s t i v i t y holes electrons - - H a ll m ob ili t y m / c m V s - H a ll c oe ff. R / c m C H S eebe ck c oe ff. S / m V / K Coordinate x /mm-2000200-400400 bca -35-10-50-4050 10 20 30 40 50 60
300 K
ABC p n
77 K
FIG. 3: (a) Room temperature Seebeck coefficient S measured along the crystal growth direction x with indication of threepoints (A,B,C) in the p–n transition region where electrical resistivity ρ (inset) was measured. (b) Hall coefficient R H measured at three temperatures T . Current and the magnetic field direction was within and perpendicular to the cleavageplane, respectively. (c) Hall mobility µ with respect to the position x along the boule. PNJ is naturally present within our crystals.Figure 3(a) shows the room-temperature Seebeck coefficient S of a Bi Te single crystal measured at variouspositions (0–60 mm) along the crystal growth direction x . It continuously and slowly changes along the crystal rodwith a sign change occurring at a position of 35 mm. A positive Seebeck coefficient in the range of 250–100 µ V/Kis obtained for the left part of the crystal boule (starting from x = 0). Further right to this region, i.e. between30 mm and about 35 mm, the Seebeck coefficient rapidly drops and eventually becomes negative, indicating that thecrystal conductivity changes from p- to n-type at the point where S ∼ µ V/K. In the region of the p–n transition( S ∼ µ V/K) a potential Seebeck microprobe demonstrates the change of sign of the Seebeck coefficient on the 200 µ m distance, limited by the resolution of the method. At x ≈
47 mm a high value S = − µ V/K is obtained.This interpretation is also corroborated by Hall measurements presented in Fig. 3(b) which were taken at thebottom and top parts of the boule. They show opposite signs of Hall coefficients R H , thereby confirming the presenceof both p- and n-type regions. The left part of the boule (0–30 mm) shows metallic behavior with small positive Hallcoefficients leading to a weakly temperature-dependent p-type carrier concentration, p = (5 − × cm − . Atthe right part a similar metallic-like behavior was observed for n-type carriers ( n (cid:39) cm − ). The large carrierconcentration and its weak temperature dependence indicate that the bottom and top parts of the boule are heavilydoped by acceptors and donors, respectively. Due to a strong reduction of the carrier concentration by two ordersof magnitude ( ∼ cm − ) in the vicinity of the p–n junction, the Hall coefficient changes sign and drops down to nm a b iii R RF F iiiv vi d I / d U / a r b . un i t s lowhigh-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.50102030405060 E - E /eV F x / n m c
110 meV nm FIG. 4: (a) Overview topographic STM image of the Bi Te crystal taken in the p–n transition zone (scan range:1.5 µ m × µ m). The surface shows two characteristic regions each of which is several hundred nm wide; while the roughsurface region (R) contains numerous islands and clusters the other, the flat region, exhibits atomically flat terraces (F) thatare separated by unit cell step edges. The zoomed in image (b) was taken within a flat region. Defects characteristic for both,n- (i-iv) and p-type (v-viii) Bi Te can be found. (c) Sequence of color-coded tunneling spectra taken along a 60 nm long line(3 nm increment) within a flat region. The energetic position of the minimum, which is marked by a red point in each spectrum,shifts from just above the Fermi level (bottom) to about 110 meV, indicative for a p–n transition with a width of about 40 nm.The hatched lines are guides to the eye. −
45 cm C − at low temperatures. These data indicate that an intrinsic binary chalchogenide TI can be obtained bya minimization of the carrier concentration through a careful control of the alloy composition between Bi and Te.Electrical resistivity temperature dependencies were measured on the samples taken from positions A, B, and C [cf.Fig 3(a)], i.e. within a distance of 2 mm of the sign change of the Hall coefficient. In all the points resistivity decreasesfor temperatures from RT to T ≈
30 K, and then remains almost constant down to the liquid helium temperatures forp- and n- type samples, while for point B resistivity increases, demonstrating semiconducting behavior. In both p- andn-parts of the crystal, the mobilities show a high average value of 10 cm V − s − at low temperatures [Fig 3(c)]. At x = 45 mm in the n-type region, i.e. close to the p–n transition, µ reaches the maximum value of about 10 cm V − s − .At the same position a maximum Seebeck coefficient S ∼ − µ V/K is observed, meaning that near the p-n transitionthere is the suppression of the defects formation.In order to investigate the minimal width of intrinsic TI p–n transitions that result from the slightly gradedstoichiometry during the crystallization process, we have performed STM/STS in the zone where the transition wasdetected in the transport measurements described above. An overview scan is displayed in Fig. 4(a). While the surfaceregion in the left and the bottom right part of the image is relatively rough (marked R) because of numerous smallislands and clusters, we also find regions with extended atomically flat terraces separated by unit-cell step edges (F).As we zoom into such an atomically flat region [Fig. 4(b)] we can observe the coexistence of defects which we identifyas being characteristic for n- and p-type Bi Te , i.e. defects (i), (ii), and (iii) as well as (v) and (vi), respectively(see above). This microscopic picture explains the drop of the Hall coefficient observed in transport, see Fig. 3 andRef. 18, 23. It is also in agreement with the mobility drop detected in the p-n junction region, which can be thus bedirectly linked to the increased number of total defects. Overall, these findings have deep implications going beyondthe scope of the present work. Minimizing bulk conductivity has always been considered as a powerful strategy to lettopological surface states dominate the scene. However, it was not clear up to know if this effect should be ascribed todefects suppression or compensation. The combination of transport measurements with atomic scale characterizationtechniques provides thus a detailed microscopic picture solving this problem, identifying defect compensation as theresponsible mechanism. Defect suppression can be safely ruled out since it would require a significant reduction ofthe defects concentration, a condition not consistent with our results.Local spectroscopic data show that these regions contain several PNJs which are very sharp. This is corroboratedby the color-coded STS data shown in Fig. 4(c) measured with a 3 nm increment along a 60 nm long line in a flatregion. The minimum position is marked by a red dot in each spectrum. Indeed we observe a shift of the minimumfrom close to the Fermi level E F at x ≤ E − E F ≈
110 meV at x ≥
50 nm, i.e. the transition from p- ton-doped Bi Te takes place over the remarkably short distance of about 40 nm, although the p- or n-type characterobserved at the two ends of the boule [cf. Fig. 2(c) and (f)] is not fully reached.Experimental data obtained on a larger length scale indicate that the doping gradient is not monotonic but exhibitssome fluctuations which lead to the coexistence of several p–n transition regions. Their distance typically amounts to ≈ µ m, a value also matching the length scale of flat (F) and rough (R) surface patches in Fig. 4(a). Both behaviorspotentially reflect the complicated convection processes that take place during the slow cooling and crystallizationprocess of the boule. Their averaged random distribution result in a Ohmic characteristic preventing spatially av-eraging technique from verifying the existence of a rectifying behavior, an aspect which may be directly tackled byusing nanoscale four probes techniques. Although their creation is not yet fully controlled, our data indicate that p–ntransition regions naturally resulting from the slightly graded stoichiometry of crystals may be very narrow. Furtherresearch may lead to optimized procedures which better control their position and width to make them useful forapplication.In summary, we demonstrated that, contrary to other materials, the ambipolar behavior of Bi Te makes possible,by taking advantage of the intrinsic defects which are unavoidably introduced during the growth process, to createtopological p–n junctions. Scanning tunneling spectroscopy reveals that their width amounts to 40 nm only. Therelated potential drop of approximately 110 meV shall allows for direct applications of PNJs in room temperaturedevices. Additionally, the defect compensation naturally achieved at the p–n interface results in a strong reductionof the bulk carriers concentration, thus paving the way to explore the PNJ topological properties. ACKNOWLEDGMENTS
The authors gratefully acknowledge stimulating discussions with Mario Italo Trioni. This work was supported bythe Deutsche Forschungsgemeinschaft within SPP 1666 (Grant No. BO1468/21-1). K.A.K. and O.E.T. acknowledgethe financial support by the RFBR (Grant nos. 13-02-92105 and 14-08-31110).
EXPERIMENTAL SECTION Bi Te crystals have been grown by the modified Bridgman technique with a temperature gradient of about 10 K/cmat the front of crystallization [24]. As-grown ingots had a single crystalline structure and were split into two partsalong the cleavage plane (0001) oriented along the growth direction (Fig. 1). One part of each crystal was cutperpendicular to the growth axis into 0.5–1 mm samples. Indium solder Ohmic contacts were used for transportmeasurements. The Hall resistance R yx and the resistance R xx were measured in the Hall bar geometry using astandard six-probe method on rectangular samples. A potential Seebeck microprobe was used to investigate theroom-temperature Seebeck coefficient in the region of p–n transition with a spatial resolution of 200 µ m. The Hallmobility was calculated from the measured conductivity and calculated carrier concentration, which was extractedfrom the measured Hall coefficient: n=1/(Rh*e). The carrier concentration dependence along the crystal was shownin [24].STM measurements were performed at a tip and sample temperature T = 4 . p < × − mbar and immediately inserted into the cryogenic STM. The PNJ was located by bringing thetip into tunneling distance from the surface and recording a local tunneling spectrum indicative for p- or n-type doping[cf. Fig. 2(c) and (f)]. After retracting the tip the sample was moved with an x − y -stage (initially by several tens tohundreds of µ m). This procedure was repeated until a spectral shift signaled that the boundary to the region governedby the other dopant had been crossed. Then a refined procedure with smaller x − y -movements was performed. THEORETICAL CALCULATION
Ab-initio theoretical calculations were performed by density functional theory, using a pseudo-potential representa-tion of the electron–ion interaction and local orbital basis sets, as implemented in the SIESTA code [27]. We used ageneralized gradient approximation (GGA) for the exchange-correlation functional [28] and a plane wave cutoff equalto 250 Ry. The experimental lattice constant for the hexagonal cell of Bi Te was used [29] and a 8 × Te in the direction orthogonal to the surface and alarge portion of vacuum of approximatively 30 ˚A. Constant–distance STM images have been simulated in the Tersoff-Hamann model [30] by calculating the Kohn-Sham local density of states (LDOS) in the energy interval between theFermi level E F and E F + eV b [31], where the applied bias eV b was fixed in agreement with the experimental setup.In particular the moderately large value of eV b , of the order of 0.4 eV, allows us to obtain numerically stable resultsand to reasonably neglect the small inaccuracy of the LDOS around the Fermi level due to the absence of spin-orbitinteraction in the calculation. A tip–surface distance equal to 3 ˚A was considered, although we verified that the resultof our STM simulations are robust against small variations of this value. The extension of the tip was taken intoaccount in the STM simulation by considering the mean value of LDOS within a radius of 4 ˚A around the tip position. ∗ corresponding author: [email protected][1] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, et al.,Science , 178 (2009), 10.1126/science.1173034.[2] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, et al., Nat. Phys. , 398(2009), http://dx.doi.org/10.1038/nphys1274.[3] Y. Jiang, Y. Wang, M. Chen, Z. Li, C. Song, K. He, L. Wang, X. Chen, X. Ma, and Q.-K. Xue, Phys. Rev. Lett. ,016401 (2012), http://link.aps.org/doi/10.1103/PhysRevLett.108.016401.[4] R. Venkatasubramanian, E. Siivola, T. Colpitts, and B. O’Quinn, Nature , 597 (2001),http://dx.doi.org/10.1038/35098012.[5] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045 (2010), http://link.aps.org/doi/10.1103/RevModPhys.82.3045.[6] P. Roushan, J. Seo, C. V. Parker, Y. S. Hor, D. Hsieh, D. Qian, A. Richardella, M. Z. Hasan, R. J. Cava, and A. Yazdani,Nature , 1106 (2009), http://dx.doi.org/10.1038/nature08308.[7] T. Zhang, P. Cheng, X. Chen, J.-F. Jia, X. Ma, K. He, L. Wang, H. Zhang, X. Dai, Z. Fang, et al., Phys. Rev. Lett. ,266803 (2009), http://link.aps.org/doi/10.1103/PhysRevLett.103.266803.[8] M. K¨onig, S. Wiedmann, C. Br¨une, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science , 096407 (2008), http://link.aps.org/doi/10.1103/PhysRevLett.100.096407.[10] X.-L. Qi, R. Li, J. Zang, and S.-C. Zhang, Science , 146802 (2010), http://link.aps.org/doi/10.1103/PhysRevLett.104.146802.[12] J. Wang, X. Chen, B.-F. Zhu, and S.-C. Zhang, Phys. Rev. B , 235131 (2012),http://link.aps.org/doi/10.1103/PhysRevB.85.235131.[13] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. , 620 (2006), http://dx.doi.org/10.1038/nphys384.[14] V. V. Cheianov, V. Fal’ko, and B. L. Altshuler, Science , 066402 (2009),http://link.aps.org/doi/10.1103/PhysRevLett.103.066402.[16] D. O. Scanlon, P. D. C. King, R. P. Singh, A. de la Torre, S. M. Walker, G. Balakrishnan, F. Baumberger, and C. R. A.Catlow, Advanced Materials , 2154 (2012), http://dx.doi.org/10.1002/adma.201200187.[17] L.-L. Wang, M. Huang, S. Thimmaiah, A. Alam, S. L. Bud’ko, A. Kaminski, T. A. Lograsso, P. Canfield, and D. D.Johnson, Phys. Rev. B , 125303 (2013), http://link.aps.org/doi/10.1103/PhysRevB.87.125303.[18] J.-L. Mi, M. Bremholm, M. Bianchi, K. Borup, S. Johnsen, M. Søndergaard, D. Guan, R. C. Hatch, P. Hofmann, andB. B. Iversen, Advanced Materials , 889 (2013), http://dx.doi.org/10.1002/adma.201203542.[19] Z. Zeng, T. A. Morgan, D. Fan, C. Li, Y. Hirono, X. Hu, Y. Zhao, J. S. Lee, J. Wang, Z. M. Wang, et al., AIP Advances , 072112 (2013), URL http://scitation.aip.org/content/aip/journal/adva/3/7/10.1063/1.4815972 .[20] Y. Zhao, C.-Z. Chang, Y. Jiang, A. DaSilva, Y. Sun, H. Wang, Y. Xing, Y. Wang, K. He, X. Ma, et al., Scientific Reports , 3060 EP (2013), URL http://dx.doi.org/10.1038/srep03060 .[21] Y. Zhao, H. Liu, X. Guo, Y. Jiang, Y. Sun, H. Wang, Y. Wang, H.-D. Li, M.-H. Xie, X.-C. Xie, et al., Nano Letters ,5244 (2014), pMID: 25102289, http://dx.doi.org/10.1021/nl502220p, URL http://dx.doi.org/10.1021/nl502220p .[22] G. Wang, X.-G. Zhu, Y.-Y. Sun, Y.-Y. Li, T. Zhang, J. Wen, X. Chen, K. He, L.-L. Wang, X.-C. Ma, et al., AdvancedMaterials , 2929 (2011), http://dx.doi.org/10.1002/adma.201100678.[23] S. Jia, H. Beidenkopf, I. Drozdov, M. K. Fuccillo, J. Seo, J. Xiong, N. P. Ong, A. Yazdani, and R. J. Cava, Phys. Rev. B , 165119 (2012), http://link.aps.org/doi/10.1103/PhysRevB.86.165119.[24] K. A. Kokh, S. V. Makarenko, V. A. Golyashov, O. A. Shegai, and O. E. Tereshchenko, CrystEngComm , 581 (2014),http://dx.doi.org/10.1039/C3CE42026D.[25] Y. Jiang, Y. Y. Sun, M. Chen, Y. Wang, Z. Li, C. Song, K. He, L. Wang, X. Chen, Q.-K. Xue, et al., Phys. Rev. Lett. , 066809 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.108.066809 .[26] P. Sessi, M. M. Otrokov, T. Bathon, M. G. Vergniory, S. S. Tsirkin, K. A. Kokh, O. E. Tereshchenko, E. V. Chulkov, andM. Bode, Phys. Rev. B , 161407 (2013), http://link.aps.org/doi/10.1103/PhysRevB.88.161407.[27] J. M. Soler, E. Artacho, G. J. D., A. Garc´ıa, J. Junquera, P. Ordej´on, and D. S´anchez-Portal, J. Phys.: Cond. Matter ,2745 (2002).[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[29] G. Wang and T. Cagin, Phys Rev. B , 075201 (2007).[30] J. Tersoff and D. R. Hamann, Phys Rev. B , 805 (1985).[31] F. Donati, P. Sessi, S. Achilli, A. Li Bassi, M. Passoni, C. S. Casari, C. E. Bottani, A. Brambilla, A. Picone, M. Finazzi,et al., Phys Rev. B79