$k \cdot p$ theory for phosphorene: Effective g-factors, Landau levels, and excitons

TL;DR

This paper develops comprehensive $k \, p$ models for phosphorene, accurately capturing its electronic, optical, and magnetic properties, including spin-orbit effects, excitons, and Landau levels, validated against ab initio data.

Contribution

The paper introduces detailed $k \, p$ Hamiltonians for phosphorene that incorporate interband spin-orbit coupling, enabling precise analysis of its anisotropic electronic and optical properties.

Findings

Predicted electron and hole g-factor correction of ~0.03.

Calculated exciton binding energy of 0.81 eV.

Validated models against experimental and ab initio data.

Abstract

Phosphorene, a single layer of black phosphorus, is a direct-band gap two-dimensional semiconductor with promising charge and spin transport properties. The electronic band structure of phosphorene is strongly affected by the structural anisotropy of the underlying crystal lattice. We describe the relevant conduction and valence bands close to the $\Gamma$ point by four- and six-band (with spin) $k \cdot p$ models, including the previously overlooked interband spin-orbit coupling which is essential for studying anisotropic crystals. All the $k \cdot p$ parameters are obtained by a robust fit to {\it ab initio} data, by taking into account the nominal band structure and the $k$-dependence of the effective mass close to $\Gamma$-point. The inclusion of interband spin-orbit coupling allows us to determine dipole transitions along both armchair and zigzag directions. The interband coupling is also key to determine the effective g-factors and Zeeman splittings of the Landau levels. We predict the electron and hole g-factor correction of $\approx 0.03$ due to the intrinsic contributions in phosphorene, which lies within the existing range of experimental data. Furthermore, we investigate excitonic effects using the $k \cdot p$ models and find exciton binding energy (0.81 eV) and exciton diameters consistent with experiments and {\it ab initio} based calculations. The proposed $k \cdot p$ Hamiltonians should be useful for investigating magnetic, spin, transport, optical properties and many-body effects in phosphorene.

Tables

Table 1. Table 1: Values of the k ⋅ p ⋅ 𝑘 𝑝 k\cdot p parameters for the different k ⋅ p ⋅ 𝑘 𝑝 k\cdot p models. The parameters P 𝑃 P and α 𝛼 \alpha have units of eV ⋅ Å ⋅ absent Å \cdot\textrm{\AA} and the parameters A 𝐴 A and B 𝐵 B are dimensionless.

	Nph4	ph4	ph6
$A_{c1}$	0.8696	0.8208	0.8405
$B_{c1}$	4.1667	0.5453	0.4800
$A_{v1}$	-0.1372	-0.1614	-0.0958
$B_{v1}$	-4.1667	-0.4588	-0.2926
$P_{c1v1}$		5.3696	5.4473
${\alpha }_{c1v1}$		-0.0195	-0.0218
$P_{c2c1}$			0.2051
${\alpha }_{c2c1}$			0.0070
$A_{c2}$			7.0622
$B_{c2}$			0.4845
$A_{c2v1}$			1.9669
$B_{c2v1}$			0.8598

Table 2. Table 2: Calculated values for the effective g-factors using the ph4 and ph6 k ⋅ p ⋅ 𝑘 𝑝 k\cdot p models. The values in parentheses indicate the calculated g-factors with reversed signs of the α 𝛼 \alpha parameters given in Table  1 .

	ph4	ph6
$g_{v1}$	2.0276 (1.9770)	2.0309 (1.9737)
$g_{c1}$	2.0276 (1.9770)	2.0295 (1.9752)
$g_{c2}$		2.0009 (2.0038)

Table 3. Table 3: Character table for the nonsymmorphic group D 2 ​ h subscript 𝐷 2 ℎ D_{2h} of phosphorene.

	$E$	$\left\{C_z\right|\overrightarrow{\tau }\}$	$C_x$	$\left\{C_y\right|\overrightarrow{\tau }\}$	$I$	$\left\{I{C}_z\right|\overrightarrow{\tau }\}$	$I{C}_x$	$\left\{I{C}_y\right|\overrightarrow{\tau }\}$	functions
${\Gamma }_{1+}$	1	1	1	1	1	1	1	1	$x^2,y^2,z^2$
${\Gamma }_{x+}$	1	-1	1	-1	1	-1	1	-1	$R_x,yz$
${\Gamma }_{z+}$	1	1	-1	-1	1	1	-1	-1	$R_z,xy$
${\Gamma }_{y+}$	1	-1	-1	1	1	-1	-1	1	$R_y,zx$
${\Gamma }_{1-}$	1	1	1	1	-1	-1	-1	-1	$xyz$
${\Gamma }_{x-}$	1	-1	1	-1	-1	1	-1	1	$x$
${\Gamma }_{z-}$	1	1	-1	-1	-1	-1	1	1	$z$
${\Gamma }_{y-}$	1	-1	-1	1	-1	1	1	-1	$y$

Table 4. Table 4: Multiplication table for the irreps of the character table given in Table 3 .

	${\Gamma }_{1+}$	${\Gamma }_{x+}$	${\Gamma }_{z+}$	${\Gamma }_{y+}$	${\Gamma }_{1-}$	${\Gamma }_{x-}$	${\Gamma }_{z-}$	${\Gamma }_{y-}$
${\Gamma }_{1+}$	${\Gamma }_{1+}$	${\Gamma }_{x+}$	${\Gamma }_{z+}$	${\Gamma }_{y+}$	${\Gamma }_{1-}$	${\Gamma }_{x-}$	${\Gamma }_{z-}$	${\Gamma }_{y-}$
${\Gamma }_{x+}$	${\Gamma }_{x+}$	${\Gamma }_{1+}$	${\Gamma }_{y+}$	${\Gamma }_{z+}$	${\Gamma }_{x-}$	${\Gamma }_{1-}$	${\Gamma }_{y-}$	${\Gamma }_{z-}$
${\Gamma }_{z+}$	${\Gamma }_{z+}$	${\Gamma }_{y+}$	${\Gamma }_{1+}$	${\Gamma }_{x+}$	${\Gamma }_{z-}$	${\Gamma }_{y-}$	${\Gamma }_{1-}$	${\Gamma }_{x-}$
${\Gamma }_{y+}$	${\Gamma }_{y+}$	${\Gamma }_{z+}$	${\Gamma }_{x+}$	${\Gamma }_{1+}$	${\Gamma }_{y-}$	${\Gamma }_{z-}$	${\Gamma }_{x-}$	${\Gamma }_{1-}$
${\Gamma }_{1-}$	${\Gamma }_{1-}$	${\Gamma }_{x-}$	${\Gamma }_{z-}$	${\Gamma }_{y-}$	${\Gamma }_{1+}$	${\Gamma }_{x+}$	${\Gamma }_{z+}$	${\Gamma }_{y+}$
${\Gamma }_{x-}$	${\Gamma }_{x-}$	${\Gamma }_{1-}$	${\Gamma }_{y-}$	${\Gamma }_{z-}$	${\Gamma }_{x+}$	${\Gamma }_{1+}$	${\Gamma }_{y+}$	${\Gamma }_{z+}$
${\Gamma }_{z-}$	${\Gamma }_{z-}$	${\Gamma }_{y-}$	${\Gamma }_{1-}$	${\Gamma }_{x-}$	${\Gamma }_{z+}$	${\Gamma }_{y+}$	${\Gamma }_{1+}$	${\Gamma }_{x+}$
${\Gamma }_{y-}$	${\Gamma }_{y-}$	${\Gamma }_{z-}$	${\Gamma }_{x-}$	${\Gamma }_{1-}$	${\Gamma }_{y+}$	${\Gamma }_{z+}$	${\Gamma }_{x+}$	${\Gamma }_{1+}$

