Spin-orbit coupling in elemental two-dimensional materials

TL;DR

This study investigates spin-orbit coupling and spin mixing in various two-dimensional elemental materials using first-principles calculations, revealing how these effects depend on atomic number and can be tuned by Fermi level adjustments.

Contribution

It provides detailed first-principles analysis of spin-orbit effects in 2D materials, including spin mixing, anisotropy, and their dependence on atomic number and symmetry.

Findings

Spin-orbit coupling scales as Z^2, with spin admixture b^2 scaling as Z^4.

Non-zero spin mixing exists in graphene despite mirror symmetry.

Spin-mixing anisotropy can be controlled by Fermi level tuning.

Abstract

The fundamental spin-orbit coupling and spin mixing in graphene and rippled honeycomb lattice materials silicene, germanene, stanene, blue phosphorene, arsenene, antimonene, and bismuthene is investigated from first principles. The intrinsic spin-orbit coupling in graphene is revisited using multi-band $k\cdot p$ theory, showing the presence of non-zero spin mixing in graphene despite the mirror symmetry. However, the spin mixing itself does not lead to the the Elliott-Yafet spin relaxation mechanism, unless the mirror symmetry is broken by external factors. For other aforementioned elemental materials we present the spin-orbit splittings at relevant symmetry points, as well as the spin admixture $b^2$ as a function of energy close to the band extrema or Fermi levels. We find that spin-orbit coupling scales as the square of the atomic number Z, as expected for valence electrons in atoms. For isolated bands, it is found that $b^2\sim Z^4$. The spin-mixing parameter also exhibits giant anisotropy which, to a large extent, can be controlled by tuning the Fermi level. Our results for $b^2$ can be directly transferred to spin relaxation time due to the Elliott-Yafet mechanism, and therefore provide an estimate of the upper limit for spin lifetimes in materials with space inversion center.

Tables

Table 1. Table 1: Calculated lattice parameters a = | v → 1 | = | v → 2 | 𝑎 subscript → 𝑣 1 subscript → 𝑣 2 a=|\vec{v}_{1}|=|\vec{v}_{2}| , δ z subscript 𝛿 𝑧 \delta_{z} [see Fig. 1 ], and kinetic energy cut-offs for the wave function ( E cut Ψ superscript subscript 𝐸 cut Ψ E_{\text{cut}}^{\Psi} ) and charge density ( E cut ρ superscript subscript 𝐸 cut 𝜌 E_{\text{cut}}^{\rho} ) applied for structural optimization. The PBEsol exhange–corellation potential and 12 × \times 12 k 𝑘 k -point grid were assumed.

Material	$a$ [Å]	${\delta }_{\text{z}}$ [Å]	$E_{\text{cut}}^{\Psi }$  [Ry]	$E_{\text{cut}}^{\rho }$ [Ry]
graphene	2.459	0	58	696
silicene	3.84	0.45	58	580
germanene	3.99	0.66	38	380
stanene	4.6	0.84	48	432
blue. phosp	3.24	1.24	58	580
arsenene	3.61	1.38	45	450
antimonene	4.12	1.634	40	480
bismuthene	4.29	1.73	42	429

Table 2. Table 2: Orbital ( Δ Δ \Delta ) and spin–orbital gaps ( Δ so K superscript subscript Δ so K \Delta_{\text{so}}^{\text{K}} ), ( Δ so Γ superscript subscript Δ so Γ \Delta_{\text{so}}^{\Gamma} ) calculated in Wien2K for PBEsol exchange–correlation potential. The character of the orbital gap, direct or indirect, is labelled by capital the letter D or I respectively. Inverted band gaps are indicated by Inv.

Material	$\Delta$ [meV]	${\Delta }_{\text{so}}^{\text{K}}$ [meV]	${\Delta }_{\text{so}}^{\Gamma }$ [meV]
graphene	0.024 (D)	0.024	9
silicene	1.48 (D)	1.48	34.6
germanene	23 (D)	23	201
stanene	25 (I)	72	461
blue phosphorus	1864 (I)	10	48
arsenene	1492 (I)	71	195
antimonene	982 (I)	174	340
bismuthene	491 (I, Inv)	702	712

Equations

$${\mathbf{H}}_{\mathbf{SO}}$$

