# Atomic-scale visualization of the p-d hybridization in III-V semiconductor surfaces doped with transition metal impurities

K. Badiane, G. Rodary, M. Amato, A. Gloter, C. David, H. Aubin, J.-C. Girard

► To cite this version:

K. Badiane, G. Rodary, M. Amato, A. Gloter, C. David, et al., Atomic-scale visualization of the p-d hybridization in III-V semiconductor surfaces doped with transition metal impurities. Physical Review B, 2022, 105 (23), pp.235443. 10.1103/PhysRevB.105.235443. hal-03763319

HAL Id: hal-03763319

https://hal.science/hal-03763319v1

Submitted on 7 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Atomic scale visualization of the *p*-*d* hybridization in III-V surface semiconductors doped with transition metal impurities

K. Badiane,¹ G. Rodary,¹ M. Amato,² A. Gloter,² C. David,¹ H. Aubin,¹ and J.-C. Girard¹

¹Centre de Nanosciences et de Nanotechnologies (C2N), CNRS,
Université Paris-Saclay, 10 Boulevard Thomas Gobert, 91120 Palaiseau, France

²Laboratoire de Physique des Solides, Université Paris-Saclay, CNRS, 91405, Orsay, France

(Dated: February 7, 2025)

The *p*-*d* hybridization of transition metal impurities in a semiconductor host is the mechanism that couples valence-band electrons and localized spins. We use scanning tunneling microscopy and spectroscopy combined with density functional theory to probe at the atomic scale hybridization of Cr single impurities with a GaAs host. Combining spatial density of states mapping and in-gap states spectroscopy of the Cr substituted at the surface of the semiconductor, we give a detailed picture of the spatial extension and the electronic structure of the strongly anisotropic wave function of Cr on GaAs(110). First principles calculations allow us to identify the electronic character and origin of each states and show that the main resonance peaks and the wave function with "drop-eye" lobes experimentally observed for 3d metal impurities in III-V semiconductors are direct local evidences of the *p*-*d* hybridization.

INTRODUCTION

The electronic properties of transition metal impuri-
ties introduced in a semiconductor host have been in-
tensively studied in the 1970s and 1980s for microelec-
tronics or optoelectronics applications (see, for instance,
Refs. [1, 2]). Interest in such systems grew during the
1990s because of the discovery of ferromagnetism in di-
luted magnetic semiconductors (DMSs). The origin of
magnetic order in these materials has been ascribed to
the carrier-mediated coupling between spins of doping
transition metal dopants, which has been described by
phenomenological models and first principles calcula-
tions [3-7]. These theoretical works have shown that the
main local cause of the ferromagnetic interaction in III-
V semiconductors between transition metal atoms is the
*p*-d hybridization between the *d* orbitals of the impurity
and the *p* orbitals of host atoms.

Scanning tunneling microscopy (STM) and spec-
troscopy (STS) contributed greatly to the understanding
of the physics of transition metal ions in semiconductors
by imaging at the atomic scale impurities deposited on
semiconductors surfaces or substituted to cations near
the surface. STM and STS studies have been carried
out on a GaAs surface for several elements: Mn [8–
15], Fe [13, 16, 17], and Co [13, 18]. Despite the dif-
ferent *d*-band fillings, the STS data show similar reso-
nances in the band-gap and the conductance maps at the
corresponding resonance voltages reveal extended wave-
functions with "crab-like" and "drop-eye" shapes. These
results have been discussed qualitatively following the
Mahadevan and Zunger approach developed for bulk ma-
terials [5], but no precise explanation of the origin of the
shapes and their relations to the contributions of the *p*
and *d* orbitals to the state has been so far given for sur-
faces. In spite of the large number of STM/STS stud-
ied showing similar results on transition metal impuri-

ties in III-V semiconductors surface, a complete descrip-
tion of the physical mechanism explaining the observed
wave function shape and resonance states is still missing.
In this paper, we aim to understand the origin of these
observations by shedding light, both experimentally and
theoretically, on the link between STS spectroscopy and
STM maps, i.e. between the electronic structure and spa-
tial distribution of the density of states. Our analysis of
the extended impurity states presents an intermediate de-
scription between the local chemical description, i.e. the
*d*-orbitals of the impurity hybridize with the *p*-orbitals of
the neighboring anions, and the physical description, i.e.
the *d*-orbitals hybridize with Bloch valence states at the
Γ point. Such a description will enable a more precise un-
derstanding of the kinetic exchange mechanism enabling
ferromagnetic ordering in those materials.