Equations

$$H_{\text{ph6}}=\left[\begin{array}[]{ccc}H_{c2}&H_{c2c1}&H_{c2v1}\\ H_{c2c1}^{\dagger}&H_{c1}&H_{c1v1}\\ H_{c2v1}^{\dagger}&H_{c1v1}^{\dagger}&H_{v1}\end{array}\right]\,,$$

$$H_a=H_{0,a}+H_{k,a}+H_{k2,a},$$

$$H_{0,a}=\left[\begin{array}[]{cc}E_{a}&0\\ 0&E_{a}\end{array}\right]\,,$$

$$H_{k,a}=\left[\begin{array}[]{cc}-iP_{a}k_{y}-\alpha_{a}k_{x}&0\\ 0&-iP_{a}k_{y}+\alpha_{a}k_{x}\end{array}\right]\,,$$

$$\vec{\Pi }=\frac{\hslash }{m_0}\vec{p}+\frac{{\hslash }^2}{4m_0^2{c}^2}\left[\vec{\sigma }\times \vec{\nabla }V(\vec{r})\right],$$

$$H_{k2,a}=\frac{\hbar^{2}}{2m_{0}}\left[\begin{array}[]{cc}A_{a}k_{x}^{2}+B_{a}k_{y}^{2}&0\\ 0&A_{a}k_{x}^{2}+B_{a}k_{y}^{2}\end{array}\right]\,,$$

$${\text{D}}_{x(y)}(\vec{k})=\sum _{c,v}{\left|⟨v,\vec{k}\left|{\Pi }_{x(y)}\right|c,\vec{k}⟩\right|}^2,$$

$$\text{D}(\theta ,E)$$

$$\times \delta \left\{E-\left[E_c(\vec{k})-E_v(\vec{k})\right]\right\},$$

$$g_n=g_0-i\frac{2m_0}{{\hslash }^2}\sum _{l\ne n}\frac{{\Pi }_x^{nl}{\Pi }_y^{ln}-{\Pi }_y^{nl}{\Pi }_x^{ln}}{E_n-E_l}$$

$$g_{v1}=g_{c1}=g_0-2\left(\frac{2m_0}{{\hslash }^2}\right)\left(\frac{P_{c1v1}{\alpha }_{c1v1}}{E_g}\right),$$

$$g_{v1}=g_0-2\left(\frac{2m_0}{{\hslash }^2}\right)\left(\frac{P_{c1v1}{\alpha }_{c1v1}}{E_g}\right),$$

$$g_{c1}=g_0-2\left(\frac{2m_0}{{\hslash }^2}\right)\left(\frac{P_{c1v1}{\alpha }_{c1v1}}{E_g}+\frac{P_{c2c1}{\alpha }_{c2c1}}{E_{c2}-E_g}\right)$$

$$g_{c2}=g_0-2\left(\frac{2m_0}{{\hslash }^2}\right)\left(\frac{P_{c2c1}{\alpha }_{c2c1}}{E_{c2}-E_g}\right),$$

$$E_{c1(v1)}^{\pm }(n)=E_{c1(v1)}+\left[M_{c1(v1)}\left(n+\frac{1}{2}\right)\pm \frac{g_0}{2}\right]{\mu }_{B}B,$$

$$\text{v}(r)=\frac{e^2}{8{\epsilon }_0{r}_0}\left[H_0\left(\frac{\epsilon }{r_0}r\right)-Y_0\left(\frac{\epsilon }{r_0}r\right)\right],$$

$$E_{c2}=⟨{\Gamma }_{x+}^{c2}\left|{\mathbf{H}}_0\right|{\Gamma }_{x+}^{c2}⟩,$$

$$E_{c1}=⟨{\Gamma }_{z-}^{c1}\left|{\mathbf{H}}_0\right|{\Gamma }_{z-}^{c1}⟩,$$

$$⟨{\Gamma }_{x+}^{v1}\left|{\mathbf{H}}_0\right|{\Gamma }_{x+}^{v1}⟩=0,$$

$$P_{c2c1}=i\frac{\hslash }{m_0}⟨{\Gamma }_{x+}^{c2}\left|p_y\right|{\Gamma }_{z-}^{c1}⟩,$$

$$P_{c1v1}=i\frac{\hslash }{m_0}⟨{\Gamma }_{z-}^{c1}\left|p_y\right|{\Gamma }_{x+}^{v1}⟩,$$

$${\alpha }_{c2c1}=\frac{{\hslash }^2}{4m_0^2{c}^2}⟨{\Gamma }_{x+}^{c2}\left|\frac{\partial V_C}{\partial y}\right|{\Gamma }_{z-}^{c1}⟩,$$

$${\alpha }_{c1v1}=\frac{{\hslash }^2}{4m_0^2{c}^2}⟨{\Gamma }_{z-}^{c1}\left|\frac{\partial V_C}{\partial y}\right|{\Gamma }_{x+}^{v1}⟩,$$

$$A_{c2}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{1-}\right]}\frac{{\left|⟨{\Gamma }_{x+}^{c2}\left|p_x\right|\beta ⟩\right|}^2}{E({\Gamma }_{x+}^{c2})-E_{\beta }},$$