$$\{\begin{array}{l}⟨{\Gamma }_5^1\left|H_{SOz}\right|{\Gamma }_5^2⟩=i{\Delta }_5\\ ⟨{\Gamma }_{6v}^1\left|H_{SOy}\right|{\Gamma }_{6v}^2⟩=i{\Delta }_{6v}\\ ⟨{\Gamma }_5^1\left|H_{SOx}\right|{\Gamma }_{1v}⟩=⟨{\Gamma }_5^2\left|H_{SOy}\right|{\Gamma }_{1v}⟩={\Delta }_{51}\\ ⟨{\Gamma }_5^1\left|H_{SOx}\right|{\Gamma }_{6v}^1⟩=-⟨{\Gamma }_5^2\left|H_{SOx}\right|{\Gamma }_{6v}^2⟩=-⟨{\Gamma }_5^1\left|H_{SOy}\right|{\Gamma }_{6v}^2⟩=-⟨{\Gamma }_5^2\left|H_{SOy}\right|{\Gamma }_{6v}^1⟩={\Delta }_{56},\end{array}$$

$$\left\{\left[|{\Gamma }_5^{-}\uparrow ⟩,|{\Gamma }_{1v}\downarrow ⟩\right],\left[|{\Gamma }_5^{+}\downarrow ⟩,|{\Gamma }_{1v}\uparrow ⟩\right],\left[|{\Gamma }_5^{+}\uparrow ⟩,|{\Gamma }_{6v}^{-}\downarrow ⟩\right],\left[|{\Gamma }_5^{-}\downarrow ⟩,|{\Gamma }_{6v}^{+}\uparrow ⟩\right],\left[|{\Gamma }_{6v}^{-}\uparrow ⟩,|{\Gamma }_{6v}^{+}\downarrow ⟩\right]\right\},$$

$$\left[\begin{array}[]{cc|cc|cc|cc|cc}\Delta_{5}&\sqrt{2}\Delta_{51}&0&0&0&0&0&0&0&0\\ \sqrt{2}\Delta_{51}^{*}&E_{1v}&0&0&0&0&0&0&0&0\\ \hline\cr 0&0&\Delta_{5}&\sqrt{2}\Delta_{51}&0&0&0&0&0&0\\ 0&0&\sqrt{2}\Delta_{51}^{*}&E_{1v}&0&0&0&0&0&0\\ \hline\cr 0&0&0&0&-\Delta_{5}&2\Delta_{56}&0&0&0&0\\ 0&0&0&0&2\Delta_{56}^{*}&E_{6v}-\Delta_{6}&0&0&0&0\\ \hline\cr 0&0&0&0&0&0&-\Delta_{5}&2\Delta_{56}&0&0\\ 0&0&0&0&0&0&2\Delta_{56}^{*}&E_{6v}-\Delta_{6}&0&0\\ \hline\cr 0&0&0&0&0&0&0&0&E_{6v}+\Delta_{6}&0\\ 0&0&0&0&0&0&0&0&0&E_{6v}+\Delta_{6}\end{array}\right]\,,$$

$$E_c$$

$$E_v$$

$$|\psi \Uparrow {⟩}_{\text{c}}$$

$$|\psi \Downarrow {⟩}_{\text{c}}$$

$$|\psi \Uparrow {⟩}_{\text{v}}$$

$$|\psi \Downarrow {⟩}_{\text{v}}$$

$$\alpha$$

$$\beta$$

$$\gamma$$

$$\lambda$$

$$\eta$$

$$\nu$$

$${\Psi }_{n,\mathbf{k}}^{\Uparrow }(\mathbf{r})$$

$${\Psi }_{n,\mathbf{k}}^{\Downarrow }(\mathbf{r})$$

$$b_{\hat{\mathbf{s}}}^2=\frac{1}{\rho (E_F)S_{BZ}}{\int }_{FC}\frac{b_{\mathbf{k}}^2(\hat{\mathbf{s}})}{\hslash |{\mathbf{v}}_F(\mathbf{k})|}dk,$$

Full text

spin-orbit coupling, spin relaxation, spintronics, graphene, 2D materials

I Introduction

Atomically thin two dimensional (2D) materials have been attracting attention of physicists for over a decade. Many materials have been successfully synthesized Feng et al. (2012); Vogt et al. (2012); Meng et al. (2013); Dávila et al. (2014); Li et al. (2014); Derivaz et al. (2015); Zhu et al. (2015a); Zhang et al. (2016); Ji et al. (2016); Reis et al. (2017); Shah et al. ; Xia et al. (2014); Castellanos-Gomez et al. (2014); Shao et al. (2018); Zhang et al. (2018) opening new routes towards novel nano-electronic and spintronic devices. Graphene, the first experimentally fabricated 2D material Novoselov et al. (2004), appears to be a perfect material for spintronics Žutić et al. (2004); Fabian et al. (2007) due to extraordinary long mean free path Novoselov et al. (2004); Tombros et al. (2007) and weak spin-orbit coupling (SOC) of carbon atoms. The ideal graphene lattice is flat [$\delta_{z}=0$ in Fig. 1 a)] and belongs to the $D_{6h}$ symmetry point group. The presence of the horizontal mirror plane of the lattice brings serious limitations to spin dynamics, as it forces the spins to be aligned perpendicularly to the graphene’s plane. This is mapped into diagonal in spin basis effective intrinsic SOC Hamiltonian in single band models Min et al. (2006). The corresponding eigenstates are therefore pure spin up and down spinors, and spin scattering is prohibited. This is, however, not the full picture. Including all $p$ orbitals in the Hamiltonian leads to coupling of $\pi$ and $\sigma$ states of the opposite spins even in the presence of the mirror symmetry of the latticeFratini et al. (2013). But even though the states are now spin mixed, there is no effective spin scattering mechanism if graphene lattice remains flat.