EXPERIMENTAL DETAILS

We explore by STM and STS the electronic states of
Cr single atoms substituted to Ga in the first layers of
a GaAs(110) surface and we combine experimental re-
sults with density functional theory (DFT) calculations.
To our knowledge, Cr atoms have yet to be explored by
STM on GaAs. It has been theoretically demonstrated
that they behave as Mn and give rise to ferromagnetic
coupling [5], but only a few DFT studies have focused
on the electronic properties of transition metal impuri-
ties near a semiconductor surface [19–21]. Experimen-
tally, while Cr and Mn induce different modifications in
the electronic structure (energy position of the acceptor
level, order of the resonance states), they both show ef-
fective magnetic coupling at low temperature [22]. Thus,
our study of Cr by STM/STS allows us to complete the
understanding of the mechanism behind the transition
metal doping of semiconductors.
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) and (b): Schematic three-dimensional view of the first three layers of the GaAs(110) surface with one Cr atom incorporated. As atoms are in red, Ga in gray and Cr in green. (c)-(f): Constant current STM image of one Cr atom substituted to a Ga in the first layer of GaAs(110) [(c) and (e)], and in the second layer [(d) and (f)]. Ball and stick atomic models are superposed to the STM images; large balls correspond to atoms in the first layer while smaller balls correspond to atoms in the second layer as defined in (a) and (b). Parameters: 3.5 nm × 3.5 nm, (c) –1.5 V, 50 pA, z range 39 pm, (d) –1.5 V, 100 pA, z range 60 pm, (e) 1.5 V, 50 pA, z range 360 pm, and (f) 1.5 V, 60 pA, z range 444 pm</figcaption></figure>

To control the substitution of Ga by Cr in the first and second layers of the surface, we developed a tip-induced manipulation method. This gives the possibility of studying magnetic impurities as a function of their coordination number: three covalent bonds with the As neighbors in the case of Cr inserted in the first layer, noted $Cr_{Ga1}$ and shown in Fig. 1(a), and four bonds in the case of Cr in the second layer, noted $Cr_{Ga2}$ and shown in Fig. 1(b), which has the same coordination number as in the bulk. For these two configurations, we image the wave function of the impurity and measure the spectra. Combining the experimental spatial and energy information with first principles total energy calculations, we obtain spatial maps of the different contributions, *s*, *p*, *d* of the wave function, which provides a detailed view of the hybridization between the *d*-orbital of the Cr atom and the *s*–*p* orbitals of the anion.

Experiments were performed with a STM system at 4.2 K in ultrahigh vacuum ($< 5 \times 10^{-11} \text{ mbar}$). The gap voltage was applied on the sample. We have cleaved *in situ* p-type highly doped ($2 \times 10^{19} \text{ cm}^{-3}$) GaAs wafers in order to expose the (110) surface and to perform so-called cross-sectional STM (X-STM) [23]. The sample

was pre-cooled in the STM cryostat and then transferred to a preparation chamber where a submonolayer quantity of Cr was deposited on the GaAs(110) surface, and replaced quickly in the STM for measurement at 4.2 K. As described by D. Kitchen et al. [10, 13, 24], transition metal adatoms deposited on GaAs(110) can be incorporated in the first surface layer. Mn, Fe and Co adatoms on GaAs(110) show different surface dynamic and incorporation behaviors. For Cr, we find that adatoms can be incorporated in the first GaAs layer using tip manipulation (Fig. 1(a)). Cr in the first layer can also be forced to move to a substitution position in the second layer (Fig. 1(b)). In each case of substitution, we observe the extraction of a Ga atom visible as an adatom: black rings also reported in Refs. [10] and [13]. The mechanism of tip-induced substitution and atoms manipulation will be described elsewhere [25].

## ELECTRONIC STRUCTURE OF Cr ATOMS IN Ga SUBSTITUTION POSITION