$$B_{c2}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{z-}\right]}\frac{{\left|⟨{\Gamma }_{x+}^{c2}\left|p_y\right|\beta ⟩\right|}^2}{E({\Gamma }_{x+}^{c2})-E_{\beta }},$$

$$A_{c2v1}=\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{1-}\right]}\frac{⟨{\Gamma }_{x+}^{c2}\left|p_x\right|\beta ⟩⟨\beta \left|p_x\right|{\Gamma }_{x+}^{v1}⟩}{E({\Gamma }_{x+}^{c2},{\Gamma }_{x+}^{v1})-E_{\beta }},$$

$$B_{c2v1}=\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{z-}\right]}\frac{⟨{\Gamma }_{x+}^{c2}\left|p_y\right|\beta ⟩⟨\beta \left|p_y\right|{\Gamma }_{x+}^{v1}⟩}{E({\Gamma }_{x+}^{c2},{\Gamma }_{x+}^{v1})-E_{\beta }},$$

$$A_{c1}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{y+}\right]}\frac{{\left|⟨{\Gamma }_{z-}^{c1}\left|p_x\right|\beta ⟩\right|}^2}{E({\Gamma }_{z-}^{c1})-E_{\beta }},$$

$$B_{c1}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{x+}\right]}\frac{{\left|⟨{\Gamma }_{z-}^{c1}\left|p_y\right|\beta ⟩\right|}^2}{E({\Gamma }_{z-}^{c1})-E_{\beta }},$$

$$A_{v1}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{1-}\right]}\frac{{\left|⟨{\Gamma }_{x+}^{v1}\left|p_x\right|\beta ⟩\right|}^2}{E({\Gamma }_{x+}^{v1})-E_{\beta }},$$

$$B_{v1}=1+\frac{2}{m_0}\sum _{\beta }^{B\left[{\Gamma }_{z-}\right]}\frac{{\left|⟨{\Gamma }_{x+}^{v1}\left|p_y\right|\beta ⟩\right|}^2}{E({\Gamma }_{x+}^{v1})-E_{\beta }},$$

Full text

I Introduction

The two-dimensional (2D) phosphorene was first synthesized in 2014Liu2014ACSNano ; Lu2014NR ; Liang2014NL ; Zhang2014ACSNano ; CastellanosGomez20142DMat ; Qiao2014NatComm and showed remarkable physical properties. For example, phosphorene is a direct band gap semiconductor with enhanced photoluminescense intensity compared to bulk black phosphorusZhang2014ACSNano . Moreover, its strong coupling to light can be varied in the far-infrared to red spectral rangeZhang2014ACSNano due to high sensitivity of the band gap to the number of monolayersKeyes1953 ; Das2014NL ; Liang2014NL ; CastellanosGomez20142DMat ; Wang2015NatNano ; Frank2019PRX . Due to the puckered crystalline structure phsophorene shows strongly anisotropic electronic propertiesLiu2014ACSNano ; Qiao2014NatComm . High carrier mobility in phosphorene allows not only observation of fundamental quantum phenomena, such as Landau levelsLi2015NatNano ; Gillgren20152DMat or quantum Hall effectLi2016NatNano , but also its potential applications to semiconductor spintronicsZutic2004RMP ; Han2014NatNano , due to the weak spin-orbit coupling (SOC) of phosphorusKurpas2016PRB . Indeed, nanosecond spin lifetimes observed in all-electrical spin injection experiments and realization of a spin valve operating at room temperatureAvsar2017NatPhys already demonstrated the robustness of spin coherence in phosphorene. Furthermore, few-layer phosphorene heterostructures are promising candidates for ultrafast switching based on optically generated surface polaritonsHuber2017NatNano .

There have already been applications of the $k\cdot p$ method to describe phosphorene.Rodin2014PRL ; Li2014PRB ; LewYanVoon2015NJP ; PereiraJr2015PRB ; Kafaei2018JAP However, to also study spin-orbit effects the effective models should fully exploit the symmetry of phosphorene by capturing the anisotropy of the interband dipole coupling. Effective g-factors illustrate this best. For monolayer transition metal dichalcogenides (TMDCs) the g-factors can be derived employing the conventional $\vec{p}$ matrix elementsKormanyos2013PRB ; Kormanyos2014PRX ; Kormanyos2015NJP (equivalent to the so called Kane matrix element in zinc-blende structuresKane1957JPCS ), in the perturbative fashion within the $k\cdot p$ frameworkRoth1959PR ; Hermann1977PRB ; LewYanVoon2009 . In phosphorene, due to its two-fold symmetry embedded in the $D_{2h}$ group, the usual Kane-like matrix elements can account only for the coupling along the armchair direction. In order to include the contribution along the zigzag direction, one needs to go beyond the conventional $\vec{p}$ matrix element and include the $k$-dependent SOC contribution. The inclusion of this term was already shown by Zhou et al. [Zhou2017PRB_gph, ] to provide the g-factor correction in phosphorene thin films with its value estimated from experimental data. As another example, the inclusion of such $k$-dependent SOC terms in III-V wurtzite semiconductors was recently shown to provide a more reliable fitting to the ab initio spin splitting and overall band structureFariaJunior2016PRB and also to add sizable corrections to the total value of the effective g-factorsFariaJunior2019PRB . Furthermore, more complete $k\cdot p$ models for phosphorene can be used as building blocks to model van der Waals heterostructuresGeim2013Nat , for instance combined with TMDCsKormanyos20152DMat and (In,Ga)Se materialsZhou2017PRB_kpInSe , overcoming computational costs of ab initio calculations.

In this work, we investigate important physical features which appear due to the inclusion of the interband SOC term in effective $k\cdot p$ models for monolayer phosphorene. We show that this additional SOC term not only provides the interband dipole coupling along zigzag direction but also allows us to predict the values of the effective g-factors from a full theoretical perspective highlighting the intrinsic contributions of monolayer phosphorene. We also analyze the Landau level (LL) spectra and show that the interband SOC term provides a correction to Zeeman splitting, consistent with the g-factor approach. Finally, we combine our effective models with the Bethe-Salpeter equation (BSE) and show that the excitonic spectra of monolayer phosphorene are in agreement with the available data in the literature. We point out that our $k\cdot p$ parameters are obtained from a systematic fitting to ab initio band structure calculations taking into account the $k$-dependence of the energy bands and the effective masses (weighting the contribution around $\Gamma$-point), a crucial point to correctly describe the linear-in-k couplings in the Hamiltonian. Besides providing a simplified and tangible understanding of the underlying physics, these effective $k\cdot p$ models might be used to investigate additional properties in phosphorene but also of more complex systems such as van der Waals heterostructures composed of several layered materials.