The in-plane components of spin can be also present when the mirror symmetry constraint is released, as it takes place in buckled honeycomb materials such as silicene or germenene ($D_{3d}$ point group symmetry). In single band model Hamiltonians, this effect is described by the, so called, intrinsic Rashba SOC Liu et al. (2011); Geissler et al. (2013), or PIA SOC in the context of functionalized graphene Gmitra et al. (2013); Kochan et al. (2017). In contrast to well known Rashba SOC due to structure inversion asymmetry, the intrinsic Rashba SOC does not remove the spin degeneracy of states, as a consequence of preserved space inversion symmetry. Nevertheless, it enables the Elliott–Yafet Elliott (1954); Yafet (1963) (E-Y) spin relaxation mechanism, which allows spin flips only accompanied with momentum scattering by non-magnetic impurities or phonons. The latter are naturally present in rippled structures due to flexular distortions of the lattice Mariani and von Oppen (2008, 2010); Morozov et al. (2008); Castro et al. (2010).

In the E-Y mechanism, the probability of spin flip follows the probability of momentum scattering, $\tau_{s}^{-1}\approx b^{2}\tau_{p}^{-1}$. Elliott (1954) The proportionality factor, $b^{2}$, is the Elliott–Yafet spin mixing (or spin admixture) parameter. It has been extensively studied for bulk materials and thin films of heavy elements Monod and Beuneu (1979); Fabian and Das Sarma (1998); Steiauf and Fähnle (2009); Zimmermann et al. (2012); Long et al. (2013a, b); Zimmermann et al. (2016), but the knowledge about $b^{2}$ in atomically thin 2D systems is very limitedKurpas et al. (2016); Avsar et al. (2017).

In this paper we perform a systematic study of the SOC and spin mixing in elemental 2D materials with a honeycomb lattice structure. We focus on materials made of elements belonging to group 14 and 15 of the periodic table. Starting from an effective multiband symmetry-based Hamiltonian, we revisit the intrinsic SOC in graphene and provide analytical solutions of the eigenstates at the K-point. We show that the expectation value of spin in the Dirac cone bands is smaller than one-half and can differ between valence and conduction band. Next, by using numerical first principles density functional theory methods we characterize intrinsic SOC and calculate spin-mixing parameter $b^{2}$ for graphene, silicene, germanene, stanene, blue phosphorene, arsenene, antimonene and bismuthene. We find, that the strength of the effective intrinsic SOC in the band structure $\lambda_{\text{so}}$ follows a quadratic dependence on the atomic number Z, as expected for valence electrons in isolated atomsShanavas et al. (2014). The spin mixing parameter $b^{2}$ also follows a scaling law, $b^{2}\sim Z^{4}$, except at spin hot spots Fabian and Das Sarma (1998). This parameter exhibits a wide range of values and giant anisotropy.

The paper is organized as follows. In Section II we briefly describe computational methods. In Section III we discuss the effective SOC Hamiltonian of graphene at the K-point and show that its eigenstates are in fact mixtures of spin up and down states. The definition of spin-mixing parameter $b^{2}$ is also given here. The two forthcoming sections, Sec. IV and Sec. V respectively, contain numerical results for graphene, silicene, germanene, stanene, blue phosphorene, arsenene, antimonene, bismuthene with discussion and conclusions.

II Methods

The structure relaxation was prformed in Quantum ESPRESSO package Giannozzi et al. (2009, 2017). For consistency, the PBEsolPerdew et al. (2008) exchange–interaction potential was used for all studied materials. The kinetic energy cut-offs for the wave function and charge density were individually adjusted for each element and are collected in Table 1. A vacuum of 15 Å was introduced to avoid spurious interactions between copies of 2D films. Scalar-relativistic pseudopotentials were used in case of graphene and silicene, whereas for heavier elements the full relativistic pseudopotentials were applied. The force and energy convergence thresholds for ionic minimization were set to $10^{-4}$ Ry/bohr and $10^{-5}$ Ry/bohr respectively. For the Brillouin zone integration a 12$\times$12 $k$-points mesh were generated using the Monkhorst–Pack scheme. The optimized unit cells have been found by minimization of the total energy with respect to the lattice constant $a$. For each value of $a$ internal forces acting on atoms were relaxed using quasi–Newton scheme as implemented in Quantum ESPRESSO. The resulting structure parameters are collected in Table 1.

The calculations of spin properties were performed using the all electron software package Wien2K Blaha et al. (2001). Spin–orbit coupling was included fully relativistically for core electrons, while valence electrons were treated within the second variational step method Singh and Nordstrom (2006). Self-consistency was achieved for a 30$\times$30 $k$-points grid with 91 $k$-points in the irreducible wedge of the Brillouin zone.