Figures 1(c)-1(f) show constant current STM images of Cr atoms in the Ga substitution position in the first and second layers of GaAs(110) for negative and positive tip-sample gap voltage. Because we control the incorporation with the tip, we are able to determine the exact Cr position before and after the incorporation processes. As expected for Ga substitution, Cr is localized in the Ga [1-10] columns (vertical lines in Fig. 1), which are imaged at positive voltage [26]. $Cr_{Ga1}$ is imaged on a Ga [001] row (horizontal lines) of the first layer (Fig. 1(e)) while $Cr_{Ga2}$ appears in between two Ga rows, i.e. in an As row of the first layer and Ga row of the second layer [Fig. 1(f)]. In other words, the Cr that substitutes a Ga of the first layer (respectively second layer) is symmetric with respect to the (1-10) plane crossing the Ga row (respectively As row) of the top surface layer. From these symmetry arguments, we can deduce unambiguously that the Cr is inserted in the second layer. In the first layer, we find similar patterns with drop-eyes (negative voltage) and crab-like (positive voltage) shape as found on other transition metal impurities [12, 13, 17, 18]; in the second layer the patterns change drastically. Such an observation of a transition metal atom in the second layer has been reported for Mn in Ref. [12], but no atomic resolution and no bias voltage dependence have been shown. It is noteworthy that in both layers the maximum of the electronic density of states is not localized exactly on the Cr atoms but distributed rather asymmetrically over several atomic lengths in a strongly anisotropic way.

Figure 2 presents differential conductance $dI(V)/dV$ spectra measured on Cr, proportional to the local electronic density of states (LDOS). They show features similar to those measured on Mn, Fe, and Co [13, 17, 18, 24]. Negative differential conductance is also observed due to
<figure><img src="image_3.png" /><figcaption>FIG. 2. STS spectra measured on top of Cr atoms incorporated (a) in the first layer of GaAs(110) surface and (b) in the second layer. Exact position of the tip during the spectra, named the center or lobe, are indicated by small squares depicted in Figs. 1(e) and (f). Peaks observed on the spectra are indicated and labeled on each peak.</figcaption></figure>

the presence of a narrow resonance peak [27]. Depending on the exact spectra position, center or lobe, large in-gap resonances are observed in the GaAs gap near the conduction band and smaller peaks in the valence band.

In order to precisely identify the spatial distribution of each state, we performed differential conductance maps on a 52×52 grid. A conductance map taken at the energy of an electronic resonance represents the probability density of the corresponding wave function. Figure 3 shows $dI(V)/dV$ maps of $Cr_{Ga1}$ (left) and $Cr_{Ga2}$ (right) at the energies of the peaks observed in Fig. 2 (except peak $p_{21}$ for which the signal is too low). For $Cr_{Ga1}$, the first state $p_{11}$ [Fig.3(a)] is weak and delocalized. In contrast, the second state $p_{12}$ (Fig.3(b)) is localized near the As neighbor atoms. The $p_{13}$ state (Fig. 3(c)) is similar to $p_{12}$ but present in addition a signal contribution in between the As atoms, closer to the Cr position. For $Cr_{Ga2}$, the states $p_{22}$ [Fig. 3(d)] and $p_{24}$ [Fig. 3(f)] present an isotropic spatial extension while $p_{23}$ [Fig. 3(e)] has an delocalized wave function. The anisotropic shape of the Cr atom in the constant current image of Fig. 1(f) can be well reproduced by the summation of the conductance images at the positive bias of Figs. 3(e) and 3(f).