The paper is organized as follows: In Sec. II we discuss the ab initio band structure of phosphorene. The effective $k\cdot p$ models with the inclusion of the interband SOC are addressed in Sec. III. In Sec. IV we investigate the behavior of the different $k\cdot p$ models under external magnetic field by calculating the effective g-factors and the LL spectra. Excitonic effects are presented in Sec. V and finally, in Sec. VI we draw our conclusions.

II Phosphorene band structure from ab initio

The initial crystal structure of phosphorene layer was taken from the bulk black phosphorusBrown1965 . New cell parameters were then found by fully relaxing a sheet of phosphorene using quasi-Newton variable-cell scheme as implemented in the Quantum Espresso packageQE2009 ; QE2017 . During this process all atoms were free to move in all directions in order to relax internal forces. The force convergence threshold and total energy convergence threshold were set to $10^{-4}$ Ry/a.u. and $10^{-5}$ Ry/a.u., respectively. A vacuum of 20$\leavevmode\nobreak\ \rm\AA$ was introduced in order to reduce spurious interactions between the periodic copies of the system. We used the norm-conserving pseudopotential with the Perdew-Burke-Ernzerhof (PBE)perdew_1996 version of the generalized gradient approximation (GGA) exchange-correlation potentials, with kinetic energy cutoffs of 70 Ry and 280 Ry for the wavefunction and charge density, respectively. The optimized cell parameters are a = 3.2986 Å  along the zigzag edge, and b = 4.6201 Å along the armchair edge. The electronic and spin properties of phosphorene were calculated using the full-potential augmented plane-wave all-electron code package WIEN2kwien2k . Self-consistency was achieved with 151 k-points in the irreducible wedge of the Brillouin zone. Tuning of the band gap was achieved by combining the undressed LDA functional with the modified Becke-Johnson potentialmBJ . The parameters of the latter were chosen to give the gap close to the recent experimentalLiang2014NL ; Wang2015NatNano and theoreticalFrank2019PRX values and to provide a realistic description of phosphorene.

The resulting ab initio band structure is shown in Fig. 1. The calculation reproduces the direct band gap of phosphorene with the value of 2.178 eV centered at the $\Gamma$ point. This value is close to the recent quantum Monte Carlo result of 2.4 eVFrank2019PRX . The energy dispersion of the top-most valence (labeled $v1$) and bottom-most conduction (labeled $c1$) bands close to the band edge displays a sizeable anisotropy with respect to the main crystallographic directions. This anisotropy is particularly pronounced for the band $v1$, which is almost flat for momenta along the $\Gamma$-X direction (along the zigzag edge of phosphorene), and very dispersive along the $\Gamma$-Y path (along the armchair edge do phosphorene). For this band, the corresponding effective masses $m_{v1,x}$ and $m_{v1,y}$ differ by more than an order of magnitude (see Section III).

Phosphorus is a light element, therefore the influence of spin-orbit coupling on the band structure is mainly limited to the removal of orbital degeneracies at certain high-symmetry points and lines, e.g., at the X-S and Y-S paths. At the S point, the splitting is 21 meV and 17 meV for the valence band $v1$ and conduction band $c1$, respectivelyKurpas2016PRB . The spin degeneracy of the energy bands is not removed by SOC, by virtue of space inversion and time reversal symmetry. We also identify in Fig. 1 the single group irreducible representations (irreps) for the energy bands at $\Gamma$ point that are used as input to the k.p Hamiltonians discussed in the next section.

III $k\cdot p$ modeling

For the effective description of the phosphorene band structure, given in Fig. 1, we focus on the first valence ($v1$), first conduction ($c1$) and second conduction ($c2$) bands, identified by the irreps of the symmetry group $D_{2h}$ provided by the ab initio calculations. This irrep identification of the energy bands is a crucial step for the development of $k\cdot p$ HamiltoniansDresselhaus1955PR ; Dresselhaus2008 ; FariaJunior2016PRB . The Hamiltonian basis set, including spin, is then given by $\left\{\left|\Gamma_{x+}^{c2}\uparrow\right\rangle,\left|\Gamma_{x+}^{c2}\downarrow\right\rangle,\left|\Gamma_{z-}^{c1}\uparrow\right\rangle,\left|\Gamma_{z-}^{c1}\downarrow\right\rangle,\left|\Gamma_{x+}^{v1}\uparrow\right\rangle,\left|\Gamma_{x+}^{v1}\downarrow\right\rangle\right\}$, in which the orbital part $\left|\Gamma_{\alpha}^{n}\right\rangle$ is written in Dirac notation, defined as $\left\langle\vec{r}|\Gamma_{\alpha}^{n}\right\rangle=u_{\alpha}^{n}(\vec{r})$. Essentially, our notation means that the Bloch function $u_{\alpha}^{n}(\vec{r})$ of the energy band $n$ transforms as the irrep $\Gamma_{\alpha}$ of the symmetry group $D_{2h}$. The vertical arrows represent the spin up and down projections, eigenvalues of the $\sigma_{z}$ Pauli matrix. In Sec. I of the Supplemental Materialsupmat we summarize the symmetry properties of the group $D_{2h}$ of phosphorene.

Within the $k\cdot p$ framework, we must compute matrix elements between different energy bands that are mediated by vector (such as $\vec{p}$) or pseudovector (such as $\vec{\nabla}V\times\vec{p}$) operators. In the language of group theory, we would write these matrix elements as the direct product $\Gamma^{n}\otimes\Gamma^{o}\otimes\Gamma^{m}$, with $\Gamma^{n(m)}$ representing the irrep of the energy band $n(m)$ and $\Gamma^{o}$ representing the irrep of the operator. If the result of such direct product contains the identity irrep (with characters equal to 1 for all symmetry operations), then the matrix element is nonzeroDresselhaus2008 . Furthermore, since the irreps of the symmetry group $D_{2h}$ are all one dimensional, the nonzero matrix elements are readily available by inspecting the multiplication table, given in Sec. I of the Supplemental Materialsupmat . We point out that a thorough and systematic investigation of the symmetry properties of phosphorene has been performed by Li and AppelbaumLi2014PRB , however without providing realistic values for the $k\cdot p$ parameters. In our study, we go beyond the derivation of $k\cdot p$ Hamiltonians based on symmetry arguments and also focus on the estimation of the Hamiltonian parameters by using the ab initio data presented in Sec. II.