III Spin–orbit coupling and spin mixing

III.1 Effective Hamiltonian of intrinsic SOC in graphene

We start the analysis from the SO interaction in graphene. Being the lightest and of the highest symmetry among all the materials considered in this paper, graphene serves as a benchmark for further discussion. In order to understand SOC effects at the Dirac point (K-point) in graphene, let us build a minimal symmetry-based HamiltonianWinkler and Zülicke (2010); Winkler (2003); Lew Yan Voon and Willatzen (2009) by analyzing the direct coupling via the SOC operator between the Dirac cone and the nearby energy bands.

The SOC term is given by

[TABLE]

with the orbital components transforming as pseudovectors, i. e., $H_{SOx}\sim R_{x}$, $H_{SOy}\sim R_{y}$, $H_{SOz}\sim R_{z}$. Considering the energy bands of graphene shown in the Fig. 1(a), we have the Dirac cone bands that belong to the irreducible representation (irrep) $\Gamma_{5}$, with states $\left|\Gamma_{5}^{1}\right\rangle\sim R_{x}$ and $\left|\Gamma_{5}^{2}\right\rangle\sim R_{y}$, the valence band $\Gamma_{6v}$ with states $\left|\Gamma_{6v}^{1}\right\rangle\sim x$ and $\left|\Gamma_{6v}^{2}\right\rangle\sim y$ and the valence band $\Gamma_{1v}$ with state $\left|\Gamma_{1v}\right\rangle\sim 1$. From the symmetry of the states and the operators within the $D_{3h}$ symmetry group of the K-point in graphene, we can find the nonzero matrix elements due to SOC, given by:

[TABLE]

with $\Delta_{5},\Delta_{6}\in\mathbb{R}$ and $\Delta_{51},\Delta_{56}\in\mathbb{C}$.

Writing the SOC Hamiltonian in the basis set

[TABLE]

with \left|\Gamma_{5(6v)}^{\pm}\right\rangle=\left(\left|\Gamma_{5(6v)}^{1}\right\rangle\pm i\left|\Gamma_{5(6v)}^{2}\right\rangle\right)\big{/}\sqrt{2}, we obtain the following block diagonal matrix

[TABLE]

with $E_{1v}<0$, $E_{6v}<0$, $\left|E_{1v}\right|\gg\left|\Delta_{5}\right|$, $\left|E_{6v}\right|\gg\left|\Delta_{5}\right|$ and $\left|E_{6v}\right|\gg\left|\Delta_{6v}\right|$.

Diagonalizing the Hamiltonian, the conduction (subscript $c$) and valence (subscript $v$) band Dirac cones have eigenvalues

[TABLE]

with eigenvectors

[TABLE]

The admixing coefficients are given by

[TABLE]

and

[TABLE]

From the group theory analysis we performed, it is possible to identify two different SOC contributions [depicted by arrows in Fig. 1(a)], the intraband SOC (an interaction of the Dirac cone with itself, couples states with the same spin) and the interband SOC (the direct coupling of the Dirac cone to the valence bands $\Gamma_{1v}$ and $\Gamma_{6v}$, couples states with opposite spins). The effect of these two SOC contributions to the Dirac cone is sketched in Fig. 1(b). If only the intraband SOC is taken into account we notice the opening of the gap and the spin projection of the conduction (identified by the label $c$) and valence (identified by the label v) bands of the Dirac cone is $|\left\langle\hat{s}_{z}\right\rangle_{c}|=|\left\langle\hat{s}_{z}\right\rangle_{v}|=\frac{1}{2}$. Notice that the different Dirac cone branches remain two-fold degenerate in spin and therefore it is enough to discuss the modulus of the spin projection. Adding the interband SOC contribution, the energy gap remains open but now the picture for the spin projection changes. Due to the mixing of the energy bands via the SOC, the eigenstates of the SOC Hamiltonian at the K-point become mixtures of spin up and down states, given in Eq. (III.1). It immediately follows that the spin projection is reduced, $|\left\langle\hat{s}_{z}\right\rangle|<\frac{1}{2}$, and the two branches acquire a slightly different value of spin projection $|\left\langle\hat{s}_{z}\right\rangle_{c}|\neq|\left\langle\hat{s}_{z}\right\rangle_{v}|$. We point out, that although $|\left\langle\hat{s}_{z}\right\rangle|<\frac{1}{2}$ the electron’s spin does not have any component along the $x$ and $y$ directions due to orthogonality of the orbital parts of the states in Eq. (III.1). The mirror symmetry of graphene is thus satisfied by the SOC Hamiltonian. The amplitudes $\beta$ and $\eta$ appear due to weak SOC, thus $|\alpha|\gg|\beta|$ , $|\lambda|\gg|\eta|$, and one can identify $|\beta|^{2}$ and $|\eta|^{2}$ in Eq. (III.1) as Elliott-Yafet spin-mixing parameters discussed below. It is important to note that to obtain the correct value of the spin-orbit gap in graphene, coupling to $d$-orbital bands is needed Gmitra et al. (2009). However, since the relevant $d$ states have spins perpendicular to the plane, this coupling does not contribute, to first order, to the spin mixing. From the symmetry point of view, these $d$ orbitals are already embedded in the $\Gamma_{5}$ states because the symmetry of the energy bands are determined from ab initio, and therefore the mixing of different orbitals are already included in the wave functions.