<figure><img src="image_4.png" /><figcaption>FIG. 3. $dI(V)/dV$ maps of Cr in the Ga substitution position in the first surface layer (top row) and in the second layer (bottom row), 2 nm × 2 nm. Images are measured at the energies of the peaks observed in Fig. 2. Schematic ball and stick atomic model are superposed to STM images. Insets of each $dI(V)/dV$ map: calculated cross section on the (110) plan of ILDOS of Cr in the first surface layer (left) and in the second surface layer (right) at energy ranges corresponding to each peaks of Fig. 4, indicated in brackets. Maps are integrated over a slice along [110] from the surface to the plan of the $Cr_{Ga2}$ atom. Isodensity contours are marked by black lines and are separated by 0.005 e/bohr$^3$. Green dots indicate the position of the Cr atom. The whole color range on each ILDOS map is (in e/bohr$^3$ unit): 0.054 in (a), 0.011 in (b), 0.029 in (c, 0.026 in (d), 0.028 in (e), and 0.031 in (f).</figcaption></figure>

## COMPUTATIONAL DETAILS

To further explain and clarify the origin of the observed peaks in Fig. 2 and their spatial distribution in Fig. 3, we have carried out first principles total energy calculations in the framework of density functional theory (DFT) as implemented in the plane-wave PWSCF code of the QUANTUM ESPRESSO distribution [28]. A local spin density approximation (LSDA) with norm-conserving pseudopotentials with the Perdew-Zunger exchange-correlation functional were employed. The GaAs (110) surface was modeled through the supercell method using a slab containing 13 atomic planes each containing 2 atoms (26 atoms for the whole structure). As for the doped surface slabs, a 2×3 supercell was built
<figure><img src="image_5.png" /><figcaption>FIG. 4. Calculated spin resolved density of states projected on the atomic orbitals for the Cr-doped (110) GaAs surface with the Cr atom in the (a) first and (b) second layer of the surface. The insets show the total density of states and the vertical black bars show the position of the maxima. The zero of the energy is set to the Fermi energy of the system. Spin up correspond to upper curves, and spin down to the lower ones.</figcaption></figure>

containing 156 atoms in order to reduce the interaction between periodic replicas of dopants [see Supplemental Material [29] and reference therein [5, 28, 30–39]]. We simulated the effect of Cr dopants in both the bulk and (110) GaAs surface. For each system the projected density of states (PDOS) was calculated by projecting the Kohn-Sham wave functions on an atomic orbitals basis set in order to evaluate the contribution of each atom to the total DOS. Experimental constant-current STM images are here discussed in light of integrated local density of states (ILDOS) calculations in the framework of the Tersoff-Hamann approximation [38, 39]. According to this assumption, the STM tunneling current is proportional to the integrated LDOS of the sample at the tip position. By employing this approach, constant-current STM images can hence be directly compared to the integrated LDOS calculated within DFT. All the computational details of our study can be found in the Supplemental Material [29].

## *p*-d HYBRIDIZATION

We calculated the electronic density of states projected on atomic orbitals for Cr-doped bulk GaAs and for a Cr atom in the first, second and third layer of the GaAs surface. In the bulk case (see SM), the Cr site has $T_d$ symmetry and produces a deep trap state within the gap that results from the hybridization of the Cr *d* orbitals and the *p* orbitals of the As atoms (see Fig. S1 and Table S1 of the SM). The tetrahedral crystal field of the zinc-blend host splits the *d* levels of the Cr atoms into *e* and $t_2$ multiplets, which are further split by the exchange interaction into spin up and spin down states. Hybridization finally splits the $t_2$ states into bonding and antibonding levels (*e* states have no orbital in the GaAs of the same symmetry to hybridize with), giving rise to six states in total (Figs. S2 and S7 of the SM). According to Ref. [5], this leads to two kinds of impurity states: the crystal field resonances (CFRs) that have *d*-orbital character and are localized on the transition metal ion and the dangling bond hybrid (DBH) that contains a large *p*-orbital component due to the *p*- *d* hybridization.

When the Cr atom substitutes a Ga atom in the third layer of the surface (Fig. S6), the electronic structure does not differ too much from what we observe in the case of Cr-doping of the bulk system. As one goes from the third to the second and first layer of the surface, the symmetry of the Cr site is lowered as evidenced by the evolution of first neighbours distances around the Cr impurity (see Table S2). For the $Cr_{Ga1}$ substitution in the first layer, the three-fold coordination and the surface buckling leads to $C_{3v}$ point-group site-symmetry. Furthermore, because of the change of the lattice parameters at the surface, a $C_S$ site-symmetry is finally obtained. As a consequence, the five-fold degeneracy of the *d*-orbital multiplet is completely lifted as shown in Fig. S7.

The insets of Fig. 4 show the total density of states for the $Cr_{Ga1}$ and $Cr_{Ga2}$ that exhibit three peaks for $Cr_{Ga1}$ and four for $Cr_{Ga2}$, corresponding to what is observed experimentally in Fig. 2. Figure 4 shows the contributions of atomic orbitals to the total density of states (*d* and *p* orbitals, while *s* orbitals have a negligible contribution), for the two spin polarizations. All the peaks observed have a strong *d* contribution coming from the chromium, as well as a an As-*p* contribution. For $Cr_{Ga1}$ (respectively $Cr_{Ga2}$), this *p* contribution is especially strong for the peaks at −0.3 eV (−0.6 eV) and 0.7 eV (0.5 eV). For the peak at 0.7 eV for $Cr_{Ga1}$ and the peaks at −0.1 eV and 0.8 eV for $Cr_{Ga2}$, the contribution of Ga-*p* orbitals is also significant as compared to the *d* contribution. This mixed character of the peaks is the consequence of *p*- *d* hybridization.

The insets of Fig. 3 show ILDOS of $Cr_{Ga1}$ (left) and $Cr_{Ga2}$ (right) at energies corresponding to the peaks of Fig. 4 and associated with the experimental maps. Here,
we represent (110) planes of ILDOS integrated over a slice along [110] from the surface to the plane of the $Cr_{Ga2}$ atom. These calculations reproduce well the $dI/dV(V)$ experimental images, considering the well-known DFT inaccuracy in reproducing experimental energy peak positions, except for in Fig. 3 (a) where the signal is too weak. For $Cr_{Ga1}$, the observed and simulated lobes are oriented towards the As first neighbors of the Cr, which can be explained by the strong As-p contribution to the density of states found in Fig. 4. Thus, the DFT results show that our STM observation of the delocalization and the anisotropy of the wave function, especially the density of states oriented towards the As neighbors giving rise to the two drop-eye lobes, is a direct demonstration of p-d hybridization observed here at the atomic scale.

Some states present wave functions with a rather compact shape centered around the Cr [insets of Fig. 3 (a), (d) and (f)], while others are more delocalized with anisotropic lobes [insets of Fig. 3 (b), (c) and (e)]. States centered on the Cr impurities have a CFR character, with a rather out-of-plane contribution, which could be described by $d_{z^2}$ orbitals. In contrast, states with lobes pointing towards the anions have a DBH character with the symmetry of $d_{xy}$ orbitals. One expects in this case a naturally strong $\sigma$ bonding with the p orbitals of the As neighbors which have the same symmetry.

Because of the lowering of site symmetry from $T_d$ in the bulk to $C_S$ at the surface, the energy levels identified by STM spectroscopy and reproduced by the DFT calculations cannot be classified according to the irreducible representations $t_2$ and $e$ of the $T_d$ point-group. Despite this, one can still identify two types of states, the localized states, essentially of d-orbital character, and the extended hybridized p-d states. Because these extended states spread over several unit cells, they should be better described as the hybridization of d-orbitals with valence Bloch states at the $\Gamma$ point [40, 41].

CONCLUSION

In conclusion, by combining STM and first principles calculations, we were able to explain the origin of the anisotropic shape of the wave function of the Cr atom substituted to Ga in the GaAs(110) surface. The lobes observed by STM in GaAs doped by transition metals are the signature of the electronic hybridization of the atom with its host. In particular, the As-p orbitals hybridize with the Cr-d orbitals giving rise to "drop-eye" lobes localized around the As atoms, which are the main distinctive attributes of the wave function of transition metal impurities in GaAs. Our DFT calculations describe the electronic origin of the resonances observed in the semiconductor gap when doped with a transition metal, and show the significant weight of the As-p, as well as the Ga-p electrons in some cases. We also present the spin

resolved contribution of these states, whose experimental confirmation could be an interesting perspective of this work.

The work was supported by the French ANR project MechaSpin (ANR-17-CE24-0024). The high-performance computing (HPC) resources for this project were granted by the Institut du developpement et des ressources en informatique scientifique (IDRIS) under the allocations A0040907723 and c2016097723 via GENCI (Grand Equipement National de Calcul Intensif).

[1] B. Clerjaud, Journal of Physics C: Solid State Physics **18**, 3615 (1985).

[2] C. A. Bates and K. W. H. Stevens, Reports on Progress in Physics **49**, 783 (1986).

[3] T. Dietl, H. Ohno, F. Matsukura, J. Cibert, and D. Ferland, Science **287**, 1019 (2000).

[4] T. Dietl, H. Ohno, and F. Matsukura, Phys. Rev. B **63**, 195205 (2001).

[5] P. Mahadevan and A. Zunger, PRB **69**, 115211 (2004).

[6] S. Sanvito, P. Ordejón, and N. A. Hill, Phys. Rev. B **63**, 165206 (2001).

[7] T. C. Schulthess, W. M. Temmerman, Z. Szotek, W. H. Butler, and G. Malcolm Stocks, Nature Materials **4**, 838 (2005).

[8] T. Tsuruoka, R. Tanimoto, N. Tachikawa, S. Ushioda, F. Matsukura, and H. Ohno, Solid State Commun. **121**, 79 (2002).

[9] A. M. Yakunin, A. Y. Silov, P. M. Koenraad, J. H. Wolter, W. Van Roy, J. De Boeck, J.-M. Tang, and M. E. Flatté, Phys. Rev. Lett. **92**, 216806 (2004).

[10] D. Kitchen, A. Richardella, J.-M. Tang, M. E. Flatte, and A. Yazdani, Nature **442**, 436 (2006).

[11] J.-M. Jancu, J.-C. Girard, M. O. Nestoklon, A. Lemaitre, F. Glas, Z. Z. Wang, and P. Voisin, Phys. Rev. Lett. **101**, 196801 (2008).

[12] J. K. Garleff, C. Celebi, W. Van Roy, J.-M. Tang, M. E. Flatte, and P. M. Koenraad, Phys. Rev. B **78**, 075313 (2008).

[13] A. Richardella, D. Kitchen, and A. Yazdani, Physical Review B **80**, 045318 (2009).

[14] C. Celebi, P. M. Koenraad, A. Y. Silov, W. Van Roy, A. M. Monakhov, J.-M. Tang, and M. E. Flatte, Phys. Rev. B **77**, 075328 (2008).

[15] A. Taninaka, S. Yoshida, K. Kanazawa, E. Hayaki, O. Takeuchi, and H. Shigekawa, Nanoscale **8**, 12118 (2016).

[16] J. Bocquel, V. R. Kortan, C. Şahin, R. P. Campion, B. L. Gallagher, M. E. Flatté, and P. M. Koenraad, Phys. Rev. B **87**, 075421 (2013).

[17] S. Mühlenberend, M. Gruyters, and R. Berndt, Physical Review B **88**, 115301 (2013).

[18] A. L. Benjamin, D. Lee, and J. A. Gupta, Journal of Vacuum Science and Technology B, Nanotechnology and Microelectronics: Materials, Processing, Measurement, and Phenomena **31**, 04D102 (2013).

[19] Z. Fang and Z. Yi, Computational Materials Science **118**, 4 (2016).
[20] Z. Fang and Z. Yi, Physica B: Condensed Matter **502**, 1 (2016).

[21] M. F. Islam and C. M. Canali, Physical Review B **85**, 155306 (2012).

[22] H. Wu, H. D. Gan, H. Z. Zheng, J. Lu, H. Zhu, Y. Ji, G. R. Li, and J. H. Zhao, Solid State Communications **151**, 456 (2011).

[23] J.-C. Girard, A. Lemaitre, A. Miard, C. David, and Z. Z. Wang, J. Vac. Sci. Technol. B **27**, 891 (2009).

[24] D. Kitchen, *Atom-by-Atom Substitution of Transition Metals in GaAs and Visualization of Hole-Mediated Interactions*, Ph.D. thesis, University of Illinois (2006).

[25] G. Rodary, to be published.

[26] R. M. Feenstra, J. A. Stroscio, J. Tersoff, and A. P. Fein, Phys. Rev. Lett. **58**, 1192 (1987).

[27] Y. Xue, S. Datta, S. Hong, R. Reifenberger, J. I. Henderson, and C. P. Kubiak, PRB **59**, R7852 (1999).

[28] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococioni, I. Dabo, A. D. 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).

