Spin-orbit coupling and spin relaxation in phosphorene: Intrinsic versus extrinsic effects

TL;DR

This study uses first-principles calculations to analyze intrinsic and extrinsic spin-orbit effects in phosphorene, revealing anisotropic spin relaxation mechanisms and their dependence on electric fields and substrates.

Contribution

It provides a detailed comparison of intrinsic and extrinsic spin relaxation mechanisms in phosphorene, highlighting the role of electric fields and substrate effects.

Findings

Intrinsic spin mixing probability is about 10^{-4} with large anisotropy.

Intrinsic spin relaxation times range from hundreds of picoseconds to nanoseconds.

Extrinsic spin relaxation via D'yakonov--Perel' mechanism can dominate under strong electric fields.

Abstract

First-principles calculations of the essential spin-orbit and spin relaxation properties of phosphorene are performed. Intrinsic spin-orbit coupling induces spin mixing with the probability of $b^2 \approx 10^{-4}$, exhibiting a large anisotropy, following the anisotropic crystalline structure of phosphorene. For realistic values of the momentum relaxation times, the intrinsic (Elliott--Yafet) spin relaxation times are hundreds of picoseconds to nanoseconds. Applying a transverse electric field (simulating gating and substrates) generates extrinsic $C_{2v}$ symmetric spin-orbit fields in phosphorene, which activate the D'yakonov--Perel' mechanism for spin relaxation. It is shown that this extrinsic spin relaxation also has a strong anisotropy, and can dominate over the Elliott-Yafet one for strong enough electric fields. Phosphorene on substrates can thus exhibit an interesting interplay of both spin relaxation mechanisms, whose individual roles could be deciphered using our results.

Tables

Table 1. Table 1: Structural parameters for relaxed phosphorene and experimental data for bulk black phosphorus.

	2D relaxed	bulk exp. (Brown and Rundqvist, 1965)
a	3.2986 Å	3.3136
b	4.6201 Å	4.376
$d_1$	2.2223 Å	2.224
$d_2$	2.2601 Å	2.244
${\alpha }_1$	95.833	96.34
${\alpha }_2$	104.085	102.09

Table 2. Table 2: Spin-orbit coupling parameters extracted from first-principles calculations for different values of the electric field.

	Valence band		Conduction band
E [V/nm]	$\gamma$ [meVÅ]	$\alpha$ [meVÅ]	$\gamma$ [meVÅ]	$\alpha$ [meVÅ]
1.0	0.46	-0.23	0.29	-0.06
2.0	0.91	-0.48	0.56	-0.11
3.0	1.36	-0.72	0.84	-0.17
5.0	2.21	-1.13	1.43	-0.29

Equations

$$b^2$$

$$b_{n,\mathbf{k}}^2$$

$${\tau }_{\mathrm{s},\mathrm{EY}}^{-1}\approx 4b^2/{\tau }_p^{-1},$$

$$H_{\mathrm{ex}}=\frac{\hslash }{2}{\mathbf{\Omega }}_{\mathbf{k}}\cdot \mathbf{\sigma },$$

$${\tau }_{\mathrm{s},\mathrm{DP}}^{-1}={\Omega }_{\perp }^2{\tau }_p,$$

$${\Omega }_x^2$$

$${\Omega }_y^2$$

Full text

Phosphorene is a monolayer of black phosphorus Brown and Rundqvist (1965); Cartz et al. (1979); Keyes (1953); Maruyama et al. (1981); Narita et al. (1983), exhibiting a direct band gap of 2 eV Castellanos-Gomez et al. (2014); Liang et al. (2014) and large anisotropic mobility Castellanos-Gomez et al. (2014); Liu et al. (2014); Qiao et al. (2014). Unlike graphene, phosphorene is a semiconductor, and unlike two-dimensional transition-metal dichalcogenides, which are semiconductors too, phosphorene is distinctly anisotropic thanks to its puckered atomic structure. The semiconductor property makes phosphorene suitable for electronic Li et al. (2014) and spintronics applications Žutić et al. (2004); Fabian et al. (2007), in particular for bipolar spin diodes and transistors Žutić et al. (2006), while the anisotropy enables directional control of the essential spin properties, such as spin-orbit coupling and spin relaxation. In contrast to graphene, whose spin properties are by now well established Han et al. (2014), there is no unified picture of the spin-orbit coupling and spin relaxation in phosphorene.