III.2 Spin-mixing parameter

Let us consider two Bloch spinors $\Psi^{\sigma}_{n,{\mathbf{k}}}({\mathbf{r}})=u_{n,{\mathbf{k}}}({\mathbf{r}})|\sigma\rangle\exp({i{\mathbf{k}}\cdot{\mathbf{r}}})$, were $n$ is the band index, $u_{n,{\mathbf{k}}}({\mathbf{r}})$ is the lattice periodic function, and $\sigma=\uparrow,\downarrow$ is the electron spin. Due to time reversal and space inversion symmetry these states are degenerate at any $\mathbf{k}$-point in the Brillouin zone (BZ), i.e., $E_{n}(\mathbf{k},\uparrow)=E_{n}(\mathbf{k},\downarrow)$. Upon the inclusion of SOC, each of $\psi_{n,\mathbf{k}}^{\sigma}$ acquires an admixture of the opposite spin component forming a new pair of degenerate Bloch states

[TABLE]

where $a_{n,{\mathbf{k}}}({\mathbf{r}})$ and $b_{n,{\mathbf{k}}}({\mathbf{r}})$ are again lattice periodic functions Elliott (1954). Usually $a_{n,{\mathbf{k}}}({\mathbf{r}})$ and $b_{n,{\mathbf{k}}}({\mathbf{r}})$ are chosen in such a way, that $b_{n,{\mathbf{k}}}({\mathbf{r}})$ stands for the coefficient of the small spin component being admixed to the large spin component which has amplitude $a_{n,{\mathbf{k}}}({\mathbf{r}})$, i.e., $|b_{n,{\mathbf{k}}}({\mathbf{r}})|^{2}\ll|a_{n,{\mathbf{k}}}({\mathbf{r}})|^{2}$. Then $\Psi^{\Uparrow}_{n,{\mathbf{k}}}({\mathbf{r}})$ is the wave function of Bloch electrons with the majority spin up and $\Psi^{\Downarrow}_{n,{\mathbf{k}}}({\mathbf{r}})$ the wave function of electrons with the majority spin down. Elliott pointed out Elliott (1954) that the probability of a spin-flip upon momentum scattering is proportional to the spin mixing parameter $b^{2}_{n,{\mathbf{k}}}=\int|b_{n,{\mathbf{k}}}({\mathbf{r}})|^{2}d^{3}{\mathbf{r}}$. The analogy of $|\beta|^{2}$ and $|\eta|^{2}$ to $b^{2}_{n,{\mathbf{k}}}$ is now transparent.

From an experimental point of view, the quantity of interest is the ensemble average of $b^{2}_{n,{\mathbf{k}}}$ for a given Fermi level rather than its value at a single ${\mathbf{k}}$-point. Therefore, it is useful to redefine the Elliott–Yafet spin mixing parameter as the Fermi contour average of $b^{2}_{n,{\mathbf{k}}}$

[TABLE]

where $\hat{\mathbf{s}}$ is the unit vector defining the spin quantization axis (SQA), $S_{BZ}$ is the area of the Fermi surface, $\rho(E_{F})$ is the density of states per spin at the Fermi level, $\mathbf{v}_{F}(\mathbf{k})$ is the Fermi velocity and the integration takes over an iso-energy contour. In electrical spin injection experiments SQA corresponds to the polarization of initial magnetization of populated electrons. Such a definition allows us to explore the anisotropy of $b^{2}$ in the band structure, similarly to what was done for selected 3D materials Zimmermann et al. (2012, 2016).

IV Results and discussion

The initial structure parameters of silicene, germanene, and stanene have been taken from Ref. [Matusalem et al., 2015]. For arsenene we used parameters from Ref. [Kamal and Ezawa, 2015], and for graphene we used the initial lattice constant 2.46 Å. Optimized lattice parameters and buckling heights are very close to the original values and are listed in Table 1.

Spin-orbit splitting.

At first we focus focus on materials from group 14. The band structure of graphene is shown in Fig. 1 and was discussed above. In the top row of Fig. 2 we show the calculated relativistic band structures of silicene, germanene, and stanene. All these materials are semimetals. The semimetalic character is manifested by the presence of a Dirac cone centered in the Brillouin zone at the K-point. Without SOC the valence and conduction branches of the cone touch at the Fermi energy, forming a zero-width band gap [dashed red line in the insets of Fig. 2 (a)].