Our strategy is to investigate different $k\cdot p$ formulations, including different energy bands and couplings, to demonstrate qualitative differences stemming from various couplings. In what follows we define the different $k\cdot p$ models we considered: ph6 (describes $c2$, $c1$, and $v1$ bands), ph4 (describes $c1$ and $v1$ bands), and Nph4 (also describes $c1$ and $v1$ bands but without any coupling between them). For the most general model, ph6, we can write the total Hamiltonian identifying the different coupling blocks

[TABLE]

with each coupling block being a $2\times 2$ matrix. To obtain the Hamiltonian for the ph4 model, $H_{\text{ph4}}$, we just remove the blocks that have contribution of the band $c2$, i. e., the 1st row and the 1st column of Eq. (1). For the model Nph4, we can obtain the Hamiltonian $H_{\text{Nph4}}$ by just removing the interaction block $H_{c1v1}$ between $c1$ and $v1$ bands from $H_{\text{ph4}}$.

Each of the Hamiltonian blocks in Eq. (1) can be written as

[TABLE]

with $H_{0,a}$ containing only $k$-independent terms, $H_{k,a}$ containing terms linear in $k$, $H_{k2,a}$ containing terms quadratic in $k$ with the subindex $a$ indicating the specific block. Specifically, taking into account the symmetry properties of phosphorene, the $k$-independent Hamiltonian block is given by

[TABLE]

which appears only for $a=\{c2,c1\}$ and the parameters $E_{a}$ indicate the energy values at $\Gamma$-point (note that $E_{c1}=E_{g}$). The Hamiltonian block with linear contribution of $k$ is given by

[TABLE]

which appears only for $a=\{c2c1,c1v1\}$ and the parameters $P_{a}$ and $\alpha_{a}$ originate from first-order perturbation theory of the operator

[TABLE]

with $P_{a}$ being the conventional Kane-like matrix elements [originated from the 1st term in Eq. (5)] and $\alpha_{a}$ being the $k$-dependent interband SOC parameter [originated from the 2nd term in Eq. (5)]. We note that due to the symmetry of energy bands involved, the $\alpha_{a}$ terms do not mix different spins. Finally, the Hamiltonian block with quadratic contribution of $k$ is given by

[TABLE]

which appears only for $a=\{c2,c1,v1,c2v1\}$ and the parameters $A_{a}$ and $B_{a}$ are the effective mass parameters obtained from the second-order perturbation theory. In the Supplemental Materialsupmat we provide the specific definitions of these parameters in Sec. II and in Sec. III we write the three different $k\cdot p$ models explicitly.

To obtain the values for the different parameters that appear in the $k\cdot p$ Hamiltonians we perform a numerical fitting to the ab initio data. For a reliable description, we considered the band structure and the effective masses (as function of the wavevector $k$) along multiple directions of the first Brillouin zone starting from the $\Gamma$-point, namely $\Gamma$-X (along $k_{x}$), $\Gamma$-Y (along $k_{y}$) and $\Gamma$-S (which combines $k_{x}$ and $k_{y}$). These different constraints (band structure, effective masses and multiple directions) are used simultaneously in the fitting procedure, implemented via the LMFIT package in Pythonlmfit ; FariaJunior2016PRB . The use of the $k$-dependent effective mass calculated from ab initio weights the points in the vicinity of the $\Gamma$-point. Furthermore, special care must be taken for the linear-in-k terms, specially the alpha parameters. From ab initio we calculate the dipole matrix elements without SOC (within the optic code in WIEN2kDraxl2006CPC ), found $P_{c1v1}$ = 5.3230 eV$\cdot\textrm{\AA}$ and $P_{c2c1}$ = 0.1757 eV$\cdot\textrm{\AA}$. For the $\alpha$ parameter, we can estimate its value by comparing the change in the effective mass with and without SOC within the ab initio calculations. This approach gives as a range of values of $\alpha_{c1v1}\sim 0.007-0.015$ eV$\cdot\textrm{\AA}$. A detailed description of the evaluation of $\alpha_{c1v1}$ can be found in in Sec. IV of the Supplemental Materialsupmat . We emphasize that for each of the three different $k\cdot p$ models (ph6, ph4 and Nph4) we performed an individual fitting of the band structure and effective masses with double precision resulting in the parameter sets shown in Table 1. From the Nph4 parameters we can readily obtain the effective masses: $m_{c1,x}=1/A_{c1}=1.15$, $m_{c1,y}=1/B_{c1}=0.24$, $m_{v1,x}=1/A_{v1}=7.29$ and $m_{v1,y}=1/B_{v1}=0.24$. For the fitted values of $P_{c1v1}$ and $P_{c2c1}$ we our values shown in Table 1 are relatively close to the calculated by ab initio, thus showing that our fitting scheme is quite reliable. We point out that this small discrepancy between the fitted and calculated values from ab initio is well known within the $k\cdot p$ framework that the parameters are slightly modified to accommodate the band structure features. We note that the signs of the $\alpha$ parameters given in Table 1 can be reversed without affecting the band structure.

The comparison between the fitted $k\cdot p$ models and the ab initio data is shown in Fig. 2. Up to 1.5 $\text{nm}^{-1}$, the limit in $k$-space used in our fitting, we find a good agreement for the band structure of all $k\cdot p$ models considered. For the effective mass analysis, we focus on $v1$ and $c1$ bands. We notice that Nph4 shows a constant dependence in all directions but agrees well very close to the $\Gamma$-point (as a parabolic dispersion should), ph4 also shows a nearly-constant dependence along $\Gamma$-X but now acquires a a dependence along $\Gamma$-Y and $\Gamma$-S directions due to the interband couplings and, finally, ph6 provides a more realistic description of the effective mass k-dependence, with additional contributions due to the quadratic coupling between to $c2$-$c1$ bands and the interband coupling terms. Furthermore, although we enforced the $k$-limit of 1.5 $\text{nm}^{-1}$ for the fitting, we also find a good agreement with ab initio at even larger $k$-values for the band structure, except for v1 band along $\Gamma$-X in ph4 and Nph4 models. Within the ph6 model it is even possible to model the anticrossing of the $c2$ and $c1$ conduction bands along $\Gamma$-Y.

Let us now turn to additional properties that can be derived from the $k\cdot p$ models. As a consequence of including the interband SOC term (in ph4 and ph6 models), it is possible to compute the dipole strength between $v1$ and $c1$ bands not only for the armchair but also for the zigzag direction. The dipole strength between the top-most valence and bottom-most conduction bands as function of the wavevector $\vec{k}$ can be written as