Phosphorene can be extracted from black phosphorus by mechanical (Xia et al., 2014; Liu et al., 2014; Li et al., 2014) or liquid (Brent et al., 2014; Kang et al., 2015) cleavage techniques. Inside phosphorene layers, each phosphorus atom is covalently bonded with three adjacent phosphorus atoms to form a puckered honeycomb structure due to $sp^{3}$ hybridization, see Fig. 1. The puckered structure can be viewed as a two-layer system in which the bonding energy is dominated by the in-plane bonds ($pp\sigma$ and $pp\pi$) that are much stronger than the bonds connecting the two sublayers ($pp\pi$) (Li and Appelbaum, 2014). Similarly to graphene, the edges of phosphorene form zig-zag (along $x$ axis) and armchair (along $y$ axis) chains [see Fig. 1b)].

Black phosphorus is described by the nonsymmorphic $D_{2h}$ point group being isomorphic with the $C_{mca}$ space group. Phosphorene shares the same point group symmetry as its bulk counterpart. Both structures have inversion symmetry leading to spin degenerate eigenstates. Spin-orbit coupling leads to the spin-mixing of the Pauli spinors—the intrinsic effect. When inversion symmetry of phosphorene is broken by an applied transverse electric field $\mathbf{E}$ or a substrate, the point group is reduced to nonsymmorphic $C_{2v}$, with the principal $C_{2}$ axis parallel to the direction of the electric field and two mirror planes $\sigma_{xz}$ and $\sigma_{yz}$ [Fig. 1 (a,b)]. In this case the spin degeneracy is lifted—the extrinsic (Rashba) effect.

From the spintronics perspective two questions are particularly important to address: (i) what is the intrinsic and extrinsic spin-orbit coupling (SOC) in phosphorene and (ii) what are the relevant spin relaxation time scales. The answer to the first question has been partially answered by Popović et al. Popović et al. (2015) who showed that the extrinsic Rashba effect, due to external electric fields, is anisotropic with respect to the two principal directions in the crystal. The second question has been addressed within $k\cdot p$ theory Li and Appelbaum (2014) for the intrinsic effects only.

Here we employ first-principles calculations to address both questions, providing state-of-the-art most realistic results for the extrinsic and intrinsic effects. First, we find that intrinsic SOC lifts degeneracy of the valence and conduction bands at the S point, by splitting the bands of about 17.5 meV and 14 meV, respectively. The extrinsic Rashba SOC is much weaker, of the order of tens of $\mu$eV close to the $\Gamma$ point (for electric fields of 1 V/nm), and is found to be significantly anisotropic for the valence band only. Second, we predict the spin lifetime in phosphorene to be hundreds of picoseconds up to nanoseconds, for the experimentally relevant mobilities. We find that for no and small electric fields up to ($\text{E}\leq 2.5$ V/nm) and carrier densities up to $8\cdot 10^{12}$ cm*-2*, the dominant spin relaxation mechanism is the Elliott–Yafet Elliott (1954); Yafet (1963). For the in-plane spin orientation the relaxation is almost twice the slower than for spins oriented out-of-plane. By increased transverse electric field the D’yakonov–Perel’ Dyakonov and Perel (1971) mechanism starts to be the most effective. For carrier density $n\approx 3\cdot 10^{12}$ cm*-2* it overtakes the Elliott–Yafet’s at $E=4$ V/nm for holes and $E=5$ V/nm for electrons. As a result, the substrates can be essential for spin dynamics in phosphorene due to an interplay between the D’yakonov–Perel’ and the Elliott–Yafet relaxation mechanisms.