The inclusion of SOC pulls them apart and introduces a spin-orbital gap $\Delta_{\text{so}}^{\text{K}}$, while the spin degeneracy is preserved by virtue of time reversal and space inversion symmetry. For graphene, silicene, and germanene, the spin-orbital gap $\Delta_{\text{so}}^{\text{K}}$ is synonymous with the fundamental band gap $\Delta$, defined as the energy distance between the valence and conduction band edges [Fig. 2 (d)]. At the $\Gamma$-point, SOC splits off the two top–most valence bands, by the energy $\Delta_{\text{so}}^{\Gamma}$ [Fig. 2 (a),(c)], and other bands lying far away from the Fermi level and being irrelevant to the discussion of low energy physics we focus on here. For graphene, silicene, and germanene the split-off bands at the $\Gamma$-point lie below the valence band maximum (VBM) at the K-point. For stanene [Fig. 2 (c)], due to strong SOC, $\Delta_{\text{so}}^{\Gamma}=200$ meV, the energy of the top-most valence band at the $\Gamma$-point overtakes the energy at the K-point by 47 meV and the band gap becomes indirect.

Contrary to group 14 semimetals, materials made of group 15 elements are semiconductors with sizable indirect band gaps [Fig. 2 (d),(e),(g),(h)]. The calculated values of the band gaps $\Delta$ are: 1.86 eV, 1.49 eV, 0.98 eV and 0.49 eV for blue phosphorene, arsenene, antimonene and bismuthene, respectively (black phosphorene was investigated in Ref. Kurpas et al., 2016). These values are consistent with other calculations Kamal and Ezawa (2015); Zhu et al. (2015b); Kecik et al. (2016), and with experimental reports Shah et al. . Except for bismuthene, the VBM is located at the $\Gamma$-point, while the conduction band minimum (CBM) lies close to the middle of $\Gamma$–M path in the Brillouin zone. The spin-orbital splittings of two top-most valence bands at the $\Gamma$-point $\Delta_{\text{so}}^{\Gamma}$ are 48 meV for blue phosphorene, 195 meV for arsenene and 340 meV for antimonene. Bismuthene [Fig. 2(h)] displays qualitatively different picture due to the inverted band gap. The gap inversion can be easily identified by comparing the band ordering for bismuthene and the remaining materials of group 15. In Fig. 2(i), we show irreducible representations of the four relevant bands of bismuthene. Without SOC the top-most valence band $\Gamma^{+}_{3}$ is two-fold degenerate.

Upon inclusion of SOC it splits off into two bands $\Gamma^{+}_{4}$ and $\Gamma^{+}_{5}+\Gamma^{+}_{6}$ separated by $\Delta_{\text{so}}^{\Gamma}\approx 700$  meV [see Fig. 2 (i)]. The latter band interchanges with the first conduction band $\Gamma_{4}^{-}$ and the gap becomes inverted, with respect to band ordering of lighter materials of group 15 [see Fig. 2 (f)]. The edge of the valence band of bismuthene lies slightly away from the $\Gamma$-point [see inset in Fig. 2(h)], with energy only 17 meV higher than the energy of the band at the $\Gamma$-point, and the character of the band gap is almost direct. For all group 15 materials the characteristic Dirac cone lies approximately 2 eV below the valence band maximum and gradually loses its linear dispersion character with an increasing atomic number Z.

Since all the studied materials have the same crystalline structure one can expect that, within the same group of periodic table, the spin-orbital gap $\Delta_{\text{so}}$ will mainly depend on the electronic configuration of the element. In the first order perturbation theory $\Delta_{\text{so}}\sim\lambda_{\text{so}}$, where $\lambda_{\text{so}}$ is the strength of SOC in the band structure. In isolated atoms, if only valence electrons are taken into account, $\lambda_{\text{so}}\sim Z^{2}$. Shanavas et al. (2014) In crystalline solids, bands close to the Fermi level are made of states of valence electrons. Therefore, one can roughly expect that $\Delta_{\text{so}}$ will also follow a quadratic dependence on Z. On the other hand, the effective SO interaction in a band is, generally, momentum dependent, and includes contributions from other bands coupled by the SO interaction. A systematic study of SO interaction would be necessary to visualize the global behavior, but this goes beyond the scope of this paper. Instead, we will focus only on the high symmetry points K and $\Gamma$. In Fig. 3 we plotted $\Delta_{\text{so}}^{\text{K}}$ and $\Delta_{\text{so}}^{\Gamma}$ as a function of atomic number Z. The values are collected in Table 2. It is seen that indeed $\Delta_{\text{so}}$ follows the $Z^{2}$ dependence very well. A deviation from the common quadratic dependence is seen for graphene and silicene at the K-point. This may be caused by the fact that the core potential is not effectively screened due to low number of core electrons. The deviation for graphene is explained by the absence of buckling and therefore of scalar coupling between $p_{z}$ and in-plane orbitals; that is, the prefactor of the scaling is drastically reduced.

