Theory of electronic and spin-orbit proximity effects in graphene on Cu(111)
We study orbital and spin-orbit proximity effects in graphene adsorbed to the Cu(111) surface by means of density functional theory (DFT). The proximity effects are caused mainly by the hybridization of graphene and copper orbitals. Our electronic structure calculations agree well with the experimentally observed features. We carry out a graphene–Cu(111) distance dependent study to obtain proximity orbital and spin-orbit coupling parameters, by fitting the DFT results to a robust low energy model Hamiltonian. We find a strong distance dependence of the Rashba and intrinsic proximity induced spin-orbit coupling parameters, which are in the meV and hundreds of eV range, respectively, for experimentally relevant distances. The Dirac spectrum of graphene also exhibits a proximity orbital gap, of about 20 meV. Furthermore, we find a band inversion within the graphene states accompanied by a reordering of spin and pseudospin states, when graphene is pressed towards copper.
Copper is an important material for graphene. Graphene–copper junctions are often encountered in technological applications.Giovannetti et al. (2008); Khomyakov et al. (2009) For example, graphene can be used to seal a copper surface to preserve its excellent plasmonic characteristics.Kravets et al. (2014) The growth of graphene via CVD by the deposition of CH on copper surfaces is amongst the most popular techniques to obtain large (poly)crystalline graphene. Gao et al. (2010) Even single layer graphene grains of millimeter size as well as pyramid-like bi- and trilayer graphene, hexagonal onion ring-like graphene grains can be grown on copper.Yan et al. (2014); Boyd et al. (2015) Important to our study, graphene produced on a copper surface exhibits a giant spin Hall effect, lilkely due to residual cupper adatoms and ad-clusters.Balakrishnan et al. (2014)
Experimentally, graphene on the Cu(111) surface has been well studied by means of angle-resolved photoemission spectroscopy (ARPES) Avila et al. (2013); Jeon et al. (2013); Marsden et al. (2013); Shikin et al. (2013); Varykhalov et al. (2010); Vita et al. (2014); Walter et al. (2011) and scanning tunneling microscopy (STM).Gao et al. (2010) The linear dispersion of graphene is found to be preserved. ARPES measurements of the graphene–Cu(111) system find, that graphene is getting electron doped,Walter et al. (2011) leading to a shift of the Dirac energy , which we define as the average energy of the graphene state energies at K with respect to the Fermi energy . Typically, is about eV with respect to the Fermi energy.Walter et al. (2011) The top of the band edge of copper begins at eV below the Fermi level. It is observed that a gap opens within the Dirac cone of graphene of about 50–180 meV.Avila et al. (2013); Marsden et al. (2013); Shikin et al. (2013); Varykhalov et al. (2010); Vita et al. (2014); Walter et al. (2011)
Spin-orbit coupling (SOC) effects in graphene on selected metal substrates were studied theoretically Abdelouahed et al. (2010); Li et al. (2011) and experimentally Rader et al. (2009); Marchenko et al. (2012) and it was noticed that substrates can induce sizeable spin-orbit effects important for spintronics applications. Žutić et al. (2004); Han et al. (2014) Spin resolved ARPES experimentsShikin et al. (2013) focused on the spin-orbit effects introduced by metallic surfaces in graphene, investigating the role of the atomic number of the substrate. It was found that the states of graphene can be split due to Rashba spin-orbit coupling by up to 100 meV in the case of Au and 10 meV in the case of Ni,Li et al. (2011) respectively. Copper substrate induced spin-orbit splittings in graphene are expected to be substantially smaller.Shikin et al. (2013) They were measured at a temperature of 40 K, which gives a resolution limit and also the upper bound for the spin-orbit effects of 3.4 meV. The mechanism introducing the spin-orbit interaction was identified to be the hybridization between substrate and graphene states. Shikin et al. (2013) Our present work agrees with this conclusion, and predicts the values of the Rashba splitting to be about 2 meV for a reasonable distance between graphene and copper, just below the stated experimental resolution of Ref. Shikin et al., 2013.
Crucial to obtain accurate graphene–metal distances is to consider van der Waals interactions. It was foundOlsen and Thygesen (2013); Vanin et al. (2010); Andersen et al. (2012), that the dispersive long-range interactions play an important role in binding, yielding graphene–copper distances of 2.91 to 3.58 Å.
Here, we focus on hybridization and proximity effects by means of DFT calculations. By the application of an effective Hubbard ,Anisimov et al. (1991) which corrects for self interaction errors, we achieve a good agreement with experiment in terms of the emission spectra and the band structure features. We carry out an analysis of the orbital composition of the band structure, giving us hints for a model Hamiltonian including spin-orbit interactions, which can be used to describe graphene in combination with many other materials that yield a or higher symmetric system. We then fit the DFT data to the model Hamiltonian and extract parameters such as the induced gap as well as spin-orbit coupling values. As the graphene–copper distance is not exactly known experimentally, and there is still a theoretical uncertainty in determining its magnitude, we carry out a distance-dependent study.
Our main finding is a strong graphene–Cu(111) distance-dependent spin-orbit coupling introduced in the graphene states. We use a model Hamiltonian to describe those states, for which we observe a Rashba spin-orbit coupling parameter which reaches values of meVs, while being absent in pristine graphene. The proximity induced intrinsic SOC is in the hundreds of eV range, a factor of ten larger than in pristine graphene. We also observe a closing of the induced gap for a graphene–copper distance of 2.4 Å. This is accompanied by a peculiar reordering of spin and pseudospin states associated with a gap inversion at small distances.
The paper is organized as follows. Sec. II deals with the computational methods used. Geometrical structure modeling is described in Sec. III.1. In Sec. III.2 we carry out the analysis of the band structure. In Sec. III.3 we introduce our model Hamiltonian and fit it to the ab-initio data. Finally, in Sec. III.4 we present our graphene–copper distance dependent study with a discussion of the proximity induced effects.
Ii Computation methods
We used DFT implemented in the plane-wave code quantum espresso.Giannozzi et al. (2009) The calculations were performed at a point sampling of if not indicated otherwise. A slab geometry was applied, where we added a minimum of 15 Å of vacuum around the structure in direction. We used the Kresse-Joubert ultrasoft (relativistic) PBEPerdew et al. (1996) projector augmented wave pseudopotentials.Kresse (1999) The plane wave energy cutoff was set to 40 Ry and the charge density cutoff to 320 Ry to ensure converged results. Van der Waals interactions were taken into account using the empirical method of Grimme.S.Grimme (2006) To cross check spin-orbit coupling calculations we also employed the all electron, full potential linearized augmented plane wave code wien2k.Blaha et al. (2001) We found that spin-orbit coupling splittings were differing at most by 10%. For processing our distance studies the atomic simulation environment (ASE)Bahn and Jacobsen (2002) was used. Hellmann-Feynman forces in relaxed structures were decreased until they were smaller than . Calculations for graphene on Cu(111) included calculations adding the Hubbard correction.Cococcioni and de Gironcoli (2005)
Iii Graphene Cu(111) study
iii.1 Choice of unit cell
The mismatch between Cu(111)’s surface lattice constant of Å Straumanis and Yu (1969) and graphene’s lattice constant of 2.46 Å is 3.8%. STM experimentsGao et al. (2010) observe regions with different moiré structures; the most observed one (30%) is a commensurate lattice configuration with a periodicity of 66 Å. Another experimentAvila et al. (2013) found that 60% of graphene grains on Cu(111) are preferentially rotated by with respect to the substrate. To account for the lattice mismatch, one would have to chose a unit cell, which is computationally very demanding, containing hundreds of atoms. We set the lattice constant of copper to be compatible with the experimental graphene lattice constant, to describe graphene as realistically as possible, following Ref. Khomyakov et al., 2009. A supporting fact to use graphene’s lattice constant is that graphene does not chemically bind to copper and its strong in-plane bonds remain intact.
In the transverse direction to the Cu(111) surface one distinguishes three non-equivalent Cu planes. We label the planes from the surface towards bulk as top, hcp and fcc, see Fig. 1. We tested three different commensurate configurations named according to which Cu layer the carbon atoms sit over. This gives rise to three possible graphene physisorbed positions named as top-fcc, top-hcp, and fcc-hcp configuration.Khomyakov et al. (2009) In Fig. 1 we show the top-fcc configuration, where one carbon atom (say from sublattice B) is on top of a copper atom of the top layer, while the other carbon (from sublattice A) is over the fcc Cu layer.
In general, the graphene sublattices have different environments. This breaks the sublattice symmetry of graphene and results in sublattice resolved spin-orbit coupling effects.Gmitra et al. (2013); Gmitra and Fabian (2015) To simulate a copper surface we used four layers of copper. We checked that the physics of the graphene low energy states does not change upon increasing the number of layers. In addition, we found good agreement of the band structure with experiment.Avila et al. (2013); Jeon et al. (2013); Marsden et al. (2013); Shikin et al. (2013); Varykhalov et al. (2010); Vita et al. (2014); Walter et al. (2011)
In our studies we first relaxed the copper slab alone without van der Waals corrections and then fixed its degrees of freedom and let just the carbon atoms relax in direction including empirical van der Waals corrections.S.Grimme (2006) To start with, the copper slab is strained in the plane such that its surface lattice constant is the same as the experimental graphene lattice constant of 2.46 Å yielding an effective bulk lattice constant of Å. This represents a compression of the copper slab by 3.8% with respect to the bulk value of 3.61 Å.Straumanis and Yu (1969) After letting the copper slab relax in direction, the distance of copper atoms from plane to plane was 2.59 Å, corresponding to an expansion of 1.7% compared to bulk copper. This compensates to some extent for the compression in the plane.
Comparing the top-fcc with the other commensurate configurations top-hcp and fcc-hcp we found slightly different graphene–Cu(111) distances of 3.10 Å, 3.11 Å, and 3.12 Å, respectively. The corrugation of the carbon atoms in direction is less than Å, expressing the weak nature of binding. The lowest energetic configuration is the top-fcc arrangement, followed by the top-hcp, which is only 2.3 meV higher in energy per unit cell. The highest one in total energy with 12.3 meV compared to top-fcc is fcc-hcp, where the nearest copper atom sits within the carbon ring. Therefore in the following study we consider the top-fcc configuration.
iii.2 Choice of methods and electronic structure
The orbital resolved electronic structure of graphene on Cu(111) is shown in Fig. 2 for DFT+ with an effective eV Cococcioni and de Gironcoli (2005) acting on the Cu 3 electrons for a copper–graphene distance of 3.09 Å. It can be seen, that the Dirac cone structure is preserved for energies higher than eV. Below this energy region the graphene states hybridize with the copper states. This can be seen by the avoided crossings if one follows the band towards the point at eV. On this way, at eV the states branch and strongly hybridize with a copper band consisting of and states. Those are states which are situated on the surfaces of the slab and whose degeneracy is broken due to the graphene potential. The graphene states starting from eV at the point are mainly unaffected. The copper and states are present in the energy region between eV and eV as well as from eV and upwards. The copper band structure obtained here is qualitatively in agreement with bulk fcc calculations. The position of the Fermi energy was converged for a dense sampling of the Brillouin zone.
There is charge transfer from the Cu(111) surface to graphene. As a result, graphene gets -doped.Khomyakov et al. (2009) The Dirac cone is shifted below the Fermi level by about eV according to experiments. We compared the effect of the relaxation of copper slabs in -direction for the relaxed and non-relaxed (bulk lattice constant of Å) cases on the doping of graphene. For non-relaxed (compressed) copper slabs the electron doping of graphene was significantly higher than for relaxed slabs due to the higher kinetic energy in non-relaxed copper. For relaxed copper slabs the Dirac energy shift is comparable to experiment,Gao et al. (2010) being 350 meV.
To account for correlation effects, we applied a Hubbard correction of 1 eV to the copper states. In this way we match the onset of the completely filled copper levels, which show up at eV below the Fermi energy in the ARPES experiments.Avila et al. (2013); Marsden et al. (2013); Shikin et al. (2013) The effect of the Hubbard correction is a rigid shift of the filled copper levels to lower energies without changing their band widths. However, we see a strong dependence of the copper level energies on the compression of the copper slab, they are 1 eV higher in energy for the compressed than for the relaxed one. The proper position of the levels is significant for the spin-orbit coupling induced proximity effect in the Dirac cone, as there can be larger hybridization, when the levels are closer to the states of interest.
All in all we find a good agreement of the band structure with experiment.Avila et al. (2013) The only shortcoming is the description of the graphene gap, which is opening at the Dirac energy . We find it to be 20 meV, which is lower than the 50 to 180 meV stated in experiments.Avila et al. (2013); Marsden et al. (2013); Shikin et al. (2013); Varykhalov et al. (2010); Vita et al. (2014); Walter et al. (2011) This deviation could be due to the limitations of semilocal and local exchange-correlation functionals.
iii.3 Model Hamiltonian
As we demonstrated above, DFT+ reasonably captures the electronic structure of graphene on the Cu(111) surface. Now we use the first-principles calculations to predict proximity induced effects of the copper surface on the spin-orbit coupling in graphene. For this purpose we study a Hamiltonian describing the low energy states of graphene on Cu(111). The Hamiltonian contains orbital and spin-orbit coupling parts and describes graphene whose symmetry point group is lowered from (pure graphene) to . Such a Hamiltonian was introduced already in the context of hydrogenated graphene Gmitra et al. (2013) in which the pseudospin symmetry gets broken explicitly by hydrogenation, but it was also found useful in graphene whose pseudospin is broken implicitly only, by placing graphene on incomensurate lattices such as MoS. Gmitra and Fabian (2015) In our case the pseudospin symmetry is broken explicitly as the pseudospin state is well defined but the two sublattices experience a different orbital environment, see Fig. 1. This proximity Hamiltonian has the form,
where is the Fermi velocity and labels the valley degree of freedom. and are the Cartesian components of the electron wave vector measured from K(K), and are the pseudospin Pauli matrices acting on the two-dimensional vector space formed by the two triangular sublattices of graphene. The first term in describes gapless Dirac states. The second term describes the effective orbital hybridization energy, which acts as a staggered potential on sublattices A and B, where is the pseudospin Pauli matrix and is the unit matrix in spin space. This Hamiltonian term leads to an orbital proximity induced gap in the Dirac spectrum of . This gap is still present even when spin-orbit coupling is turned off. A consequence of the pseudospin inversion asymmetry is the sublattice-resolved intrinsic spin-orbit coupling. As intrinsic spin-orbit coupling is a next-nearest neighbor hopping, it acts solely on a given sublattice. We describe it with parameters and for sublattice A and B, respectively. We denote by the spin Pauli matrix and by the unit matrix acting on the pseudospin space. If , the spin degeneracy gets lifted already by this intrinsic term, reflecting the loss of space inversion symmetry. The space inversion asymmetry itself gives rise to Rashba type spin-orbit coupling whose strength is measured by , which is a nearest-neighbor spin-flip hopping, contributing further to the spin splitting of the low energy bands.
The four eigenvalues of the model Hamiltonian at the K point () read
where and for compactness. We ordered the eigenvalues by decreasing energies, where we assumed . The eigenstates and always have spin- expectation values of and , and pseudospin- expectation values of and and are localized on sublattice B and A, respectively. The eigenstates with and in general are mixtures of sublattices and spin directions, but have almost , and , under the assumption that , . In the model Hamiltonian there are four unknown parameters. To construct a set of independent equations we also take into account the spin- expectation value for the first eigenstate denoted by . The model parameters thus can be expressed as follows
We note that special care has to be taken when associating the order of the DFT eigenvalues with respect to the model Hamiltonian eigenvalues. For every state we compared the sublattice localization and values for both the DFT and model calculations.
In Fig. 3 we compare low energy graphene bands calculated from DFT and model, for a distance of . The fitted model parameters are , meV, meV, and .
The proximity effects, both, the orbital and spin-orbit coupling ones are significant. The hybridization gap dominates the energy scale. It yields a gap value of meV. The Rashba spin-orbit coupling parameter of 1.2 meV indicates a very strong effect of the space inversion asymmetry, which would correspond to a transverse electric field of 240 V/nm for bare graphene.Gmitra et al. (2009) The intrinsic spin-orbit coupling parameters have opposite sign and their amplitudes are significantly enhanced in comparison to the tens of eV in bare graphene.Gmitra et al. (2009) The Fermi velocity is (equivalent to a nearest neighbor hopping of 2.55 eV). We see that the band structure is isotropic in this range of points and the model description agrees very well with the DFT data. We observed a good agreement up to energies eV away from the Dirac energy.
We also compare the band spin splittings of the valence and conduction bands, see Fig. 4.
It can be seen, that by construction, the splittings at K are described exactly. The model reproduces very well the narrowing of the band splittings for points up to away from the K point even though only information from the K point enters. As the model does not include spin-orbit coupling terms dependent on , both the valence and conduction band splittings from the model calculations saturate at a common value for larger due to the Rashba SOC. To include dependent terms one needs to consider terms such as pseudospin inversion asymmetry (PIA)Gmitra and Fabian (2015); Gmitra et al. (2013); Irmer et al. (2015) which can capture the dependence of the splittings. In the DFT calculations we observed that the splittings for valence (conduction) bands increase (decrease) with larger distances from K as the interaction with copper levels increases (decreases) and the induced spin-orbit effects are stronger (weaker).
iii.4 Distance study
Standard DFT can not account for dispersive forces. Different methods dealing with van der Waals effects often yield inconsistent resultsVanin et al. (2010); Olsen and Thygesen (2013); Andersen et al. (2012) when trying to treat graphene on metal surfaces. Therefore, we conduct calculations of electronic properties for different graphene–Cu(111) distances. We used the Hubbard correctionCococcioni and de Gironcoli (2005) with eV for Cu electrons. The relative coordinates of the atoms within the copper slab and within graphene were fixed and the graphene–copper distance was varied. We apply the same analysis as in Sec. III.3 for each distance configuration and extract the total energy of the structure, the Dirac energy shift , the hybridization gap , the Rashba and intrinsic spin-orbit coupling parameters as well as spin- expectation values of the graphene states at the K point.
In Fig. 5(a) we show the total energy as a function of the graphene distance from the Cu(111) surface. The curve is shifted with respect to the minimal total energy at the distance of 3.09 Å. The energy dependence has a rather shallow minimum where the energy increases by just 0.5 eV when graphene is pushed to a distance of 2 Å.
Fig. 5(b) visualizes the shift of the Dirac energy with respect to the Fermi level. We see that graphene stays -doped for distances smaller than 3.5 Å, and the curve has two regimes. For larger distances down to 2.5 Å there is a linear behavior with a positive slope, the more graphene is pushed towards the Cu(111) surface, the more -doped it gets. For distances smaller as 2.5 Å the slope reverses its sign and is more shallow. This means that there occurs a significant charge transfer from the copper slab to the graphene sheet, which saturates at smaller distances.
Figure 5(c) shows the values for the proximity induced potential and the Rashba spin-orbit parameter . The Rashba parameter is increasing steadily with decreasing distance. We also plot the derivative of the Rashba parameter with respect to the distance . One sees that the Fermi level shift and the change in the Rashba parameter are correlated by comparing the derivative of the Rashba parameter to the Fermi level shift. Both curves change their trend at 2.5 Å. We can see that the origin of the Rashba spin-orbit coupling is due to charge doping (determined by the Fermi energy shift ), leading to a built-in electric field, and due to the positioning of the graphene sheet in the electrostatic potential of the Cu(111) surface. At the distance of 2.5 Å the charge doping stops, and therefore the Rashba spin-orbit coupling increases at a lower pace. It remains increasing though, as the graphene sheet resides in a potential which becomes steeper as it gets closer to the nuclei of copper. It is surprising that the , which first increases from larger to smaller distances, decreases, becomes zero at 2.4 Å and then inverts its sign. We will discuss this in more detail later. We estimate the pressure one would have to exert on graphene to reach this distance as
where is the energy difference between the lowest energetic state and the state where the transition happens, their distance difference, and is the area of the unit cell. The bulk modulus of copper for comparison is 184 GPa.Wyckoff (1976)
The amplitudes of the intrinsic spin-orbit coupling parameters and strongly increase as graphene is pushed towards the Cu(111) surface, see Fig. 5(d). For large distances both parameters tend to values comparable in size as in pure graphene. For smaller distances the sublattice asymmetry transfers to the parameters and is much stronger affected due to the specific graphene sublattice positioning on Cu(111). reaches values up to 7 meV, whereas stays smaller than 1 meV for all tested distances and tends to saturate at 1 meV when reaching a small distance of 1.8 Å.
The last two panels in Fig. 5(e) and (f) show the spin- expectation values at K for the eigenvalues , where 1 labels the lowest energy and 4 the highest energy state. From Fig. 5(e) we can see that the outermost expectation values represent spin states of mixed spin, as values of spin 1/2 are only reached, when graphene is well separated from copper. The spin expectation values of states 2 and 3 are pure states and are always quantized in direction. When the hybridization gap closes, at 2.4 Å, we observe, that the signs of all spin expectation values change abruptly. This behavior is exemplified in Fig. 6. When the distance of graphene to copper is decreased, the spin as well as pseudospin signs change.
In Fig. 7 we show the topology of the bands obtained from DFT calculations around K for distances of 2.2 and 2.5 Å, with the corresponding spin- expectation values. The plot is consistent with Fig. 5(e) and (f), for 2.5 Å the band structure resembles the one in Fig. 3 and has spin up-down-up-down sequence, where the inner eigenstates have pure components. The spin- character within the bands stays the same. The band structure topology for 2.2 Å is different. At the K point the inner eigenstates again have pure spin, but all signs are reversed. Furthermore, the spin- character is not preserved within the bands. There is evidence for a band inversion for the inner bands with a significant spin mixing to outermost bands. The spin reversal is accompanied by a change of the pseudospin character of the states. The valence states become localized on the A sublattice and conduction bands on sublattice B. We note that similarly to the spin mixing for the outermost bands, the states are also sublattice mixed, which is also depicted in Fig. 6. Our model is able to reproduce the spin- behavior of Fig. 7 (not shown here).
We have shown that the electronic band structure measured by ARPES is reasonably described by DFT+ calculations. We are able to correctly describe the Fermi level position and copper band onset. Based on this good orbital description we predict the spin-orbit coupling effects by analyzing the low energy graphene states using a robust model Hamiltonian. We show, that our Hamiltonian is able to describe the spin-orbit induced band splittings even away from the K point. We extracted spin-orbit coupling parameters as well as spin expectation values dependent on the graphene–copper distance and found a strong distance-dependent behavior of spin-orbit coupling parameters and a reordering of the spin and pseudospin structure at the Dirac point at =2.4 Å. At low distances the Dirac band structure gets inverted due to the overlap of opposite spin valence and conduction bands. Our findings are experimentally verifiable with techniques such as ARPES, by increasing the resolution to resolve the meV and sub meV spectral ranges.
Acknowledgements.This work was supported by the DFG SFB Grant No. 689 and GRK Grant No. 1570, and by the EU Seventh Framework Programme under Grant Agreement No. 604391 Graphene Flagship.
- Giovannetti et al. (2008) G. Giovannetti, P. Khomyakov, G. Brocks, V. Karpan, J. van den Brink, and P. Kelly, Phys. Rev. Lett. 101, 026803 (2008).
- Khomyakov et al. (2009) P. Khomyakov, G. Giovannetti, P. C. Rusu, G. Brocks, J. van den Brink, and P. J. Kelly, Phys. Rev. B 79, 195425 (2009).
- Kravets et al. (2014) V. G. Kravets, R. Jalil, Y.-J. Kim, D. Ansell, D. E. Aznakayeva, B. Thackray, L. Britnell, B. D. Belle, F. Withers, I. P. Radko, Z. Han, S. I. Bozhevolnyi, K. S. Novoselov, A. K. Geim, and A. N. Grigorenko, Sci. Rep. 4, 5517 (2014).
- Gao et al. (2010) L. Gao, J. R. Guest, and N. P. Guisinger, Nano Lett. 10, 3512 (2010).
- Yan et al. (2014) Z. Yan, Z. Peng, and J. Tour, Acc. Chem. Res. 47, 1327 (2014).
- Boyd et al. (2015) D. Boyd, W.-H. Lin, C.-C. Hsu, M. Teague, C.-C. Chen, Y.-Y. Lo, W.-Y. Chan, W.-B. Su, T.-C. Cheng, C.-S. Chang, C.-I. Wu, and N.-C. Yeh, Nat. Commun. 6, 6620 (2015).
- Balakrishnan et al. (2014) J. Balakrishnan, G. K. W. Koon, A. Avsar, Y. Ho, J. H. Lee, M. Jaiswal, S.-J. Baeck, J.-H. Ahn, A. Ferreira, M. a. Cazalilla, A. H. C. Neto, and B. Özyilmaz, Nat. Commun. 5, 4748 (2014).
- Avila et al. (2013) J. Avila, I. Razado, S. Lorcy, R. Fleurier, E. Pichonat, D. Vignaud, X. Wallart, and M. C. Asensio, Sci. Rep. 3, 2439 (2013).
- Jeon et al. (2013) C. Jeon, H.-N. Hwang, W.-G. Lee, Y. G. Jung, K. S. Kim, C.-Y. Park, and C.-C. Hwang, Nanoscale 5, 8210 (2013).
- Marsden et al. (2013) A. J. Marsden, M.-C. Asensio, J. Avila, P. Dudin, A. Barinov, P. Moras, P. M. Sheverdyaeva, T. W. White, I. Maskery, G. Costantini, N. R. Wilson, and G. R. Bell, Rap. Res. Lett. 7, 643 (2013).
- Shikin et al. (2013) A. M. Shikin, A. G. Rybkin, D. Marchenko, A. a. Rybkina, M. R. Scholz, O. Rader, and A. Varykhalov, New J. Phys. 15, 013016 (2013).
- Varykhalov et al. (2010) A. Varykhalov, M. R. Scholz, T. K. Kim, and O. Rader, Phys. Rev. B 82, 121101(R) (2010).
- Vita et al. (2014) H. Vita, S. Böttcher, K. Horn, E. N. Voloshina, R. E. Ovcharenko, T. Kampen, A. Thissen, and Y. S. Dedkov, Sci. Rep. 4, 5704 (2014).
- Walter et al. (2011) A. L. Walter, S. Nie, A. Bostwick, K. S. Kim, L. Moreschini, Y. J. Chang, D. Innocenti, K. Horn, K. F. McCarty, and E. Rotenberg, Phys Rev. B 84, 195443 (2011).
- Abdelouahed et al. (2010) S. Abdelouahed, A. Ernst, J. Henk, I. V. Maznichenko, and I. Mertig, Phys. Rev. B 82, 1 (2010).
- Li et al. (2011) Z. Y. Li, Z. Q. Yang, S. Qiao, J. Hu, and R. Q. Wu, J. Phys. Condens. Matter 23, 225502 (2011).
- Rader et al. (2009) O. Rader, A. Varykhalov, J. Sánchez-Barriga, D. Marchenko, A. Rybkin, and A. Shikin, Phys Rev. Lett. 102, 057602 (2009).
- Marchenko et al. (2012) D. Marchenko, A. Varykhalov, M. R. Scholz, G. Bihlmayer, E. I. Rashba, A. Rybkin, a. M. Shikin, and O. Rader, Nat. Commun. 3, 1232 (2012).
- Žutić et al. (2004) I. Žutić, J. Fabian, and S. D. Sarma, Rev. Mod. Phys. 76, 323 (2004).
- Han et al. (2014) W. Han, R. K. Kawakami, M. Gmitra, and J. Fabian, Nat. Nano 9, 794 (2014).
- Olsen and Thygesen (2013) T. Olsen and K. Thygesen, Phys. Rev. B 87, 075111 (2013).
- Vanin et al. (2010) M. Vanin, J. J. Mortensen, A. K. Kelkkanen, J. M. Garcia-Lastra, K. S. Thygesen, and K. W. Jacobsen, Phys. Rev. B 81, 081408(R) (2010).
- Andersen et al. (2012) M. Andersen, L. Hornekaer, and B. Hammer, Phys. Rev. B 86, 085405 (2012).
- Anisimov et al. (1991) V. Anisimov, J. Zaanen, and O. Andersen, Phys. Rev. B 44 (1991).
- Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matter 21, 395502 (2009).
- Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Kresse (1999) G. Kresse, Phys. Rev. B 59, 1758 (1999).
- S.Grimme (2006) S.Grimme, J. Comp. Chem. 27, 1787 (2006).
- Blaha et al. (2001) P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, WIEN2k, An Augmented Plane Wave+Local Orbitals Program for Calculating Crystal Properties, edited by K. Schwarz, Vol. 1 (Technische Universität Wien, Austria, 2001).
- Bahn and Jacobsen (2002) S. Bahn and K. W. Jacobsen, Comp. Sci. Eng. 4, 56 (2002).
- Cococcioni and de Gironcoli (2005) M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105 (2005).
- Straumanis and Yu (1969) M. E. Straumanis and L. S. Yu, Acta Cryst. Sect. A 25, 676 (1969).
- Gmitra et al. (2013) M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. Lett. 110, 246602 (2013).
- Gmitra and Fabian (2015) M. Gmitra and J. Fabian, Phys. Rev. B 92, 155403 (2015).
- Gmitra et al. (2009) M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and J. Fabian, Phys. Rev. B 80, 235431 (2009).
- Irmer et al. (2015) S. Irmer, T. Frank, S. Putz, M. Gmitra, D. Kochan, and J. Fabian, Phys. Rev. B 91, 115141 (2015).
- Wyckoff (1976) R. W. G. Wyckoff, Crystal Structure, 2nd ed. (Interscience, New York, 1976).