Intrinsic phosphorene belongs to the family of centrosymmetric crystals for which the dominant spin relaxation mechanism is the Elliott–Yafet scattering Elliott (1954); Yafet (1963). The strength of the scattering can be quantified by the spin admixture coefficient $b^{2}$ Elliott (1954), which is a signature of intrinsic spin-orbit coupling. The time-reversal and space-inversion symmetries require two Bloch states $\Psi_{n,k}^{\uparrow}(\mathbf{r})$ and $\Psi_{n,k}^{\downarrow}(\mathbf{r})$ of the same band $n$ and momentum $k$ to be degenerate (Kramer’s doublets). Due to spin-orbit coupling these states are mixtures of spin up $|\uparrow\rangle$ and spin down $|\downarrow\rangle$ Pauli spinors: $\Psi_{n,\mathbf{k}}^{\uparrow}(\mathbf{r})=\left[a_{n,\mathbf{k}}(\mathbf{r})|\uparrow\rangle+b_{n,\mathbf{k}}(\mathbf{r})|\downarrow\rangle\right]e^{i\mathbf{k}\mathbf{r}}$, $\Psi_{n,\mathbf{k}}^{\downarrow}(\mathbf{r})=[a_{n,\mathbf{-k}}^{\ast}(\mathbf{r})|\downarrow\rangle-b_{n,\mathbf{-k}}^{\ast}(\mathbf{r})|\uparrow\rangle]e^{i\mathbf{k}\mathbf{r}}$. For a generic Bloch state the modulation functions $a_{n,\mathbf{k}}(\mathbf{r})$ and $b_{n,\mathbf{k}}(\mathbf{r})$ are selected to diagonalize the spin magnetic moment along the chosen direction, corresponding to the injected spin in experiment; for weak spin-orbit coupling, $b_{n,\mathbf{k}}(\mathbf{r})$ stands for the small spin component being admixed to the large spin component $a_{n,\mathbf{k}}(\mathbf{r})$, i. e., $|a_{n,\mathbf{k}}(\mathbf{r})|^{2}\gg|b_{n,\mathbf{k}}(\mathbf{r})|^{2}$. The Eliott–Yafet scattering parameter $b^{2}$ is defined as the Fermi surface average of the unit cell integrated admixture coefficient $b^{2}_{n,\mathbf{k}}$,

[TABLE]

where $0\leq b^{2}\leq 0.5$, $\rho(\varepsilon_{\rm F})$ is the density of states per spin at the Fermi level, $v_{\rm F}$ is the Fermi velocity, and $S_{\rm BZ}$ is the area of the Brillouin zone. If the scattering potential is spin-independent (scalar impurities and phonons), the intrinsic SOC leads to spin-flip scattering. The Elliott–Yafet mechanism gives for the spin relaxation rate, Elliott (1954); Fabian and Das Sarma (1998)

[TABLE]

where $\tau_{p}^{-1}$ is the momentum relaxation rate.

Extrinsic effects appear once the space inversion symmetry gets broken, e.g., by a substrate or external fields. The spin degeneracy gets lifted and another spin relaxation mechanism appears: D’yakonov–Perel’ Dyakonov and Perel (1971). This mechanism can be viewed as a motional narrowing of the spin precession in a fluctuating (due to momentum scattering) emerging spin-orbit field $\mathbf{\Omega}_{\mathbf{k}}$, which is related to the spin splitting as

[TABLE]

where $\boldsymbol{\sigma}$ is the vector of Pauli matrices. In the relevant limit of small correlation times ($\Omega\tau_{p}\ll 1$), i. e., when the precession angle between the scattering events is small, the spin relaxation rate becomes, Dyakonov and Perel (1971)

[TABLE]

where $\Omega_{\perp}^{2}=\langle\Omega^{2}_{\mathbf{k},\perp}\rangle$ is the Fermi contour average of the squared spin-orbit field projected to the plane perpendicular to the spin orientation.

In our first-principles calculations we used the initial crystal structure parameters from Ref. Brown and Rundqvist (1965) for bulk black phosphorus. A sheet of phosphorene was placed in vacuum of 20 Å and fully relaxed using quasi–Newton variable–cell scheme as implemented in the Quantum Espresso (Giannozzi et al., 2009) package. Positions of atoms have been relaxed in all directions with the force convergence threshold $10^{-4}$ Ry/a.u. and total energy convergence condition $10^{-5}$ Ry/a.u.. The norm-conserving pseudopotential, with kinetic energy cutoffs of 70 Ry and 280 Ry for the wavefunction and charge density respectively, has been used along with the PBE exchange-correlation functional Perdew et al. (1996). Obtained structural parameters are summarized in Ref. see Suplementary information .