[TABLE]

in which $x(y)$ refers to the zigzag (armchair) direction, the summation for the indices $c(v)$ takes into account both spin components of $c1(v1)$ bands and the $\vec{\Pi}$ operator is given in Eq. (5). For the angular dependence of the dipole strength as a function of the transition energy, we can define the following quantity

[TABLE]

where $\theta$ is the angle defined as zero with respect to the armchair axis (y direction) and the $\vec{k}$ dependence in the summation take into account all possible transitions at a given transition energy of $E_{c}(\vec{k})-E_{v}(\vec{k})$. To compute such dipole strengths we use the values of $P$ and $\alpha$ given in Table 1.

We discuss in Fig. 3 the dipole strength features of the ph4 and ph6 $k\cdot p$ Hamiltonians. In Figs. 3(a-c) we show the calculated dipole strength ratio between armchair and zigzag directions [$\text{D}_{y}$ / $\text{D}_{x}$, obtained from Eq. (7)], as function of the transition energy, mapping the wavevector $\vec{k}$ to the transition energy $E_{c}(\vec{k})-E_{v}(\vec{k})$, along different directions of the Brillouin zone. We found that this ratio is quite large (4 to 5 orders of magnitude) but it is nonetheless nonzero, a feature that can only be achieved with the inclusion of the interband SOC terms (the $\alpha$ parameters). Moreover, this ratio for ph4 and ph6 are slighly different but both models (with their respective parameter sets) provide nearly the same trends. A complementary, and perhaps more instructive, way to visualize the dipole strength is to look at its angular dependence for a fixed transition energy using Eq. (8). For a transition energy of 5 meV above the band gap, we show in Fig. 3(d) the angular behavior of the dipole strength assuming the armchair direction at $\theta=0$. The 2-fold behavior of this angular dependence clearly reflects the symmetry of the phosphorene lattice. From the experimental perspective, recent photoluminescence measurements by Wang et al.Wang2015NatNano reports that the emission along zigzag is consistently less than 3% of the emission along armchair while in the study by Xu et al.Xu2016AdvMat , the photoluminescence intensity as function of the polarization angle [Fig. 3(c)] suggests that the ratio between armchair and zigzag directions is two orders of magnitude or more. Based on our findings of $\text{D}_{y}$ / $\text{D}_{x}$ which take into account the intrinsic selection rules in phosphorene, the reduced dipole ratio observed in experiment might be associated to more complex phenomena of many-body interactions and carrier relaxation.

It is also interesting to investigate the behavior of the effective $k\cdot p$ models by calculating the density of states (DOS) and carrier density, physical quantities that play an important role in transport experiments. We show in Figs. 4(a-b) the DOS for conduction and valence bands, respectively. As we noted previously for the band structure and effective masses shown in Fig. 2, the main differences among the $k\cdot p$ models appear in the valence band. Specifically, Nph4 provides a constant DOS dispersion, ph4 a nearly-constant behavior and ph6 shows a higher value close to the band edge which decreases slightly as the energy increases, also observed in the ab initio calculations in Fig.1(a) of Ref. [Rudenko2015PRB, ]. This is a typical feature observed in systems with reduced dimensionality and large effective mass, for instance in the valence band of conventional quantum wells based on zinc-blende GaAsFariaJunior2015PRB and wurtzite GaNChuang1997SST . By integrating the DOS we can obtain the carrier density as function of the energy, and similar trends of the DOS can also be seen in Figs. 4(c-d) for conduction and valence band, respectively.

IV Effective g-factors and Landau levels

In this section we incorporate the influence of external magnetic fields within the $k\cdot p$ models to investigate the effective g-factors and the LL spectra. Let us first start with the calculation of the effective g-factors, following the conventional perturbative approach within the $k\cdot p$ frameworkRoth1959PR ; Hermann1977PRB ; LewYanVoon2009 . We focus on magnetic fields that are oriented out of the monolayer plane, i. e., along $z$ direction, following the coordinate system shown in the inset of Fig. 1. Under these conditions, the effective g-factor can be generally written as

[TABLE]

in which $g_{0}$ is the bare electron g-factor, $n$ is the band of interest, $l$ runs over the other bands in the $k\cdot p$ model, the energy values in the denominator are the values at $\Gamma$-point, the $\vec{\Pi}$ operator is defined in Eq. 5 and $\Pi^{nl}_{x(y)}=\left\langle n\left|\Pi_{x(y)}\right|l\right\rangle$, which follows the same form of the interband dipole coupling shown in Eq. (7).

Evaluating Eq. (9) specifically for the ph4 $k\cdot p$ model, the top-most valence and bottom-most conduction band g-factors read

[TABLE]

while for the ph6 $k\cdot p$ model, the g-factors are given by

[TABLE]

[TABLE]

and

[TABLE]

with the parameters $P$ and $\alpha$ (discussed in Sec. III) with values given in Table 1 were used to compute the g-factor values. We emphasize that without including the interband SOC term (given by the $\alpha$ parameters) there is no correction to the effective g-factors from the bare electron g-factor, $g_{0}$.

Evaluating the g-factors given in Eqs. (10), (11), (12) and (13) using the parameters given in Table 1, we show in Table 2 our predicted values for the effective g-factors within the ph4 and ph6 models. Due to the small value of $\alpha$, the corrections to the bare electron g-factor are of the order of $10^{-2}$. Furthermore, the values obtained for $g_{v1}$ and $g_{c1}$ from the different $k\cdot p$ models are consistent with each other. For the second conduction band, $g_{c2}$ can be only be accounted within the ph6 model and our predicted value is nearly the bare electron g-factor. The available g-factors experimentally determined for few-layer phosphorene via transport experiments for the top-most valence band are $g_{v1}=$ 1.8 - 2.7 by Gillgren et al.Gillgren20152DMat and $g_{v1}=2.0\pm 0.1$ by Li et al.Li2016NatNano . From a theoretical perspective, the value of $g_{v1}=2.14$ for monolayer black phosphorus has been determined by Zhou et al.Zhou2017PRB_gph , however, the interband SOC parameter $\alpha_{c1v1}$ was not theoretically obtained but estimated from the experimental data of Ref. [Long20162DMat, ], leading to a value of $\alpha_{c1v1}\sim 0.45\;\textrm{eV}\cdot\textrm{\AA}$. For the conduction band, recent experiments by Yang et al.Yang2018NL in few-layer black phosphorus found surprisingly large g-factors of 5.7$\pm$0.7 for filling factor $\nu=1$ due to electron-electron interaction. The g-factor value drops to 2.5$\pm$0.1 for filling factor $\nu=7$ once screening effects take place. Given the uncertainty in the experimentally determined effective g-factors obtained in few-layer phosphorene samples, certainly with the influence of many-body interactions, our predicted values have the advantage that they were obtained from a full theoretical approach taking into account the intrinsic contributions of bare phosphorene monolayers.