Spin-mixing.

Let us now discuss the spin–mixing parameter $b^{2}$. In Fig. 4 we show the calculated Fermi contour averaged spin–mixing parameter for group 14 materials. For out-of-plane spin polarization (SQA=Z), $b^{2}$ is almost independent of the position of the Fermi energy $\text{E}_{\text{F}}$, both for the valence and conduction band. We relate it to the fact, that around the K-point the two bands forming the Dirac cone are well separated from the others and the effective SOC in the valence and conduction bands near the K-point is almost momentum independent (within the range of doping considered here). The intraband SOC, involving the valence and conduction branch of the cone, does not contribute to $b^{2}$, as was shown by our effective model. The intrinsic Rashba (PIA) SOC vanishes at the K-point and grows linearly with momentum Liu et al. (2011); Geissler et al. (2013); Gmitra et al. (2013); Kochan et al. (2017). Therefore its contribution around the K-point is small. For stanene [Fig. 4 (a)] the valence band edge is at the $\Gamma$-point. Initially very small value of $b^{2}$ rapidly increases with doping, due to interaction with the lower valence band. At $\text{E}_{\text{F}}\approx-47$ meV the band around the K-point starts contributing to the Fermi contour and a discontinuous increase of $b^{2}$ is observed. We stress, that the finite value of $b^{2}$ for graphene, of the order of $10^{-7}$, does not imply spin-flip scattering by scalar impurities. Such scattering is prohibited by the mirror symmetry of the lattice, i.e. $\langle\Gamma_{5}^{\pm}|V_{\rm imp}|\Gamma_{1v}\rangle=\langle\Gamma_{5}^{\pm}|V_{\rm imp}|\Gamma_{6v}^{\mp}\rangle=0$, if the impurity potential $V_{\rm imp}$ is even upon mirror reflection.

For spins polarized in-plane (SQA=X/Y, bottom row in Fig. 4) $b^{2}$ is almost one-half for $\text{E}_{\text{F}}=0$ due to the spin hot spot Fabian and Das Sarma (1998) at the K-point (only $z$ component of spin is allowed). With increasing doping it starts to decrease towards the values similar to SQA=Z. Again, $b^{2}$ in the valence band of stanene is an exception. For the whole doping range it does not go below $b^{2}=0.2$, and spins remain almost fully mixed. Even for high doping, $\text{E}_{\text{F}}=-100$ meV, $b^{2}$ for in-plane spin orientation is two orders of magnitude greater than for out-of-plane spins.

For group 15 materials [Fig. 5] $b^{2}$ displays more diversity due to more complicated band structures around the band gap. Nevertheless, similar trends as for group 14 materials can be identified: (i) for SQA=Z, if the VBM is centered at the $\Gamma$-point, $b^{2}$ grows exponentially when moving away from the Brillouin zone center. This happens in the valence band of blue phosphorene, arsenene, and antimonene as shown in Fig. 5 (a), or in the conduction band of bismuthene [Fig. 5 (b)]. (ii) When the band edge is away from the high symmetry points and the band is relatively well separated from the others (no spin hot spots due to accidental anticrossings occur), $b^{2}$ exhibits a small variation with doping. Such behavior is observed in the conduction band of phosphorene, arsenene, and antimonene [Fig. 5 (b),(d)]. The reduced symmetry in k-space also results in weak anisotropy of $b^{2}$. The values of $b^{2}$ in the conduction and valence bands of phosphorene, arsenene and antimonene, are almost identical. (iii) A discontinuous change of $b^{2}$ takes place when a next band crosses the Fermi level. The contribution of this band to the averaged $b^{2}$ is far from being trivial. It depends on the form and strength of SOC in the band at a given $\mathbf{k}$-point, and on the number of states contributing to the Fermi contour. For example, for stanene $b^{2}$ decreases when the energy band around the K-point starts contributing to the total average [Fig. 4 (a)], and for phosphorene decreases when another valence band crosses the Fermi level [Fig. 5 (a),(c)].

We have also checked how $b^{2}$ scales with the atomic number Z. Within first order non-degenerate perturbation theory, the admixture amplitude in Eq. (10), $b_{n,{\mathbf{k}}}$ is proportional to $\lambda_{\text{so}}$. Taking $\lambda_{\text{so}}\sim\text{Z}^{2}$ one can expect that $b^{2}$ should follow $\text{Z}^{4}$ dependence. In Fig. 6 we plot average $b^{2}$ in the conduction band versus the atomic number Z (corresponding to a given material) for iso-energy contour at $\text{E}_{\text{F}}=60$ meV and SQA=Z [see grey vertical lines in Figs 4 (b) and 5 (b)]. The Fermi energy was chosen such that the spin-mixing parameter is not strongly influenced by the vicinity of a spin hot spot and reflects the pure SOC in the band. As can be seen, our results agree well with the estimate given by perturbation theory, though small deviation from $Z^{4}$ are observed. This can be attributed to a rather complex nature of spin-orbit coupling in many electron crystalline solids.