Further electronic structure calculations have been performed using the full-potential linearized augmented plane-wave method as implemented in all-electron code package Wien2k Blaha et al. (2001). Self-consistency has been achieved for $16\times 12\times 1$ Monkhorst-Pack $k$-point grid with 151 $k$–points in the irreducible wedge of the Brillouin zone. SOC has been included fully relativistically for core electrons while five valence electrons have been treated within second variational step method Singh and Nordstrom (2006). For the calculations with the transverse electric field we considered vacuum size of 25 Å. It is known that standard DFT methods underestimate the bandgap of semiconductors. Theoretical bandgaps of phosphorene spread between $0.7$ eV and $2.2$ eV depending on the method of calculationsRodin et al. (2014); Qiao et al. (2014); Liu et al. (2014); Hu et al. (2014); Rudenko and Katsnelson (2014); Tran et al. (2014). On the other hand, recent experimental reports suggest the bandgap of phosphorene of about 2 eV Liang et al. (2014); Wang et al. (2015). It has been reported for standard semiconductors Chantis et al. (2006) that the underestimation of the bandgap impairs significantly SOC effects. Therefore, to consider a realistic bandgap we perform the calculations with undressed LDA functional along with the modified Becke–Johnson potential Tran and Blaha (2009) parametrized to give the bandgap of $2.17$ eV.

The calculated band structure of phosphorene is shown in Fig. 2. We get a direct gap at the zone center. Some DFT calculations Popović et al. (2015); Rodin et al. (2014); Ziletti et al. (2015) report a nearly indirect bandgap, with a somewhat displaced valence band maximum. The valence band in the vicinity of the $\Gamma$ point along $k_{x}$ is nearly dispersionless while it is very dispersive in the $k_{y}$ direction. Similar, but substantially smaller dispersion anisotropy is seen in the conduction band. Close to the $\Gamma$ point both, the valence and the conduction band have mainly $p_{z}$ orbital character, the latter having a small admixture of $p_{y}$ orbitals see Suplementary information . The next conduction band minimum appears at energy $50$ meV above the global conduction band minimum (in direction towards the X point) and consists mainly of $p_{x}$ and $p_{y}$ orbitals with an admixture of $d_{x^{2}-y^{2}}$ electrons. The small distance of this band to the conduction band minimum is reflected as an increase of the slope in the carrier density $n(\varepsilon_{\rm F})$ shown in Fig. 2(e) as a function of the Fermi level $\varepsilon_{\rm F}$. For the valence band the carrier density is a smooth quadratic function of $\varepsilon_{\rm F}$. Empirical fits for the dependence $n(\varepsilon_{F})$, which should be useful for interpreting experiments, are given in the caption to Fig. 2.

Intrinsic effects.

We first discuss the intrinsic SOC and the Elliott–Yafet spin relaxation.

The intrinsic SOC in phosphorene is relatively strong but does not modify substantially the bandstructure close to the bandgap. The orbital degeneracy of the fourfold degenerate bands, see inset to Fig. 2(c), is split into two pairs of spin degenerate bands. The splitting is maximal at the S point, $17.5$ meV and $14$ meV for valence and conduction band respectively, gradually decreasing towards the time reversal points X and Y. The states at the Brillouin zone edges sticking together due to nonsymmorphicity of the $D_{2h}$ group Dresselhaus et al. (2008).

The important effect of the intrinsic SOC is the spin mixing, quantified by $b^{2}_{\mathbf{k}}$. In Fig. 3(a,b) we show the distribution of spin mixing parameter $b^{2}_{\mathbf{k}}$ in the first Brillouin zone of phosphorene for the spin quantization axis oriented perpendicular to the 2D plane. The other spin orientations are discussed in Ref. see Suplementary information . For momenta corresponding to anticrossings and at the BZ edges (except the points X and Y for which $b^{2}_{\mathbf{k}}$ is zero) the values of $b^{2}_{\mathbf{k}}$ are close to $\frac{1}{2}$. The Bloch eigenstates here are fully spin-mixed, forming spin hot spots Fabian and Das Sarma (1998); Fabian and Sarma (1999). At the zone center $b^{2}_{\mathbf{k}}$ is about $10^{-4}$. Perturbation theory gives that $b$ is roughly the ratio of the intrinsic spin-orbit coupling (order 10 meV) and the band gap (order 1 eV), thus $b\sim 0.01$, matches well to the calculated first-principles value of $b^{2}\sim 10^{-4}$. For comparison, the $k\cdot p$ theory gives $b^{2}$ between $10^{-6}$ to $10^{-5}$ Li and Appelbaum (2014).