Now let us turn to the LL spectra of the system, going beyond the effective g-factor approach by considering the envelope function approximation, combined with the minimal coupling and the Zeeman termLewYanVoon2009 ; Pryor2006PRL ; vanBree2012PRB ; FariaJunior2019PRB . First, we notice that because of the lack of interband coupling in the Nph4 model, it is possible to find analytical solutions for the LLsZhou2015SREP , given by

[TABLE]

with the subindex $c1(v1)$ denoting the bottom-most conduction (top-most valence) band, the superindex $\pm$ indicating the positive and negative Zeeman split LLs, $n=0,1,2,\ldots$ indicating the LL index, $E_{c1}=E_{g}$, $M_{c1}=2\sqrt{A_{c1}B_{c1}}$, $E_{v1}=0$ and $M_{v1}=-2\sqrt{A_{v1}B_{v1}}$. For the ph4 and ph6 models, we employ the numerical technique of the finite differencesChuang1997SST to obtain the LL spectra. We also apply this numerical approach to the Nph4 model for comparison. We assumed $B=B\hat{z}$ and the vector potential given by $\vec{A}=Bx\hat{y}$. For the numerical discretization we considered the a system size of $L=200\;\text{nm}$ with 401 points and hard-wall boundary conditions.

In Fig. 5 we summarize the LL spectra and the ZS for the different $k\cdot p$ models with parameters from Table 1. Let us start with the Nph4 model, with LLs shown in Fig. 5(a) and ZSs in Figs. 5(d-e). We found a very good agreement between our numerical calculations and the analytical approach of Eq. (14), however, due to the lack of interband coupling, the ZS is given only by the bare electron g-factor. For ph4 and ph6 models LL spectra, shown in Figs. 5(b-c) respectively, the situation is quite similar with just different energy separations between LLs arising due to the different coupling in the Hamiltonians. Turning to the ZS, shown in Figs. 5(f-i), we observe that both ph4 and ph6 models slightly deviate from the bare electron g-factor case, as expected and already noticed in the effective g-factor case. These different features in the LLs and ZS available in the ph4 and ph6 models could be investigated as signatures in magneto-transport spectraZhou2015SREP ; Tahir2015PRB .

V Excitons

Since phosphorene is a direct band gap semiconductor with interesting optical properties, it is also important to check the consistency of our effective $k\cdot p$ models and parameters by calculating the excitonic spectra. Using the effective BSERohlfing2000PRB ; Scharf2016PRB ; Tedeschi2018PRBr , we focus on direct excitons at zero temperature assuming the electron-hole interaction to be mediated by the Rytova-Keldysh potentialRytova1967 ; Keldysh1979 ; Cudazzo2011PRB , given by

[TABLE]

in which $H_{0}$ is the zeroth-order Struve function, $Y_{0}$ is the zeroth-order Bessel function of the second kind, $r=\sqrt{x^{2}+y^{2}}$, $e$ is the electron charge, $\varepsilon_{0}$ is the vacuum permittivity, $\varepsilon=(\varepsilon_{t}+\varepsilon_{b})/2$ is the effective dielectric constant given by the average of top, $\varepsilon_{t}$, and bottom, $\varepsilon_{b}$, dielectric constants, and $r_{0}=2\pi\zeta$ with $\zeta$ being the 2D dielectric susceptibility with typical values in the literature of $\zeta=4.1\;\textrm{\AA}$Rodin2014PRB , $\zeta=3.85\;\textrm{\AA}$Seixas2015PRB and $\zeta=3.8\;\textrm{\AA}$Prada2015PRB . In our calculations we assume $\zeta=4\;\textrm{\AA}$, leading to a value of $r_{0}\approx 25\;\textrm{\AA}$.

We analyze the effect of the dielectric environment considering different values of $\varepsilon$ from 1 to 5, in order to cover a reasonable range of experimental realizations of phosphorene. For instance, we identify three important cases that fall within the range of $\varepsilon$ values we considered: (i) freestanding phosphorene is equivalent to $\varepsilon=1$ ($\varepsilon_{t}=\varepsilon_{b}=1$); (ii) phosphorene on SiO${}_{\text{2}}$ substrate is equivalent to $\varepsilon=2.45$ ($\varepsilon_{t}=1$ and $\varepsilon_{b}=3.9$Berkelbach2013PRB ); and (iii) encapsulated phosphorene with boron nitride is equivalent to $\varepsilon=4.5$ ($\varepsilon_{t}=\varepsilon_{b}=4.5$Stier2018PRL ). For all the $k\cdot p$ models we solve the BSE numerically using a 2D $k$-grid of -0.5 to 0.5 $\textrm{\AA}^{-1}$ in $k_{x}$ (zigzag) and -0.3 to 0.3 $\textrm{\AA}^{-1}$ in $k_{y}$ (armchair) with a total discretization of $(2\text{N}_{x}+1)\times(2\text{N}_{y}+1)$ with $(\text{N}_{x},\text{N}_{y})=\left[(70,42),(60,36)\right]$ The relative error between these two mesh discretizations is $\sim 1\%$. The final values of exciton biding energies and diameters are then obtained using a linear extrapolation of the values calculated in these two grid sizes. We define the exciton diameter as the full width at half maximum of the exciton probability density in real space.