[29] See Supplemental Material at [URL] for computational details.

[30] X. Duan, M. Peressi, and S. Baroni, Phys. Rev. B **72**, 085341 (2005).

[31] X. Duan, S. Baroni, S. Modesti, and M. Peressi, Appl. Phys. Lett. **88**, 022115 (2006).

[32] X. Duan, M. Peressi, and S. Baroni, Phys. Rev. B **75**, 035338 (2007).

[33] J. L. A. Alves, J. Hebenstreit, and M. Scheffler, Phys. Rev. B **44**, 6188 (1991).

[34] Z. Yi, Y. Ma, and M. Rohlfing, J. Phys. Chem. C **115**, 23455 (2011).

[35] F. Tilley, M. Roy, P. Maksym, P. Koenraad, C. Krammel, and J. Ulloa, Phys. Rev. B **93**, 035313 (2016).

[36] S. Tong, A. Lubinsky, B. Mrstik, and M. Van Hove, Physical Review B **17**, 3303 (1978).

[37] J. Thakur, H. S. Saini, M. Singh, A. Reshak, and M. K. Kashyap, Physica E: Low-dimensional Systems and Nanostructures **78**, 35 (2016).

[38] J. Tersoff and D. Hamann, Phys. Rev. B **31**, 805 (1985).

[39] F. Giustino, *Materials modelling using density functional theory: properties and predictions* (Oxford University Press, 2014).

[40] B. Koiller, X. Hu, and S. D. Sarma, Physical Review Letters **88**, 027903 (2001).

[41] B. Voisin, J. Bocquel, A. Tankasala, M. Usman, J. Salfi, R. Rahman, M. Y. Simmons, L. C. L. Hollenberg, and S. Rogge, Nature Communications **11**, 10.1038/s41467-020-19835-1 (2020).