Knowing $b^{2}$ we now calculate the Elliott–Yafet spin relaxation rates using Eq. (3). For the momentum relaxation we take the typical experimental value of $\tau_{p}=100$ fs. The results can be easily rescaled for the actual experimental mobilities. Calculated $\tau_{\rm s,EY}^{-1}$ as a function of carrier density $n(\varepsilon_{\rm F})$, for valence and conduction bands and different spin quantization axes, are shown in Fig. 3(c,d). The relaxation rates are almost independent of $n$, which follows $b^{2}$ since we use a constant momentum relaxation time. The monotonicity of $\tau_{\rm s,EY}^{-1}$ is then unambiguously determined by $b^{2}$. The spin relaxation rates of holes are greater than those of electrons. Most striking is the strong anisotropy. The largest spin relaxation is for out-of-plane spins, which relax roughly twice as fast as the in-plane spins. We predict the longest spin lifetimes for armchair-oriented spins in the conduction band, and zigzag-oriented spins in the valence band. In the recent $k\cdot p$ theory Li and Appelbaum (2014), the estimated ratio between the spin relaxation rate for out-of-plane to in-plane spins was $\sim 4$, which is an overestimation in view of our first-principles results, but is in a qualitative agreement. Similar anisotropies in Elliott–Yafet spin lifetimes were also predicted for anisotropic bulk materials and thin metallic films Zimmermann et al. (2012); Long et al. (2013).

Extrinsic effects.

In realistic situations phosphorene sits on a substrate or is studied in a gating electric field which breaks space inversion symmetry ($D_{2h}\longrightarrow C_{2v}$). An extrinsic Rashba spin-orbit field emerges, lifting the spin degeneracy, $\varepsilon_{\mathbf{k}\uparrow}\neq\varepsilon_{\mathbf{k}\downarrow}$, according to Eq. (4), except at time-reversal invariant points. Emerged spin-orbit fields give rise to spin relaxation due to the D’yakonov–Perel’ mechanism, which competes with the Elliott–Yafet spin-flip scattering. Here we model the symmetry breaking by applying a transverse electric field, all within the first-principles calculations, of 1 V/nm. In Fig. 4(a,b) we plot the spin-orbit field $\Omega_{\mathbf{k}}$ magnitude over the first Brillouin zone. Similarly to $b^{2}_{\mathbf{k}}$, the values of $\Omega_{\mathbf{k}}$ are peaked at the bands anticrossings and at the BZ edges. At the time-reversal points the $\Omega_{\mathbf{k}}$ is zero.

We note, that for a bare PBE exchange-correlation functional Perdew et al. (1996) (band gap $E_{g}\approx 1$ eV), the Rashba spin-orbit coupling due to external electric fields exhibits a strong anisotropy in the valence as well as in the conduction band Popović et al. (2015). Our calculations show that increasing the bandgap to the experimental value $\approx 2$ eV removes the anisotropy from the conduction band, while it is preserved for the valence band.