In Fig. 6 we summarize our findings for the excitonic spectra focusing on the lowest 3 excitonic states. We show in Figs. 6(a)-(c) the exciton biding energy, $\text{E}_{\text{B}}$ (measured from the single-particle band gap), the exciton diameter along the zigzag direction, $\text{L}_{x}$, and the exciton diameter along the armchair direction, $\text{L}_{y}$, as function of the effective dielectric constant, respectively. As an initial remark, all the models are in reasonable quantitative agreement with each other and thus it suffices to discuss the general behavior of $\text{E}_{\text{B}}$, $\text{L}_{x}$ and $\text{L}_{y}$ as function of $\varepsilon$ irrespective of the model. Increasing the value of $\varepsilon$, we show that the magnitude of $\text{E}_{\text{B}}$ decreases in a nonlinear fashion whereas $\text{L}_{x}$ and $\text{L}_{y}$ increase linearly, with increasing slope for higher exciton levels. This behavior is in qualitative agreement with 2D excitons obtained from the Wannier equation and the typical Coulomb potential ($\text{E}_{\text{B}}\sim\varepsilon^{-2}$ and $\text{R}\sim\varepsilon$)Haug2004 . In Figs. 6(d)-(f) and Figs. 6(g)-(i) we show the probability density of the exciton wavefunction in real space for $\varepsilon=1$ and $\varepsilon=5$, respectively. We can identify the first exciton to a 1s state, the second exciton to a 2p${}_{\text{x}}$ and the third exciton to a 2s state. This exciton ordering is in agreement with recent ab initio calculationsQiu2017NL . We emphasize that although Nph4 provides a consistent description of excitonic effects, it does not allow any further calculations involving dipole coupling [see Fig. 3] and therefore one should rely on either the ph4 or the ph6 models. From the exciton perspective we revisit the dipole ratio between armchair and zigzag directions for the first exciton level and found the values of $6.4\times 10^{4}$ for ph4 and $5.1\times 10^{4}$ for ph6, thus showing the same order of magnitude as the values calculated in Fig. 3.

Finally, let us compare our calculated exciton biding energies with the literature. For a freestanding mononalyer phosphorene we find $\text{E}_{\text{B}}$ $\sim$0.81 eV, in good agreement with the range found in the literature 0.75 - 0.86 eV Chaves2016PRB ; Prada2015PRB ; Rodin2014PRB ; Tran2014PRB ; Tran20152DMat ; Hunt2018PRB ; Frank2019PRX . For a monolayer phosphorene on SiO${}_{\text{2}}$ substrate we obtained $\text{E}_{\text{B}}$ $\sim$0.44 eV, consistent with the reported values of 0.38 eVCastellanosGomez20142DMat , 0.4 eVRodin2014PRB and 0.46 eVHunt2018PRB . And for encapsulated phosphorene with boron nitride, we found $\text{E}_{\text{B}}$ $\sim$0.26 eV, also in good agreement with reported values of $\sim$0.22 eVRodin2014PRB and $\sim$0.26 eVHunt2018PRB .

VI Conclusions

In summary, we have developed effective four- and six-band $k\cdot p$ Hamiltonians exploiting the full symmetry of monolayer phosphorene by including the interband spin-orbit coupling interaction, a term previously neglected in the literature and of crucial importance for the treatment of anisotropic 2D materials. Specifically, the inclusion of such spin-orbit term allows the calculation of two important features in phosphorene: (i) the proper description of the interband dipole interaction not only along armchair but also along zigzag directions and (ii) the estimation of the effective g-factors and the Zeeman splittings of Landau levels. To obtain the $k\cdot p$ parameters, we have performed a systematic fitting of the band structure and the $k$-dependence of the effective masses to reliable ab initio calculations. Our $k\cdot p$ models highlight the intrinsic characteristics of monolayer phosphorene and were investigated in light of different physical aspects showing a correct description of the dipole selection rules, effective g-factors, stable behavior of the Landau level spectra and consistent description of the excitonic spectra (in good agreement with reported values). Furthermore, our findings suggest that one must be careful in comparing the calculated intrinsic g-factors to the available experimental data since many-body effects could have an important contribution. Finally, the presented $k\cdot p$ Hamiltonians and parameters can be directly applied to investigate many-body effects, transport phenomena in the presence of magnetic field, optical properties including excitonic effects and can be possibly coupled to other available Hamiltonians of 2D materials to investigate novel van der Waals heterostructures.

Note added in proof. After this manuscript was completed we became aware of Ref. [Li2019arXiv, ] which also calculates the interband spin-orbit coupling parameter $\alpha$ using a different theoretical approach.

Acknowledgments

The authors acknowledge financial support from the Alexander von Humboldt Foundation, Capes (grant No. 99999.000420/2016-06), DFG SFB 1177 (A09, B05), National Science Centre under the contract DEC-2018/29/B/ST3/01892 and VVGS-2019-1227. We acknowledge anonymous referees for valuable feedback to this manuscript. We thank Pengke Li for helpful discussions.

Supplemental Material for the paper “$k\cdot p$ theory for phosphorene: effective g-factors, Landau levels, and excitons”

I. Symmetry properties of phosphorene

In this section we briefly discuss the symmetry properties of phosphorene. The unit cell in direct space and the coordinate system is shown in Fig. 7. In Table 3 we present the character table for the nonsymmorphic group $D_{2h}$ of phosphorene. The symmetry operations are as follows: $E$ is the identity, $C_{a}$ is a rotation of $\pi/2$ along the $a=x,y,z$ axis, $I$ is the inverstion and $\tau$ is the translation vector given by $\vec{\tau}=\frac{a}{2}\hat{x}+\frac{b}{2}\hat{y}$. The multiplication table of the irreducible represenations are given in the Table 4. Please refer to the study of Li and AppelbaumLi2014PRB …

II. Definition of k.p parameters

In this section we present the definitions of the k.p parameters that appear in the Hamiltonians in the main text. The unperturbed terms are given by

[TABLE]

[TABLE]

[TABLE]

in which $\mathbf{H_{0}}$ is the unperturbed term of the Hamiltonian.

For the terms that appear linear in k, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

in which $V_{C}$ is the crystal potential.

For the terms that appear with the quadratic k, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with the irreductible representations that contribute to nonzero matrix elements indicated in the square brackets above the summation.

III. Hamiltonians

Explicity matrix form of the Hamiltonian ph6:

[TABLE]

[TABLE]

Explicity matrix form of the Hamiltonian ph4:

[TABLE]

[TABLE]

Explicity matrix form of the Hamiltonian Nph4:

[TABLE]

[TABLE]

To shorten the notation, the $M_{a}$ terms are defined as

[TABLE]

IV. Estimating the interband spin-orbit coupling

The contribution to the effective mass in conduction bands from $\alpha_{c1v1}$ in the ph4 model can be written as:

[TABLE]

The value $\frac{\alpha_{c1v1}^{2}}{\left(\frac{\hbar^{2}}{2m_{0}}\right)E_{g}}$ is of the order of the spin-orbit couping (SOC) correction to the effective mass since no other parameters in the Hamiltonian would contribute to that. Therefore, by knowing the effective mass with and without SOC, we can estimate $\alpha_{c1v1}$ as:

[TABLE]

[TABLE]

Fitting of $1/m^{*}$ using the ab initio calculations without and with SOC in double precision:

[TABLE]