Finally, we have calculated spin mixing anisotropy, which is a measure of spin relaxation anisotropy. The ratio $b^{2}_{\text{SQA=X}}/b^{2}_{\text{SQA=Z}}$ (in-plane to out-of-plane spin polarization) for materials of group 14 is shown in Fig. 7. All materials display giant (when compared to corresponding anisotropies of 3D materials) and doping dependent anisotropy. For most materials, the anisotropy is driven by the spin hot spot for in-plane polarized spins at the K or $\Gamma$ points. The highest anisotropy at large $\text{E}_{\text{F}}$ is observed for stanene in the valence band, between $10^{2}$ and $10^{3}$, and results in strong spin mixing for in-plane spin polarization. For graphene, we find that $b^{2}_{\text{SQA=X}}/b^{2}_{\text{SQA=Z}}\approx 0.5$ for $E_{F}>30$ meV. Similar trends are observed for materials of group 15. In the valence band [Fig. 8 a)] a spin hot spot at the $\Gamma$ point ($E_{F}=0$ meV) results in a huge anisotropy, which decreases when moving away from the high-symmetry point. This happens for blue phosphorene, arsenene, and antimonene. Bismuthene displays almost no anisotropy of $b^{2}$ in the valence band due to strong spin mixing for all spin polarizations. The picture is oposite in the conduction band [Fig. 8 b)]. Anisotropic behavior of $b^{2}$ is found for blue phosphorene, arsenene, and antimonene, while $b^{2}$ for bismuthene shows doping dependent anistropy. Essentially, $b^{2}$ exhibits a strong anisotropy if the BZ wedge defined by the doping range contains spin hot spots or spin hot regions, while the anisotropy is not well pronounced otherwise.

The above calculated spin admixture can be used to obtain realistic estimates of spin relaxation times. Indeed, all the studied elemental 2D materials have space inversion symmetry and are thus expected to exhibit spin relaxation according to the Elliott-Yafet mechanism Elliott (1954); Yafet (1963). The only other input needed is the momentum relaxation time $\tau_{p}$ which can be obtained from electrical transport, for example. The link to spin relaxation is provided by the Elliott relation $\tau_{s}^{-1}\approx b^{2}\tau_{p}^{-1}$, which should be valid in regions where $b^{2}\lesssim 0.2$, where perturbation theory holds. Two connected basic assumptions are needed: the spin-orbit coupled bands should be spectrally separated by more than is the spin-orbit coupling matrix element between them, and the spin expectation value is close to 1/2.Elliott (1954)

At spin hot spots, which occur at $K$ and $\Gamma$ points for our studied materials, these two assumptions are in general violated and the mechanism needs to be modified, see for example Ref. [Boross et al., 2013; Szolnoki et al., 2017]. At these points, for the spin direction at which $b^{2}$ is of order 1, spin relaxation and and momentum relaxation times become comparable Fabian and Das Sarma (1999).

V Conclusions

We have performed a systematic study of spin-orbit coupling in elemental two-dimensional materials of group 14 and 15 of the periodic table. Starting from symmetry arguments we have formulated an effective multiband symmetry-based SOC Hamiltonian for graphene at the K-point. We have shown that even if the mirror symmetry of the lattice protects the spin in graphene from acquiring the $x$ and $y$ components, spin mixing due to the intrinsic SOC is still possible but does not lead to spin relaxation. Using first principles numerical methods we analyzed intrinsic SOC and calculated the Elliott-Yafet spin-mixing parameter $b^{2}$ for graphene and other honeycomb lattice materials. We have shown that spin-orbit coupling in the band structure scales as a square function of the atomic number Z. Away from spin hot spots the spin-mixing parameter also follows the exponential scaling power law, $b^{2}\sim\text{Z}^{4}$. We identified three main factors having the strongest influence on the overall average value of $b^{2}$, namely, the strength of the intrinsic SOC, the shape of the Fermi contour, and the presence of spin hot spots inside or close to the contour. For almost all materials $b^{2}$ shows substantial and doping dependent anisotropy. Our results for $b^{2}$ can be translated into spin relaxation times, once the momentum relaxation time is known. Therefore they provide valuable information on the potential application of those materials in spintronics.

Acknowledgments

This work was supported by the National Science Centre under the contract DEC-2018/29/B/ST3/01892, Alexander von Humboldt Foundation, Capes (grant No. 99999.000420/2016-06), SFB 1277 (A09 and B05) in part by PAAD Infrastructure co-financed by Operational Programme Innovative Economy, Objective 2.3, and by the Scientific Grant Agency of Ministry of Education of Slovak Republic (Grant No. 1/0531/19) and VVGS-2018-887. The authors gratefully acknowledge the Gauss Center for Supercomputing e.V. for providing computational resources on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Center.