To obtain the spin relaxation rates for the D’yakonov–Perel’ mechanism, we resolve the coordinate components of the vector spin-orbit fields $\Omega$ which lie in the phosphorene plane: $\Omega_{x}$ along $x$ (zigzag) and $\Omega_{y}$ along $y$ (armchair) directions. We extract these components by fitting an effective $C_{2v}$ symmetric spin-orbit coupling Hamiltonian to the first-principles data see Suplementary information . In Fig. 4(c,d) we show the calculated spin relaxation rates using Eq. (5), assuming $E=1$ V/nm and $\tau_{p}=100$ fs. The spin lifetime is exceptionally long, of a few ns, and exceeds the lifetime from the Elliott–Yafet mechanism. The relaxation rates for conduction electrons are twice smaller than for the valence electrons. With a growing electric field, the D’yakonov–Perel’ mechanism becomes more significant. For the valence band, it surpasses the Elliot–Yafet’s for $n\gtrsim 6\cdot 10^{12}$ cm*-2*, $n\gtrsim 3\cdot 10^{12}$ cm*-2* and $n\gtrsim 2\cdot 10^{12}$ cm*-2* for electric fields $E=3$ V/nm, $E=4$ V/nm and $E=5$ V/nm respectively. For the conduction band the transitions happen for slightly higher carrier densities. For electric fields $E=5$ V/nm the D’yakonov–Perel’ spin relaxation rates are $\approx 23$ times bigger than for $E=1$ V/nm. Similarly to $b^{2}$, the spin-orbit field $\Omega^{2}$ reveals a strong anisotropy. As a result the in-plane spins relax about $1.5-3$ times slower than the out-of-plane ones.

In summary, we have studied intrinsic and extrinsic spin-orbit coupling and spin relaxation mechanisms in phosphorene. The Elliott–Yafet spin relaxation gives spin lifetimes less than nanoseconds for experimentally relevant samples. The D’yakonov–Perel’ mechanism matters at large electric fields. The lifetimes exhibit a large anisotropy for in-plane and out-of-plane spin orientations.

We acknowledge funding from DFG SPP 1538, SFB 689, NCN DEC-2013/11/B/ST3/00824 and the EU Seventh Framework Programme under Grant Agreement No. 604391 Graphene Flagship.

Supplementary information

Structural parameters

Structure parameters for phosphorene obtained by structural relaxation Giannozzi et al. (2009) and experimental values for black phosphorus Brown and Rundqvist (1965) are given in Tab. 1. The lattice constants, interlayer bond distances and bond angles are sketched in Fig. 5.

Spin–mixing parameters

In Fig. 6 we plot the distribution of the spin-mixing parameter $b^{2}_{\mathbf{k}}$ in the first Brillouin zone for phosphorene and spin quantization axis along the zigzag [Fig. 6(a,b)] and the armchair direction [Fig. 6(c,d)]. The differences between the two spin quantization axes are remarkable. For the zigzag spin quantization axis the spin hot spots are located both at the bands anticrossings and at the BZ edges, whereas for the armchair one the hot regions are only at the anticrossings. Close to the zone center, around the maximum of the valence band and the minimum of the conduction band, the differences in $b^{2}_{\mathbf{k}}$ are not so evident. The corresponding values of $b^{2}$, calculated versus carrier density, are similar for the zigzag and armchair spin quantization axes, as shown in Fig. 7. The parameter $b^{2}$ takes the largest values for the spins oriented out-of-plane, giving the shortest spin lifetimes. The spin quantization axis with the minimal values of $b^{2}$ (and the longest spin lifetime) is along the zigzag and the armchair edges of phosphorene for the valence and the conduction band respectively. The ratio between maximal and minimal $b^{2}$ is roughly 2.25 for holes and 1.8 for electrons.

Effective spin–orbit Hamiltonian and spin–orbit fields

The extrinsic SOC effects close to the $\Gamma$ point can be described by the effective spin-orbit Hamiltonian preserving $C_{2v}$ symmetry, $H_{\rm soc}=\alpha\left(\sigma_{y}k_{x}-\sigma_{x}k_{y}\right)+\gamma\left(\sigma_{y}k_{x}+\sigma_{x}k_{y}\right)$, where $\alpha$ and $\gamma$ are the SOC parameters, akin for Rashba and Dresselhaus SOC in semiconductor heterostructures Bychkov and Rashba (1984), $k_{x}$ and $k_{y}$ are the components of the wave vector, and $\sigma_{x}$ and $\sigma_{y}$ are the Pauli matrices. Diagonalizing $H_{\rm soc}$ we can express the energy spin splitting $\Delta\varepsilon(\mathbf{k})=2\sqrt{(\gamma+\alpha)^{2}k_{x}^{2}+(\gamma-\alpha)^{2}k_{y}^{2}}$ along principal axes in the form $\Delta\varepsilon=2|\alpha+\gamma|k_{x}$ and $\Delta\varepsilon=2|\alpha-\gamma|k_{y}$. The calculated spin-orbit coupling parameters $\alpha$ and $\gamma$ are listed in Tab. 2 for several experimentally relevant values of transverse electric field.

