Spin splitting in 2D monochalcogenide semiconductors

We report ab initio calculations of the spin splitting of the uppermost valence band (UVB) and the lowermost conduction band (LCB) in bulk and atomically thin GaS, GaSe, GaTe, and InSe. These layered monochalcogenides appear in four major polytypes depending on the stacking order, except for the monoclinic GaTe. Bulk and few-layer ε-and γ -type, and odd-number β-type GaS, GaSe, and InSe crystals are noncentrosymmetric. The spin splittings of the UVB and the LCB near the Γ-point in the Brillouin zone are finite, but still smaller than those in a zinc-blende semiconductor such as GaAs. On the other hand, the spin splitting is zero in centrosymmetric bulk and even-number few-layer β-type GaS, GaSe, and InSe, owing to the constraint of spatial inversion symmetry. By contrast, GaTe exhibits zero spin splitting because it is centrosymmetric down to a single layer. In these monochalcogenide semiconductors, the separation of the non-degenerate conduction and valence bands from adjacent bands results in the suppression of Elliot-Yafet spin relaxation mechanism. Therefore, the electron- and hole-spin relaxation times in these systems with zero or minimal spin splittings are expected to exceed those in GaAs when the D’yakonov-Perel’ spin relaxation mechanism is also suppressed.

Potential applications in spin-dependent electronics and optoelectronics have driven the search for materials capable of exhibiting a high degree of spin polarization and long spin relaxation time 1,2 . However, optical generation of electron and hole spin polarization and resulting polarized luminescence are typically limited by the mixing of degenerate valence bands in most semiconductors 2 . Recent reports of valley polarization in atomically thin transition metal dichalcogenides (TMDs) [3][4][5][6][7] suggest potential exploitation of both spin and valley degrees of freedom for electronics and optoelectronics. In an experimental study 8 , we demonstrated the high generation and preservation of optical spin polarization and dynamics in a group-III monochalcogenide, GaSe, under nonresonant optical pumping. The observed near unity optical spin polarization 9,10 is attributed to suppressed electron and hole spin relaxation rates resulting from reduced valence-band mixing. However, the microscopic spin relaxation mechanisms in GaSe and related monochalcogenides are not fully understood.
In metals and semiconductors, the major spin relaxation mechanisms-including Elliott-Yafet (EY) 11,12 and D'yakonov-Perel' (DP) [13][14][15] mechanisms-are associated with the spin-orbit interaction (SOI) and the spin-orbit-induced spin splitting, ∆ ( ) = | ( , ↑) − ( , ↓)| 16 . Considering spin-relaxation with a four-state (two bands with spin) model Hamiltonian in the absence of an external magnetic field, one can relate the spin relaxation rate of electrons (holes) when the Fermi energy (corresponding to Fermi vector k F ) is away from the conduction (valence) band edge with the following equation 16  where τ Γ = / ħ between the adjacent bands with energy separation Δ g . In GaSe, the p z -like uppermost valence band (UVB) is well isolated from the lowermost conduction band (LCB) (~2 eV) and the adjacent , p x y -like valence bands, and as a result L/Δ g ≈ 0.02-0.04 17,18 . The hole-spin relaxation due to the EY mechanism is thus expected to be much smaller than the momentum relaxation rate Γ p . The spin relaxation caused by the DP mechanism can be seen as being due to the precession of spins in an effective magnetic field associated with Δ s (k) 2,15,16 . The DP spin relaxation rate is proportional to the spin splitting, where τ p is the momentum relaxation time. Therefore, when the spin relaxation is dominated by the DP mechanism, the smaller the spin splitting, the longer the spin relaxation time τ = /Γ ħ s s for the same momentum relaxation rate Γ p . To understand the spin relaxation, one first needs the momentum (  k)-dependent ∆ ( )  k s of the bands near the fundamental gap. In the absence of magnetic fields, ∆ ( )  k s is zero in centrosymmetric crystals because of the constraints of time-reversal symmetry [ ( ,↑) = (− ,↓)   E k E k , Kramers degeneracy] and spatial inversion symmetry When the inversion symmetry is broken in crystals (bulk inversion asymmetry (BIA)) 19 or heterostructures (structural inversion asymmetry (SIA)) 20 in GaAs and other zinc-blende semiconductors has been a subject of considerable interest [22][23][24][25][26][27] since the seminal work of Dresselhaus 19 . Ab initio calculations, such as LDA (or GGA) and self-consistent GW methods, of ∆ ( )  k s in bulk GaAs and two-dimensional GaAs-based superlattices and heterostructures have improved the understanding of the spin splitting [23][24][25][26] . A few theoretical calculations of ∆ ( )  k s in TMDs also have been reported 28,29 . In this study, we report ab initio calculations of ∆ ( )  k s of the uppermost valence band (UVB) and the lowermost conduction band (LCB) in GaSe and related group-III monochalcogenides, including GaS, GaTe, and InSe. 3 ) are noncentrosymmetric. Bulk ε-, γ-, and δ-MX crystals, which appear in an AB, ABC, and ABCD stacking order, are noncentrosymmetric, while β-MX crystals are centrosymmetric with an AB stacking order. Additionally, there are two exceptions: (1) an atomically thin β-MX crystal with even-number layers is centrosymmetric, and (2) a bilayer δ-MX crystal can be identical to either a bilayer ε-MX crystal (noncentrosymmetric) or a bilayer β-MX crystal (centrosymmetric) depending on which two layers are isolated from a bulk δ-MX crystal. GaTe appears as a distorted form of the MX structure, where one out of three Ga-Ga bonds lies in the a − b plane (Fig. 1c,d). In contrast to MX crystals, GaTe crystals belong to the monoclinic lattice system (space group C h 2 3 ), and are centrosymmetric down to a single layer 30 . Bulk and few-layer εand γ-type, as well as odd-number few-layer β-type GaS, GaSe, and InSe crystals exhibit finite spin splittings, while bulk and even-number few-layer β-type GaS, GaSe, and InSe as well as GaTe crystals exhibit zero spin splitting. The difference is due to the constraints of the aforementioned time-reversal and spatial inversion symmetry (or the lack of it).