In Fig. 8 we plot the spin splitting $\Delta\varepsilon$ in the vicinity of the $\Gamma$ point for the valence and conduction bands for several values of the transverse electric field. For the valence band the splittings along $\Gamma$–X are twice smaller that for $\Gamma$–Y path. The anisotropy is much less pronounced in the conduction band.

Having the parameters $\alpha$ and $\gamma$ of the effective Hamiltonian $H_{\rm soc}$ one can calculate averages of squared spin–orbit field components: $\Omega^{2}_{x}$, $\Omega^{2}_{y}$, and $\Omega^{2}=\Omega^{2}_{x}+\Omega^{2}_{y}$

[TABLE]

where $k_{x}$ ($k_{y}$) is the momentum of carriers moving along zigzag (armchair) direction of phosphorene; the average is taken over the Fermi contour. The spin-orbit fields are plotted in Fig. 9 as a function of the carrier doping. The calculated values of $\Omega^{2}$ for the valence band are roughly twice as large as for the conduction band. A significant anisotropy with respect to main crystal axes is observed giving the ratio between the maximal ($\Omega^{2}$) and the minimal ($\Omega^{2}_{x}$) values about 2.5 for the valence and 3 for the conduction band.

Elliott–Yafet versus D’yakonov–Perel’

We compare the spin scattering rates from the Elliott–Yafet and the D’yakonov–Perel’ mechanisms for spins pointing along zigzag direction in Fig. 10. For electric fields $E\leq 2$ V/nm the Elliott–Yafet dominates spin scattering of electrons and holes.

With increasing electric field the contribution of the D’yakonov–Perel’ mechanism to spin relaxation increases and depends on carrier type and density. For $E=5$ V/nm it overtakes the Elliott–Yafet’s for carrier densities $n\geq 2\cdot 10^{12}$ cm*-2* for holes and $n\geq 3\cdot 10^{12}$ cm*-2* for electrons. Increasing carrier concentration the relaxation rates grow up to 25 ns*-1*, see Fig. 11(a,b), and the spin lifetime becomes mainly limited by the D’yakonov–Perel’ mechanism for all spin directions and carrier types.

Orbital resolved bandstructure

In figure 12 we show calculated $s$, $p$ and $d$ orbital resolved band structure plots along the high symmetry lines in the first Brillouin zone. The valence and conduction band edges near the zone center are formed by $p_{z}$ orbitals. The first valence band preserves its $p_{z}$ orbital character within bandwidth of about 2 eV and then further mix with other bands of $p_{x}$ and $p_{y}$ orbital character. In detail we see that the first valence band along all high symmetry lines is almost purely of $p_{z}$ character with a small contribution from the $s$ electrons, see Fig. 12(c,d). Exception is along the $\Gamma$-$Y$ line when the $p_{y}$ band rises in energy with increasing momentum $k$. At the $Y$ point the two $p_{y}$ bands merge and stick together at the zone edge towards the $S$ point. Along this way the stuck bands change their character from $p_{y}$ to $p_{z}$. At the energies about $2$ eV below the Fermi level, a contribution from $d_{xy}$ orbital character is found which is an admixture to the $p_{x}$ band, Fig. 12(a) and (f).

For the first conduction band the situation is similar to the first valence band. An important difference is noted along the $\Gamma$-$X$ line where another parabolic-like band of mixed $p_{x}$ and $p_{y}$ character crosses the $p_{z}$ band. The energy offset of the band minimum is about 50 meV above the conduction band edge at the zone center. Presence of this band influences charge carrier concentration as shown in Fig. 2(e) in the main text. Contrary to the valence band a contribution from the $d$ electrons in the conduction band manifold is found at much lower energies. At about 0.5 eV from the conduction band edge there is a significant contribution from the $d_{x^{2}-y^{2}}$ electrons to the second conduction band, see Fig. 12.