Results
Band structure: bulk versus single-layer. In Fig. 2, we show the electronic band structures of bulk and monolayer β-GaS, ε-GaSe, and ε-InSe, which are the most naturally abundant. The general features of the electronic band structures, except the spin splitting, are nearly polytype-independent, owing to the weak inter-layer interactions. The lowermost conduction band (LCB) has s-like symmetry, whereas the two uppermost valence bands (UVBs) have p z -like symmetry. The , p x y -like valence bands appear ~1 eV below the UVB as a result of the crystal field and SOI. The calculated band structures for ε-GaSe show a nearly direct band gap at the Γ-point of the Brillouin zone (BZ), where a valley appears in the UVB. The energy of the LCB at the Γ-point is ~0.5 eV lower than that at the M point, consistent with the hybrid density functional calculations 31 . On the contrary, tight binding calculations show that the energy of the LCB at the M point for GaSe is ~10 meV below that at the Γ-point in the BZ 32 .
The band gap is seen to decrease with increasing atomic number (Ga → In or S → Se). The calculated band gaps are 2.0 eV, 1.3 eV, and 0.71 eV for β-GaS, ε-GaSe, and ε-InSe, respectively, which are each smaller than the experimental values (~3.1 eV, 2.0 eV, and 1.3 eV) 30 . The band-gap underestimation can be remedied with, for example, the HSE06 hybrid functional 31,33,34 . In the absence of SOI, the , p x y states are doubly degenerate at the Γ-point. On the other hand, the SOI lifts this energy degeneracy with a spin-orbit splitting Δ SO ≈ 0.09 eV, 0.34 eV, and 0.31 eV in GaS, GaSe, and InSe, respectively. Δ SO in GaSe and InSe are similar in magnitude, but a factor of three smaller in GaS, agreeing with previously reported calculations [35][36][37][38] . Δ SO in GaS is minimal, as expected from the weak SOI in the lighter S anions which govern the characteristics of the few uppermost valence bands of GaS. Monolayer GaS, GaSe, and InSe have very similar band structures (Fig. 2). We note two different features in the band structures of monolayer MXs in comparison with their bulk counterparts: (1) the quantum confinement along the c-axis increases the band gap to 2.36 eV, 1.78 eV, and 1.4 eV for GaS, GaSe and InSe, respectively, and (2) the band gap becomes indirect as the valley at the Γ -point becomes wider in momentum ( =  k k ) and deeper in energy (E).
The band structure of GaTe (bulk) has also been calculated with GGA 39,40 , showing a direct band gap of ~1 eV. The inclusion of SOI causes negligible changes in the UVB and LCB of GaTe. Monolayer GaTe shows a direct band gap of 1.4 eV (Fig. 3), with LCB having two nearly degenerate minima at the Γ and C points. At the C point, LCB has s-like symmetry while the UVB has p y -like symmetry. SOI removes the , p x y degeneracy of valence bands at Γ, with Δ SO ≈ 0.2 eV. Δ SO is smaller in GaTe than in GaSe despite Te being heavier than Se. The reduction in the strength of Δ SO is due to the quenching of orbital angular momentum in the lower symmetry crystalline structure, as demonstrated by a sizable Δ SO ≈ 0.7 eV calculated for a hypothetical β-type GaTe (space group D h  In contrast to ε-GaSe, bulk and even-number few-layer β-GaSe crystals have zero spin splitting (Fig. 4c,d), obeying the constraint of spatial inversion symmetry. The ∆ ( ′) k s v and ∆ ( ′) k s c in odd-number few-layer β-GaSe crystals are finite, but diminish rapidly with increasing layers. In trilayer β-GaSe, the UVB spin splitting is less than 1 meV, and LCB spin splitting is smaller by a factor of five compared to that of the monolayer. The thickness dependent spin splitting in β-GaSe presented here are consistent with those reported in MoS 2 29 , which has the same symmetry as β-GaSe. Bulk γ-GaSe has similar spin splittings as bulk ε-GaSe, with decreasing spin splittings as the number of layers increase.
In Fig. 5, we compare ∆ ( ′) k s v and ∆ ( ′) k s c in monolayer GaS, GaSe, and InSe (group-III monochalcogenides) and bulk GaAs (a representative zinc-blende III-V semiconductor). Among the monolayer group-III monochalcogenides, overall spin splittings decrease from GaSe, to InSe, and then to GaS. The spin splittings typically increase with the increasing atomic number of constituent atoms as result of the    enhanced SOI in the heavier atoms. However, other details of the band structure such as the band gap also contribute to the spin splittings.
The valence band of GaAs consists of a heavy hole (HH), a light hole (LH), and a split-off (SO) band 22,27 . The calculated HH spin splitting is close to that in the UVB of GaSe. However, the spin splittings in the LH and SO bands are at least a factor of two larger than that in the UVB of GaSe. The calculated overall LCB spin splitting in GaAs is also larger than that in GaSe. The magnitude at k′ = 0.15 is more than two times larger in GaAs than in GaSe. The spin splitting of the heavy-hole band is reduced by about a factor of two when the GW method is used in lieu of the GGA method.

Discussion
The spin splittings discussed above concern mainly the overall spin splitting up to k′ = 0.15. To understand the spin relaxation mechanisms, we need to identify the k-dependence of the spin splitting in the vicinity of the Γ-point. At small k, the ⋅ k p theory predicts that, in noncentrosymmetric zinc-blende and wurtzite structures, the k-dependence of the spin splitting contains both a linear and a cubic term when the core levels are considered 19,26,41,42 . To illustrate the k-dependence of the spin splitting, we fit the calculated ∆ ( ′) k s in ε-GaSe with the function ∆ ( ′) = ′ + ′ k A k B k s 3 for k′ < 0.05 (Table 1). The energy scales for the coefficients A (meV) and B (eV) are consistent with those determined from GW calculations for GaAs 23,26 . Although B is three to four orders of magnitude larger than A, there exists a crossover value of ′ = / ∼ − k A B 10 2 below which the linear term dominates. In contrast to the GaAs case where the linear term is negligible for the LCB, we find a sizable linear term for the LCB in GaSe. The cubic coefficient B for the LCB is ~4-5 eV for the monolayer to the bulk, with the bilayer case being slightly different. In contrast to the UVB, there appears to be an odd-even-layer effect: B values for the odd layers (1 and 3) are larger, but are close to the bulk values for the even layers (2 and 4). For the LCB, the A values are similar for all the layers, except for the bilayer (A = 2.0 meV) and the bulk (A = 0.2 meV). Note that bilayer GaSe has an unusually large A value, and for the UVB, there appears to be an odd-even effect like that in the LCB. The A value for the bilayer is nearly three times that for the monolayer, whereas A for the four-layer is two times that for the trilayer. As pointed out in the case of GaAs, these subtle differences are due to the characteristics of the UVB and LCB energy values and wave functions, and their mixing with other bands including the core levels 19,26,41 .
The DFT-based theories such as LDA (or GGA) underestimate band gaps and do not give accurate effective masses, resulting in overestimated ∆ ( )  k s 23,25,26 . GW calculations reproduce more accurate band parameters, such as the band gap and effective mass, but are computationally more intensive than LDA (GGA) calculations. For simplicity, in this work, we have used GGA to calculate ∆ ( )  k s . The GGA calculation underestimates the GaAs band gap by a factor of ten; however, the spin splitting only deviates from that determined by the GW calculation by a factor of two. The band gaps are underestimated by the GGA for GaSe and related monochalcogenides by a factor of approximately two, which is significantly less than that for GaAs. For example, in GaSe, the GGA calculation gives a band gap of about 1 eV, which is off from the GW/HSE06 calculation 34 and the measured band gap (~2 eV) 35 by a factor of two. Therefore, we expect the GGA calculation to produce spin splittings close to the value obtained with the GW calculation. We also expect similar variations of ∆ s with k from one conduction/valence band to another and from bulk to atomically thin layers.

Conclusion
We present a systematic study of spin-orbit-induced spin splittings bulk and atomically thin group-III monochalcogenides MX′ (M = Ga, In; X′ = S, Se, Te). The spin splitting vary with anion element and crystal symmetry. Centrosymmetric crystals, including bulk β-type GaS, GaSe, and InSe, as well as monoclinic GaTe down to the monolayer, have zero spin splitting, as anticipated from the constraints of spatial inversion symmetry and time-reversal symmetry. Among the monolayer group-III monochalcogenides, overall spin splittings decrease from GaSe, to InSe, and then to GaS. The calculated spin splitting in the UVB of GaSe is close to that of the HH, but is at least a factor of two smaller than those in the LH and SO bands in GaAs. The calculated overall LCB spin splitting in GaSe is also smaller than that in GaAs. The magnitude at k′ = 0.15 is more than two times smaller in GaSe than in GaAs. In these monochalcogenide semiconductors, the separation of the non-degenerate conduction and valence bands from other adjacent bands results in suppression of Elliot-Yafet spin relaxation mechanism. Therefore, the electron and hole spin relaxation times in these systems with zero or minimal spin splittings and reduced valence-band mixing are expected to be longer than those in a zinc-blende semiconductor (eg., GaAs 22,27,43 ), owing to the suppression of D'yakonov-Perel' and Elliot-Yafet spin relaxation mechanisms.

Methods
We compute the band structures and ∆ ( )  k s of valence and conduction bands with the projector augmented wave method as implemented in the VASP [44][45][46][47][48] package and the full-potential (linearized) augmented plane-wave as implemented in the WIEN2k 49,50 package. The band structures are calculated with the WIEN2k package, with the optimized crystal structures determined by minimizing the total energy with all electrons (including core electrons) with VASP. In all calculations, exchange-correlation energies are determined by the Perdew-Burke-Ernzerhof (PBE) 51 generalized gradient approximation (GGA) 52 , which systematically underestimates the band gaps and produces ( )  E k dispersions (effective masses) different from experimental values. These shortcomings of the GGA also limit the accuracy of the calculated ∆ ( )  k s . The spin-orbit interaction (SOI) is included in our calculations of the overall band structure and the spin splitting of a given band in a self-consistent manner using a second variation approach [53][54][55][56][57][58]  , where m e is the electron mass, c the speed of light, L and  S the orbital and spin momentum vectors, and ( ) V r an effective single particle local potential seen by the electron. This form of H so is correct as long as ( ) V r is local and isotropic. In Hartree approximation and LDA, the effective potential is indeed local, though it is not always isotropic. The isotropic approximation is valid because the dominant contribution to H so is from regions near the nucleus. However, local approximations do not give correct band structure near the band gap. The accuracy of the band gap can be improved with hybrid models such as the Heyd-Scuseria-Ernzerhof (HSE06) 33 (a mixture of non-local and local exchange) or GW-like theories.
To model a few-layer thin film, we create a supercell (supercell method 52 ) containing one few-layer structure and a 15-25 Å thick vacuum spacer, which is large enough to suppress interactions arising from the artificial periodicity present in the supercell method. The crystalline c-axis of the supercell is set perpendicular to the crystalline a-b plane. In this way, one can distinguish the effects of intra-and inter-layer interactions on the electronic structures in few-layer structures. The number of atoms in a unit cell is as follows: eight for εand β-MX, twelve for γ-MX, sixteen for δ-MX, and twelve for monoclinic GaTe. To obtain an energy accuracy of 0.1 meV in self-consistent calculations, we use Γ-centered Monkhorst-Pack 59  k-meshes of 24 × 24 × 4 and 24 × 24 × 1 for bulk and few-layer GaSe-type structures, respectively. For GaTe, we use meshes of 16 × 6 × 8 and 16 × 6 × 2 for bulk and few-layer GaTe, respectively.