Limited accuracy of conduction band effective mass equations for semiconductor quantum dots

TL;DR

This paper systematically derives and evaluates conduction-band effective mass equations for semiconductor quantum dots, highlighting their limitations and the importance of including non-parabolicity and off-diagonal terms for accurate modeling.

Contribution

It provides a hierarchical classification of effective mass equations derived from 8-band kp theory and assesses their accuracy in quantum dot spectral and spin properties.

Findings

Effective mass equations have limited accuracy without non-parabolicity corrections.

Including off-diagonal valence band terms improves modeling accuracy.

The Roth-Lax-Zwerdling formula effectively predicts g-factors in nanostructures.

Abstract

Effective mass equations are the simplest models of carrier states in a semiconductor structures that reduce the complexity of a solid-state system to Schr\"odinger- or Pauli-like equations resempling those well known from quantum mechanics textbooks. Here we present a systematic derivation of a conduction-band effective mass equation for a self-assembled semiconductor quantum dot in a magnetic field from the 8-band kp theory. The derivation allows us to classify various forms of the effective mass equations in terms of a hierarchy of approximations. We assess the accuracy of the approximations in calculating selected spectral and spin-related characteristics. We indicate the importance of preserving the off-diagonal terms of the valence band Hamiltonian and argue that an effective mass theory cannot reach satisfactory accuracy without self-consistently including non-parabolicity corrections and renormalization of kp parameters. Quantitative comparison with the 8-band kp results supports the phenomenological Roth-Lax-Zwerdling formula for the g-factor in a nanostructure.

Tables

Table 1. Table 1: Material parameters used in the calculations. [ 33 , 3 ]

	GaAs	InAs	Interpolation for ${\mathrm{In}}_x{\mathrm{Ga}}_{1-x}\mathrm{As}$
$E_v$	$0.0$ eV	$0.21$ eV	linear
$E_0$	1.519 eV	0.417 eV	$0.417x+1.519\left(1-x\right)-0.477x\left(1-x\right)$
$E_0^{\prime }$	4.488 eV	4.390 eV	linear
$E_Q$	17.535 eV	18.255 eV	linear
$m^{\ast }$	$0.0665m_0$	$0.0229m_0$	$\left[0.0229x+0.0665\left(1-x\right)-0.0091x\left(1-x\right)\right]m_0$
$\Delta$	0.341 eV	0.39 eV	$0.39x+0.341\left(1-x\right)-0.15x\left(1-x\right)$
${\Delta }_0^{\prime }$	0.171 eV	0.25 eV	linear
$a_c$	-7.17 eV	-5.08 eV	$-5.08x-7.17\left(1-x\right)-2.61x\left(1-x\right)$
$a_v$	1.16 eV	1.00 eV	linear
$b_v$	-2.0 eV	-1.8 eV	linear
$d_v$	-4.8 eV	-3.6 eV	linear
${\gamma }_1$	6.98	20.0	$1/\left[\left(1-x\right)/6.98+x/20.0\right]$
${\gamma }_2$	2.06	8.5	$1/\left[\left(1-x\right)/8.5+x/2.06\right]$
${\gamma }_3$	2.93	9.2	$1/\left[\left(1-x\right)/9.2+x/2.93\right]$
$P^{\prime }$	$4.780i$	$0.873i$	linear
$Q$	8.165	8.331	linear

Equations

$$H=\left(\begin{array}[]{ccc}H_{\mathrm{6c6c}}&H_{\mathrm{6c8v}}&H_{\mathrm{6c7v}}\\ H_{\mathrm{8v6c}}&H_{\mathrm{8v8v}}&H_{\mathrm{8v7v}}\\ H_{\mathrm{7v6c}}&H_{\mathrm{7v8v}}&H_{\mathrm{7v7v}}\\ \end{array}\right),$$

$$H_{6c6c}$$

$$H_{8v8v}$$

$$+a_{\mathrm{v}}\mathrm{Tr}\eta -b_{\mathrm{v}}\left[\left(J_x^2-\frac{1}{3}{J}^2\right){\eta }_{xx}+c.p.\right]-\frac{d_{\mathrm{v}}}{\sqrt{3}}\left[\{J_x,J_y\}{\eta }_{xy}+c.p.\right]-i\frac{{\hslash }^2}{m_0}\left[k_{[x}{\kappa }^{\prime }{k}_{y]}{J}_z+k_{[x}q{k}_{y]}{J}_z^3+c.p.\right],$$

$$H_{7v7v}$$

$$H_{6c8v}$$

$$+\frac{\sqrt{3}}{2}\left[(T_{xx}-T_{yy})\left(\frac{2}{3}{k}_z{B}_{8v}^{-}{k}_z-\frac{1}{3}{k}_x{B}_{8v}^{-}{k}_x-\frac{1}{3}{k}_y{B}_{8v}^{-}{k}_y\right)-T_{zz}(k_x{B}_{8v}^{-}{k}_x-k_y{B}_{8v}^{-}{k}_y)\right],$$

$$H_{6c7v}$$

$$H_{8v7v}$$

$$-3b_{\mathrm{v}}\left(T_{xx}^{†}{\eta }_{xx}+c.p.\right)-\sqrt{3}{d}_{\mathrm{v}}\left(2T_{xy}^{†}{\eta }_{xy}+c.p.\right)-i\frac{3{\hslash }^2}{2m_0}\left[k_{[x}{\kappa }^{\prime }{k}_{y]}{T}_z^{†}+c.p.\right]-3\left({\mu }_B{B}_z{T}_z^{†}+c.p.\right).$$

$$A_c^{\prime }\equiv \frac{m_0}{m^{\prime }}=\frac{m_0}{m^{\ast }}-\frac{2}{3}\frac{E_P}{E_0}-\frac{1}{3}\frac{E_P}{E_0+{\Delta }_0},g^{\prime }=2,{\kappa }^{\prime }=-\frac{1}{3}\left({\gamma }_1^{\prime }-2{\gamma }_2^{\prime }-3{\gamma }_3^{\prime }+2\right);$$

$${\tilde{H}}^{(\mathrm{d})}=\cosh \hat{S}{H}^{(\mathrm{d})}+\sinh \hat{S}{H}^{(\mathrm{od})},{\tilde{H}}^{(\mathrm{od})}=\sinh \hat{S}{H}^{(\mathrm{d})}+\cosh \hat{S}{H}^{(\mathrm{od})}.$$

$$H^{(\mathrm{od})}=-\sinh \hat{S}{\tilde{H}}^{(\mathrm{d})}.$$

$$H^{(\mathrm{od})}=-\hat{S}{\tilde{H}}^{(\mathrm{d})}=[S,{\tilde{H}}^{(\mathrm{d})}].$$

$${\tilde{H}}_{\mathrm{c}}^{\prime }{S}_{\mathrm{cv}}-S_{\mathrm{cv}}({\tilde{H}}_{\mathrm{v}}^{\prime }-\hat{\Delta })=-H_{\mathrm{cv}}.$$

$$S_{\mathrm{cv}}=X^{(0)}+\sum _j{k}_j{X}_j^{(1)}+\sum _{jl}{k}_j{X}_{jl}^{(2)}{k}_l,H_{\mathrm{cv}}=-\sum _j{k}_j{C}_j^{(1)}-\sum _{jl}{k}_j{C}_{jl}^{(2)}{k}_l,$$

$${\tilde{H}}_{\mathrm{c}}^{\prime }{X}_{jl}^{(2)}-X_{jl}^{(2)}({\tilde{H}}_{\mathrm{v}}^{\prime }-\hat{\Delta })$$

$${\tilde{H}}_{\mathrm{c}}^{\prime }{X}_j^{(1)}-X_j^{(1)}({\tilde{H}}_{\mathrm{v}}^{\prime }-\hat{\Delta })$$

$${\tilde{H}}_{\mathrm{c}}^{\prime }{X}^{(0)}-X^{(0)}({\tilde{H}}_{\mathrm{v}}^{\prime }-\hat{\Delta })$$

$$\frac{m_0}{\tilde{m}}=\frac{m_0}{m^{\prime }}+\frac{2}{3}\frac{E_P}{E_g^{\prime }}+\frac{1}{3}\frac{E_P}{E_g^{\prime }+{\Delta }_0},{\tilde{\gamma }}_1={\gamma }_1^{\prime }+\frac{E_P}{3E_g^{\prime }+{\Delta }_0},{\tilde{\gamma }}_{2,3}={\gamma }_{2,3}^{\prime }+\frac{1}{2}\frac{E_P}{3E_g^{\prime }+{\Delta }_0},$$

$$X_{jl}^{(2)}=C_{jl}^{(2)}{\mathcal{D}}^{-1},X_j^{(1)}=C_j^{(1)}{\mathcal{D}}^{-1}+C_{jl}^{(2)}\left[k_l,{\mathcal{D}}^{-1}\right],X^{(0)}=C_j^{(1)}[k_j,{\chi }_{\mathrm{c}}^{\prime }]{\mathcal{D}}^{-2}-\left[k_l,[k_j,{\chi }_{\mathrm{c}}^{\prime }]C_{jl}^{(2)}\right]{\mathcal{D}}^{-2}.$$

$${\tilde{H}}^{(2)}=\sum _{jl}{k}_j{C}_j^{(1)}{\mathcal{D}}^{-1}{C}_l^{(1)†}{k}_l+\frac{1}{2}\sum _{jl}\left([k_j,{\chi }_{\mathrm{c}}^{\prime }]C_j^{(1)}{\mathcal{D}}^{-2}{C}_l^{(1)†}{k}_l+h.c.\right)$$

$${\tilde{H}}^{(3)}=\sum _{jln}{k}_j{C}_{jl}^{(2)}{k}_l{\mathcal{D}}^{-1}{C}_n^{(1)†}{k}_l+\frac{1}{2}\sum _{jln}\left\{-\left[k_l,[k_j,{\chi }_{\mathrm{c}}^{\prime }]C_{jl}^{(2)}\right]{\mathcal{D}}^{-2}{C}_n^{(1)†}{k}_n+k_l{C}_{jl}^{(2)}{k}_j{\mathcal{D}}^{-2}{C}_n^{(1)†}[k_n,{\chi }_{\mathrm{c}}^{\prime }]+h.c.\right\}.$$

$${\tilde{H}}^{(2)}=\sum _{jl}{k}_{j}P{\tilde{\mathcal{T}}}_j{\mathcal{D}}^{-1}{\tilde{\mathcal{T}}}_l^{†}P{k}_l+\frac{1}{2}\sum _{jl}\left(P{\tilde{\mathcal{T}}}_j{\mathcal{D}}^{-1}[k_j,{\chi }_{\mathrm{c}}^{\prime }]{\mathcal{D}}^{-1}{\tilde{\mathcal{T}}}_l^{†}P{k}_l+h.c.\right),$$

$${\tilde{H}}^{(2)}=\frac{1}{2}\sum _{jl}\left\{k_j{k}_l,P{\tilde{\mathcal{T}}}_j{\mathcal{D}}^{-1}{\tilde{\mathcal{T}}}_l^{†}P\right\}+\frac{1}{2}\sum _{jl}\left(P{\tilde{\mathcal{T}}}_j{\mathcal{D}}^{-1}[k_j,H_{\mathrm{v}}^{\prime }]{\mathcal{D}}^{-1}{\tilde{\mathcal{T}}}_l^{†}P{k}_l+h.c.\right),$$

$${\left({\tilde{\mathcal{T}}}_j{\tilde{\mathcal{D}}}^{-n}{\tilde{\mathcal{T}}}_l^{†}\right)}_{\mathrm{s},\mathrm{as}}=\frac{1}{2}\left({\tilde{\mathcal{T}}}_j{\tilde{\mathcal{D}}}^{-n}{\tilde{\mathcal{T}}}_l^{†}\pm {\tilde{\mathcal{T}}}_l{\tilde{\mathcal{D}}}^{-n}{\tilde{\mathcal{T}}}_j^{†}\right).$$

$${\left({\tilde{\mathcal{T}}}_j{\tilde{\mathcal{D}}}^{-n}{\tilde{\mathcal{T}}}_l^{†}\right)}_{\mathrm{s}}=\frac{{\hslash }^2}{2m_0{P}^2}\left({\delta }_{jl}-2{\eta }_{jl}\right)\frac{F_j^{(n)}+F_l^{(n)}}{2}{\mathbb{I}}_{2\times 2},$$

$${\tilde{H}}^{(\mathrm{s},1)}=\sum _j\frac{{\hslash }^2{k}_j^2}{2m_j^{\ast }},\text{mbox}where\frac{m_0}{m_j^{\ast }}=\frac{m_0}{m^{\prime }}+F_j^{(1)}.$$

$${\tilde{H}}^{(\mathrm{s},2)}=-\frac{{\hslash }^2}{4m_0}\sum _{jl}\left[k_l,[k_j,{\chi }_{\mathrm{c}}^{\prime }]\left(\frac{1}{2}{\delta }_{jl}-{\eta }_{jl}\right)\left(F_j^{(2)}+F_l^{(2)}\right)\right].$$

$${\tilde{H}}^{(\mathrm{as},1)}=\frac{1}{4}\sum _{jl}\left\{k_j,\left[k_l,{\mathcal{A}}_{jl}\right]\right\},\text{mbox}where{\mathcal{A}}_{jl}=\frac{i{\hslash }^2{g}^{\prime }}{4m_0}\sum _m{\epsilon }_{jlm}{\sigma }_m+P{\left({\tilde{\mathcal{T}}}_j{\mathcal{D}}^{-1}{\tilde{\mathcal{T}}}_l^{†}\right)}_{\mathrm{as}}P.$$

$${\left({\tilde{\mathcal{T}}}_j{\tilde{\mathcal{D}}}^{-n}{\tilde{\mathcal{T}}}_l^{†}\right)}_{\mathrm{as}}=-\frac{i}{2E_{\mathrm{P}}}\sum _m\left[{\epsilon }_{jlm}-\sum _n({\eta }_{jn}{\epsilon }_{nlm}+{\eta }_{nl}{\epsilon }_{jnm})\right]\delta g_m^{(n)}{\sigma }_m,$$

$$\delta g_x^{(n)}=\delta g_y^{(n)}=(2/3)E_{\mathrm{P}}(1/E_{\mathrm{lh}}^n-1/E_{\mathrm{so}}^n),\delta g_z^{(n)}=E_{\mathrm{P}}[-1/(3E_{\mathrm{lh}}^n)+1/E_{\mathrm{hh}}^n-2/(3E_{\mathrm{so}}^n)].$$

$${\tilde{H}}^{(\mathrm{as},1)}\approx \frac{1}{2}{\mu }_{\mathrm{B}}\mathbf{B}\hat{g}\mathbf{\sigma }.$$

$$g^{(\mathrm{Roth})}=g^{\prime }-\frac{2E_{\mathrm{P}}}{3}\left(\frac{1}{E_{\mathrm{g}}}-\frac{1}{E_{\mathrm{g}}+{\Delta }_{\mathrm{SO}}}\right).$$

$$H^{(\mathrm{as},2)}=\frac{e{P}^2}{4E_{\mathrm{P}}}\sum _{jlm}{\mathcal{E}}_j{\epsilon }_{jlm}\delta g_m^{(2)}{\sigma }_{\mathrm{m}}{k}_l+h.c.$$

Full text

Introduction

The $\bm{k}\!\cdot\!\bm{p}$ method is a well-established approach to calculating the electronic and magnetic properties of bulk semiconductors [1]. It has also been applied to semiconductor nanostructures within the envelope function approximation [2, 3], in which a carrier state is described as a superposition of contributions from different bands, with local amplitudes smoothly varying in space and referred to as the envelope wave functions. Its accuracy is in many cases sufficient in many applications throughout mesoscopic physics, even when subtle details of the spectrum are considered [4, 5, 6]. At the same time, it offers much lower computiational costs, better scalability to larger systems and much higher transparency than the more exact atomistic approaches [7, 8, 9]. The applicability of the $\bm{k}\!\cdot\!\bm{p}$ method for modulated systems (nanostructures) has been confirmed by providing a rigorous derivation from the full Schrödinger equation [10, 11]. Recent developments bring the $\bm{k}\!\cdot\!\bm{p}$ methods down to the atomistic level[12]. Currently, the 8-band $\bm{k}\!\cdot\!\bm{p}$ model is a widely tested and generally trusted standard for calculating carrier states (including magnetic effects) in mesoscopic structures, in particular in strained, self-assembled systems[13, 14, 15, 16].

In many cases, when the focus is on the electron states, a single-band description is desired, which would offer further reduction of computational costs, make it possible to approximate the problem by a simple, analytically solvable model, and provide still more straight-forward interpretation referring to the well-known properties of the usual Schrödinger or Pauli equations. Such models, known as effective mass equations, have been used for decades to describe smoothly varying perturbations of semiconductor systems, like shallow impurities [17] and later to semiconductor heterostructures [18, 19]. In the case of quantum dots (QDs), an example of a simple approach based on the effective mass approximation is the “particle in a box model” [20] that can further be approximated by the celebrated (and analytically solvable) Fock-Darwin model [21, 22], which once offered general understanding of QD properties [23] and is still widely used, at least as the first approximation to many problems.

In a periodic system, the reduction of the multi-band $\bm{k}\!\cdot\!\bm{p}$ model to an effective mass theory is achieved via a quasi-degenerate perturbation theory that has its roots in quantum chemistry and is referred to as Löwdin partitioning [24]. By means of this procedure, one gets for the conduction band (cb) in the lowest (second) order a Schrödinger- (or Pauli-) like equation for the envelope function with a constant effective mass and Landé factor expressed via well-known formulas in terms of the band-edge energies. Two problems are encountered in an attempt to heuristically generalize such an equation to a modulated system (a nanostructure): First, it is generally accepted that (in analogy to the multi-band Hamiltonian), the components of the wave vector $\bm{k}$ should be replaced by derivatives (momentum operators in the coordinate representation). However, at the same time, the band edges and other parameters become position-dependent via their dependence on local composition and strain. This leads to a Schrödinger-like equation with a position-dependent effective mass and Landé factor and the ordering between the momenta and these two parameters starts to play a role. This ordering is arbitrary if the effective mass model is introduced as a heuristic extension of the bulk formulas. Arguments based on solvable models of abrupt interfaces [25, 26] lead to the conclusion that the only physically correct ordering in a non-homogenous medium is the most common “$k(1/m^{*})k$”. This conclusion has been formally verified by rigorous derivation from the exact multi-band envelope function equations [10, 2].

In this work we would like to address a more practical question of the achievable quantitative accuracy of the cb effective mass theories for quantum dot (QD) systems. Thus, the goal of this paper is to rigorously derive a family of effective mass Hamiltonians for a cb electron in a strained, inhomogeneous system by applying a systematic series of approximations to the Löwdin partitioning of the 8-band $\bm{k}\!\cdot\!\bm{p}$ Hamiltonian in the envelope function approximation and to validate their predictions against the results of the 8-band model with respect to the energy spectrum and Landé factors. By effective mass Hamiltonians we understand cb Hamiltonians for a carrier in an external magnetic field, obtained by a unitary transformation (partitioning) that eliminates the coupling between the conduction and valence bands of the original Hamiltonian, that contain quadratic (in momentum $\bm{k}$) kinetic energy terms and may contain cubic spin-orbit couplings. We will show that reproducing the results from the 8-band theory by an effective mass equation is possible with limited precision and the form of the equation required to achieve the best available accuracy is much more complicated than the usual Schrödinger- or Pauli-like equation for the envelope function. On the other hand, the Landé $g$ factor of the nanostructure can be reasonably estimated by the simple semi-phenomenological Roth-Lax-Zwerdling formula.

The starting point: 8-band $\bm{k}\!\cdot\!\bm{p}$ model

The 8-band $\bm{k}\!\cdot\!\bm{p}$ Hamiltonian in the envelope function approximation is defined in the block notation as[3]

[TABLE]

where the blocks refer in the standard way to the cb (6c), the $j=3/2$ valence band (vb, 8v) and the $j=1/2$ (spin-orbit split-off) vb (7v) and are explicitly given by[3, 14]

[TABLE]

Here $\{\mathcal{O}_{1},\mathcal{O}_{2}\}=\mathcal{O}_{1}\mathcal{O}_{2}+\mathcal{O}_{2}\mathcal{O}_{1}$, $k_{\{i}\mathcal{O}k_{j\}}=k_{i}\mathcal{O}k_{j}+k_{j}\mathcal{O}k_{i}$, $k_{[i}\mathcal{O}k_{j]}=k_{i}\mathcal{O}k_{j}-k_{j}\mathcal{O}k_{i}$ for any operators $\mathcal{O}$, $\mathcal{O}_{1}$, $\mathcal{O}_{2}$; $E_{\mathrm{c}}$ and $E_{\mathrm{v}}$ are the cb and vb edges, respectively ($E_{0}=E_{\mathrm{c}}-E_{\mathrm{v}}$ is the fundamental band gap in a bulk crystal); $\eta$ is the strain tensor corresponding to the static deformation due to the lattice mismatch; $\bm{k}=-i\nabla+e\bm{A}/\hbar$, where $\bm{A}$ is the vector potential of the magnetic field $\bm{B}$; $\tilde{\bm{k}}=\bm{k}(\mathbb{I}-\eta)$; $V_{\mathrm{p}}$ is the piezoelectric potential; $m_{0}$ is the free electron mass; $A^{\prime}_{c}$, $g^{\prime}$ and $\kappa^{\prime}$ are given by[3]

[TABLE]

$P=\hbar(E_{P}/2m_{0})^{1/2}$ (see below for the definition of $E_{P}$); $\gamma_{i}^{\prime}$ and $\kappa^{\prime}$ are the Luttinger parameters with removed contributions from the $\Gamma_{6}$ cb, $\gamma^{\prime}_{1}=\gamma_{1}-E_{P}/(3E_{0}+\Delta_{0})$, $\gamma^{\prime}_{2,3}=\gamma_{2,3}-E_{P}/(6E_{0}+2\Delta_{0})$, $\mu_{B}$ is the Bohr magneton; $q$ is another parameter of the Luttinger Hamiltonian; $B_{7\mathrm{v}}=(P^{\prime}Q/i)[1/(E_{0}-E_{0}^{\prime}-\Delta_{0}^{\prime})-1/(\Delta_{0}+E_{0}^{\prime}+\Delta_{0}^{\prime})]$, $B_{8\mathrm{v}}^{\pm}=(P^{\prime}Q/2i)[\pm 1/(E_{0}-E_{0}^{\prime}-\Delta_{0}^{\prime})\mp 1/(E_{0}^{\prime}+\Delta_{0}^{\prime})+1/(E_{0}-E_{0}^{\prime})-1/E_{0}^{\prime}]$, where $P^{\prime}$ and $Q$ are couplings to higher conduction bands; $\sigma_{i}$ are Pauli matrices; $J_{i}$ are matrices of the $j=3/2$ representation of angular momentum; $T_{i}$ are matrix representations of a vector operator between $j=1/2$ and $j=3/2$ states, i.e., $T_{x}=(T^{(1)}_{-1}-T^{(1)}_{+1})/\sqrt{2}$, $T_{y}=-(T^{(1)}_{-1}+T^{(1)}_{+1})/\sqrt{2}$, $T_{z}=T^{(1)}_{0}$, with the matrix elements of the spherical components $T^{(1)}_{q}$ given in terms of the Clebsch-Gordan coefficients $\langle j_{1}j_{2};m_{1}m_{2}|jm\rangle$ by the Wigner-Eckart theorem, $\langle m|T^{(1)}_{q}|m^{\prime}\rangle=-\sqrt{2/3}\langle 3/2,m^{\prime};1,q|1/2,m\rangle$, for $m=\pm 1/2$, $m^{\prime}=-3/2,\ldots,3/2$; and $T_{ij}=T_{i}J_{j}+T_{j}J_{i}$. The system is placed in an axial magnetic field. In numerical calculations we use gauge-invariant discretization scheme[13] for the covariant derivative.

The material parameters used in our $\bm{k}\!\cdot\!\bm{p}$ calculations are given in Table 1. In order to avoid $A^{\prime}_{c}<0$, which would break the ellipticity condition, we rescale $E_{P}$ to obtain $A^{\prime}_{c}=1$, which gives[27] $E_{P}=(m_{0}/m^{*}-1)E_{0}(E_{0}+\Delta_{0})/(E_{0}+2\Delta_{0}/3)$. Due to inconsistency of the reported values[3, 28], we calculate $q$ using the perturbative formula[29] $q=(2/9)E_{Q}[1/E^{\prime}_{0}-1/(E^{\prime}_{0}+\Delta^{\prime}_{0})]$, where $E_{Q}$, $E^{\prime}_{0}$ and $\Delta^{\prime}_{0}$ are $14$ band $\bm{k}\!\cdot\!\bm{p}$ parameters[3]. We account for the strain within a continuous elasticity approach[30]. Piezoelectric field in the system is calculated up to the second order in the polarization[31] with the parameters taken from Ref. [32].

Derivation of the effective mass Hamiltonian

The essence of the method [24, 34] is to perturbatively decouple the group of states of interest from all the other states by using a canonical transformation $T=e^{S}$, with an anti-hermitian operator $S$, in order to obtain a transformed Hamiltonian $\tilde{H}=THT^{{\dagger}}$, in which the inter-band terms (treated as a perturbation) are approximately eliminated. We will use a modified version of the van Vleck quasi-degenerate perturbation theory [35] as presented in Ref. [36]. The Hamiltonian is split into its block-diagonal and block-off-diagonal parts (coupling states within a single group of states and between the two groups of states, respectively), $H=H^{(\mathrm{d})}+H^{(\mathrm{od})}$. The operator $S$ is required to have null matrix elements within the groups (that is, it has to be block-off-diagonal). The group of states of interest here are the cb states. Since we start from the 8-band $\bm{k}\!\cdot\!\bm{p}$ theory, the other group of states are vb states.

We define the superoperator $\hat{S}$ representing the adjoint action of $S$ on the algebra of operators: $\hat{S}\mathcal{O}=[\mathcal{O},\hat{S}]$, $\hat{S}^{2}\mathcal{O}=[[\mathcal{O},\hat{S}],\hat{S}]$, etc. for any operator $\mathcal{O}$. Functions of $\hat{S}$ are defined via their power series expansion. Then, from the Campbell-Baker-Hausdorff expansion, $\tilde{H}=e^{S}He^{-S}=e^{\hat{S}}H=\cosh\hat{S}H+\sinh\hat{S}H$. Even functions of $\hat{S}$ transform block-diagonal operators into block-diagonal operators and block-off-diagonal operators into block-off-diagonal operators, and the opposite holds for odd functions of $\hat{S}$. Therefore, the block-diagonal and block-off-diagonal parts of $\tilde{H}$ are

[TABLE]

We require that the resulting Hamiltonian is block-diagonal, hence $\tilde{H}^{(\mathrm{od})}=0$. Inverting Eqs. (4) then yields

[TABLE]

In our case $H^{(\mathrm{od})}$ consists of $H_{\mathrm{6c8v}},H_{\mathrm{6c7v}}$ [Eq. (2d) and Eq. (2e)] and their hermitian conjugates, $H^{(\mathrm{od})}=H_{\mathrm{cv}}+\mathrm{h.c.}=H_{\mathrm{6c8v}}\oplus H_{\mathrm{6c7v}}+\mathrm{h.c.}$, hence it contains terms linear and quadratic in $\bm{k}$. Since $\tilde{H}^{(\mathrm{d})}$ contains $k$-independent terms, $S$ must be $O(k)$. According to the first of Eqs. (4), the leading order terms in $S$ yield quadratic corrections to $\tilde{H}^{(\mathrm{d})}$. We therefore take into account only the linear term in Eq. (5) and write

[TABLE]

The neglected corrections are $O(k^{3})$ and would lead to quartic terms in $\tilde{H}^{(\mathrm{d})}$, which is beyond the effective mass approximation, even with spin-orbit terms. Quantitatively, the linear truncation according to Eq. (6) amounts to neglecting corrections of relative magnitude $\Delta E/E_{0}$, where $\Delta E$ is the interband energy separation due to confinement. Taking $\Delta E$ as the excitation energy in the direction of strongest confinement ($\sim 200$ meV) one gets a rough estimate of 20% for the error due to truncation.

The remaining, block-diagonal part of the transformed Hamiltonian is now written as $\tilde{H}^{(\mathrm{d})}=H_{0}+\tilde{H}_{\mathrm{c}}^{\prime}\oplus\tilde{H}_{\mathrm{v}}^{\prime}$, where $H_{0}$ is assumed to be diagonal in a certain basis $|\alpha,i\rangle$ (where $\alpha$ denotes a group of states and $i$ are individual states within these groups), and is selected in such a way that the remaining parts are in some sense small. $\tilde{H}_{\mathrm{c}}^{\prime}$ and $\tilde{H}_{\mathrm{v}}^{\prime}$ denote the cb and vb blocks of $\tilde{H}$, respectively, with the corresponding parts of $H_{0}$ subtracted. In the problem at hand, we choose for $H_{0}$ a diagonal Hamiltonian which is constant and proportional to unity within each of the 6c, 8v and 7v bands and approximately represents the band edges in a strained QD. We will denote the respective energy values by $\overline{E}_{\alpha}$, $\alpha=\mathrm{6c,8v,7v}$. One can understand them as the average band edges corrected by hydrostatic strain (while the splitting between heavy and light holes within the 8v band is not included). The operator $S$ is written as $S=S_{\mathrm{cv}}+\mathrm{h.c.}$, where $S_{\mathrm{cv}}$ denotes one of the two off-diagonal blocks (a $2\times 6$ matrix in the standard $\bm{k}\!\cdot\!\bm{p}$ matrix notation). The effective band gaps are denoted by $\overline{E}_{\mathrm{6c}}-\overline{E}_{\mathrm{8v}}=E_{\mathrm{g}}$ and $\overline{E}_{\mathrm{6c}}-\overline{E}_{\mathrm{7v}}=E_{\mathrm{g}}+\Delta_{\mathrm{SO}}$. Finally we define a diagonal operator $\hat{\Delta}=\mathrm{diag}(E_{\mathrm{g}},E_{\mathrm{g}},E_{\mathrm{g}},E_{\mathrm{g}},E_{\mathrm{g}}+\Delta_{\mathrm{SO}},E_{\mathrm{g}}+\Delta_{\mathrm{SO}})$, where the entries correspond to the 6 valence bands of the 8-band $\bm{k}\!\cdot\!\bm{p}$ model. Then, from Eq. (6), one finds

[TABLE]

Note that the arbitrariness of choosing the diagonal Hamiltonian $H_{0}$ is removed here, as the subtracted energies (the operator $\hat{\Delta}$) are added back to the remaining part of the Hamiltonian.

Eq. (7) has a structure of a Sylvester equation but the operators appearing here are not finite-dimensional matrices. One can treat this equation as a matrix one, in the sense of the block notation over the subbands, but then the problem of non-commutativity of the matrix elements appears (due to non-commutativity of $\bm{k}$ with position-dependent quantities), precluding the application of standard algebraic methods for solving this equation. In order to overcome this difficulty, we expand the operators $S_{\mathrm{cv}}$ and $H_{\mathrm{cv}}$ in powers of $\bm{k}$,

[TABLE]

where the coefficients $C^{(1)}$ and $C^{(2)}$ are defined by comparison with the explicit form of Eq. (2d) and Eq. (2e). Then, upon rearrangement of terms one gets from Eq. (7) in the subsequent (formal) orders in $\bm{k}$

[TABLE]

The non-commutativity problem persists since $\tilde{H}_{\mathrm{c}}^{\prime}$ and $\tilde{H}_{\mathrm{v}}^{\prime}$ contain $k$-dependent terms, while $X^{(n)}$ are position-dependent. Returning to Eq. (6) one can see that the terms of $\tilde{H}^{(\mathrm{d})}$ linear and quadratic in $k$ generate corrections to $S$ on the order of $k^{2}$ and $k^{3}$, respectively. According to Eq. (4), these corrections generate terms $O(k^{3})$ and $O(k^{4})$ in the effective mass Hamiltonian for the cb. The latter are beyond the usual effective mass approximation, while the former correspond to spin-orbit terms but (by a simple perturbation argument) appear with a coefficient $\sim P^{2}C_{k}/E_{\mathrm{g}}^{2}\sim 0.6$ nm${}^{3}\cdot$meV (using InAs parameters), which is two orders of magnitude smaller than the Dresselhaus coefficient for InAs, $\alpha_{\mathrm{D}}=27$ nm${}^{3}\cdot$meV. It appears, therefore, that the kinetic part of $\tilde{H}_{\mathrm{v}}^{\prime}$ can be discarded in the derivation of a usual effective mass equation in the parabolic band approximation, that is, one with $O(k^{2})$ kinetic terms and the relevant $O(k^{3})$ spin-orbit corrections. However, as we will see below, including non-parabolicity effects at this point improves accuracy of the modeling of a self-assembled QD. Therefore, it seems reasonable to keep the kinetic part of $\tilde{H}_{\mathrm{c}}^{\prime}$ and $\tilde{H}_{\mathrm{v}}^{\prime}$. In order to obtain a solvable system of equations, we therefore propose to self-consistently replace the $\bm{k}$-dependent terms by their averages in the eigenstate of interest.

In our derivation, Eqs. (8a)–(8c) depend on the blocks of the transformed Hamiltonian $\tilde{H}$ rather than on the initial Hamiltonian $H$ hence, together with Eq. (4), they form a system that cannot be solved in a closed form. Clearly, in the leading order one could replace $\tilde{H}_{\mathrm{c}}^{\prime}$ and $\tilde{H}_{\mathrm{v}}^{\prime}$ by the original blocks $H_{\mathrm{c}}^{\prime}$ and $H_{\mathrm{v}}^{\prime}$ (the correction is $O(k^{2})$, yielding corrections that are formally $O(k^{4})$ in the resulting cb Hamiltonian). However, as we will see, including at least some corrections to these blocks improves the accuracy of the equation. It is known that the major corrections to the cb and vb Hamiltonians resulting from the decoupling procedure are the renormalization of the cb electron mass and of the Luttinger parameters, respectively. Therefore, we propose to take into account these strong effects only and to use, in place of $\tilde{H}_{\mathrm{c}}^{\prime}$ and $\tilde{H}_{\mathrm{v}}^{\prime}$, the cb and vb blocks of the original 8-band Hamiltonian but with the position-dependent renormalized parameters,

[TABLE]

where $E_{\mathrm{g}}^{\prime}=E_{0}+(a_{\mathrm{c}}-a_{\mathrm{v}})\operatorname{Tr}\eta$ is the local band gap, including the hydrostatic strain-induced shift. Note that these band-decoupling corrections are $O(k^{2})$, so the distinction between $\tilde{H}_{\mathrm{c,v}}^{\prime}$ and $H_{\mathrm{c,v}}^{\prime}$ is only important if the $k$-dependent terms are included self-consistently, as proposed above.

In the approximations proposed here, Eqs. (8a)–(8c) are a system of usual Sylvester equations that can be solved iteratively. A solution can be obtained in a closed, analytical form using the general method of Ref. [37]. However, the form of the solution simplifies considerably if one discards the contribution of the spin-dependent part of $\tilde{H}_{\mathrm{c}}^{\prime}$ to the operator $S$ (these terms remain included to the leading order in $H_{\mathrm{6c6c}}$ that is part of $H^{(\mathrm{d})}$, see Eq. (2a) and Eq. (4)). These terms are very small compared to any other energy scales in the problem, hence their contribution to $S$ is negligible. Within this approximation, one has $\tilde{H}_{\mathrm{c}}^{\prime}=\chi_{\mathrm{c}}^{\prime}\mathbb{I}_{2\times 2}$, where $\chi_{\mathrm{c}}^{\prime}$ is a scalar function of position and $\mathbb{I}_{2\times 2}$ is a $2\times 2$ unit matrix. Then, the solution to the above system of equations can be obtained trivially. Denoting $\mathcal{D}=\hat{\Delta}+\chi_{\mathrm{c}}^{\prime}\mathbb{I}_{6\times 6}-\tilde{H}_{\mathrm{v}}^{\prime}$ one has

[TABLE]

From Eqs. (4), with the condition $\tilde{H}^{\mathrm{od}}=0$, one finds [36] $\tilde{H}^{\mathrm{d}}=H^{\mathrm{d}}+\tanh(\hat{S}/2)H^{\mathrm{od}}$ or, in the leading order, $\tilde{H}^{\mathrm{d}}\approx H^{\mathrm{d}}+(1/2)[H^{\mathrm{od}},S]$. Then, the correction to the cb Hamiltonian up to the order $k^{3}$ can be decomposed into two parts. The first one, which we will denote by $\tilde{H}^{(2)}$, is formally quadratic in $\bm{k}$ and is proportional to $C_{j}^{(1)}C_{l}^{(1)}$ in our notation. This contribution yields corrections to the electron effective mass and Landé factor. The second part, denoted $\tilde{H}^{(3)}$, is of third order in $\bm{k}$ and contains terms proportional to $C_{jl}^{(2)}C_{n}^{(1)}$. It includes the Dresselhaus spin-orbit term. Thus, the resulting transformed Hamiltonian can be written as $\tilde{H}=H_{\mathrm{c}}+\tilde{H}^{(2)}+\tilde{H}^{(3)}$. From Eq. (4) in the linear approximation one finds

[TABLE]

and

[TABLE]

The first term in Eq. (11) leads to the usual Dresselhaus spin-orbit coupling. The other terms are linear or quadratic in $k$ and yield a small correction to the kinetic and Zeeman terms of the effective mass Hamiltonian.

Interpretation, verification and discussion

In this Section, we first present various approximations to the effective mass Hamiltonian and relate the equations resulting from some of these approximations to the common form of the effective mass Hamiltonian, written in terms of the effective mass tensor and the $g$-factor given by the Roth formula. Next we define a series of approximation for which we present quantitative comparison of the predictions from the effective mass theory with the results from the 8-band $\bm{k}\!\cdot\!\bm{p}$ Hamiltonian.

Interpretation of the effective mass equation

We restrict the discussion to the quadratic term [Eq. (10)] that determines the fundamental properties of the energy spectrum.

In order to simplify notation, we define the $2\times 6$ matrices $\mathcal{T}_{i}=\sqrt{3}T_{i}\oplus(-1/\sqrt{3})\sigma_{i}$, where the two components of the direct sum correspond to the $j=3/2$ (hh and lh) and $j=1/2$ (spin-orbit split-off) subbands of the vb. Then, by direct inspection of Eq. (2d) and Eq. (2e) one finds $C_{j}^{(1)}=-P\widetilde{\mathcal{T}}_{j}$, where $\widetilde{\mathcal{T}}_{j}=\sum_{n}(\delta_{jn}-\eta_{jn})\mathcal{T}_{n}$. The second-order correction to the effective mass Hamiltonian is then

[TABLE]

where we use the fact that $[k_{j},\chi_{\mathrm{c}}^{\prime}]$ is a number and commutes with the matrices $\widetilde{\mathcal{T}}_{i}$ and $\mathcal{D}$.

It might seem that the above Hamiltonian confirms the correctness of the particular (BenDaniel-Duke[38]) ordering of the operators in the kinetic term. This is not the case: a simple manipulation of the terms in Eq. (12) under assumption $P=\mathrm{const}$ allows one to rewrite the second-order correction in the equivalent form

[TABLE]

where a different (Gora-Williams-Bastard[39, 40, 41]) ordering appears in the kinetic term. Interestingly, the correcting terms in Eq. (12) and Eq. (13) involve only spatial derivatives of cb and vb parameters, respectively. It is hence clear that the two orderings are, in a sense, dual and each of them corresponds to neglecting terms in the Hamiltonian, related to spatial modulation of either the conduction or valence band. As we will see explicitly below, these additional terms are essential to correctly reproduce the Rashba spin-orbit interaction.

Further analytical insight into the somewhat unusual form of our effective mass Hamiltonian $\tilde{H}$ is hindered by the need to invert the “band gap operator” $\mathcal{D}$, which produces rather intransparent and intractable formulas (and, in practice, is performed numerically). One might expect that neglecting the off-diagonal elements of $\mathcal{D}$, from which this difficulty stems, is a good approximation, since these elements are rather small compared to the band gap. The operator $\mathcal{D}$ in this approximation will be denoted by $\tilde{\mathcal{D}}$. Since it is diagonal its inverse powers are found trivially. Let us denote its diagonal elements by $(E_{\mathrm{hh}},E_{\mathrm{lh}},E_{\mathrm{lh}},E_{\mathrm{hh}},E_{\mathrm{so}},E_{\mathrm{so}})$. They can be interpreted as the offset of the local edges of the three valence subbands with respect to the cb edge at a given point in space. We split the $\mathcal{T}_{j}\tilde{\mathcal{D}}^{-n}\mathcal{T}_{l}^{{\dagger}}$ matrix into a symmetric and asymmetric part,

[TABLE]

Using the explicit forms of the matrices $\mathcal{T}_{j}$ one finds up to the linear order in $\eta$

[TABLE]

where $F_{x}^{(n)}=F_{y}^{(n)}=E_{\mathrm{P}}(2/E^{n}_{\mathrm{so}}+1/E^{n}_{\mathrm{lh}}+3E^{n}_{\mathrm{hh}})/6$, $F_{z}^{(n)}=E_{\mathrm{P}}(1/E^{n}_{\mathrm{so}}+2/E^{n}_{\mathrm{lh}})/3$. Thus, the symmetric part is spin-independent and contributes to the kinetic part of the effective mass Hamiltonian in the well-known way via the electron effective mass tensor; indeed, the first term in Eq. (12) combined with the kinetic term in $H_{\mathrm{6c6c}}$ [Eq. (2a)] is

[TABLE]

The second term yields

[TABLE]

Thus, the symmetric component of the first term in Eq. (12) yields the usual kinetic energy of the effective mass theory (Eq. (15)) in the commonly used ordering “$k(m^{*})^{-1}k$” (BenDaniel-Duke[38] ordering). This is, however, corrected by the second term (Eq. (16)) and one could as well start from Eq. (13) and arrive at an equivalent Hamiltonian with Gora-Williams-Bastard[39, 40, 41] ordering and with an $H_{\mathrm{v}}$-dependent correction term instead of Eq. (16).

The asymmetric part of $H^{(2)}$ can be interpreted most easily starting form the alternative form of the Hamiltonian given in Eq. (13). The asymmetric contribution to the first term in Eq. (13), combined with the Zeeman term of $H_{\mathrm{6c6c}}$ [Eq. (2a)] can be written in the form

[TABLE]

In the diagonal approximation, $\mathcal{D}\approx\tilde{\mathcal{D}}$, the asymmetric part yields up to order $O(\eta)$

[TABLE]

where

[TABLE]

Using the relation $\left[k_{j},k_{l}\right]=-\frac{ie}{\hbar}\varepsilon_{jlm}B_{m^{\prime}}$, and defining the (diagonal) Landé tensor $\hat{g}_{mm^{\prime}}=g^{\prime}\delta_{mm^{\prime}}-\left[\delta_{mm^{\prime}}(1+\operatorname{Tr}\eta)-\eta_{mm^{\prime}}\right]\delta g_{m}^{(1)}$, one gets the usual Zeeman Hamiltonian

[TABLE]

The form of the Landé tensor $\hat{g}$ appearing in this equation is close to the widely used isotropic Roth-Lax-Zwerdling[42] formula

[TABLE]

Here, $\hat{g}$ differs from $g^{{}^{(\mathrm{Roth})}}$ by including strain effects as well as anisotropy resulting from the hh-lh splitting in the nanostructure. The Roth-Lax-Zwerdling formula has the advantage that the fundamental band gap $E_{\mathrm{g}}$, instead of being calculated at some level of approximation to the $\bm{k}\!\cdot\!\bm{p}$ theory, can be taken from experiment[43, 44, 45]. While the standard Roth-Lax-Zwerdling formula is linked to the 8-band $\bm{k}\!\cdot\!\bm{p}$ Hamiltonian via a diagonal and parabolic approximations, the full effective mass Hamiltonian accounts also for non-parabolicity corrections and the full structure of the vb.

The asymmetric part of the second term in Eq. (13) is a generalized Rashba term. To see this, assume that the spatial variation of $H_{\mathrm{v}}^{\prime}$ results from an external electric field $\bm{\mathcal{E}}$, hence $[k_{j},H_{\mathrm{v}}^{\prime}]=-i\nabla H_{\mathrm{v}}^{\prime}=-ie\bm{\mathcal{E}}\mathbb{I}_{6\times 6}$. Then, taking only the strain-independent terms in the second term in Eq. (13) and using Eq. (17) one finds in the diagonal approximation

[TABLE]

Neglecting the anisotropy due to the lh-hh splitting in Eq. (18) one finds

[TABLE]

which is the standard Rashba spin-orbit term [3]. The full Eq. (13) yields the Rashba term generalized to arbitrary inhomogeneity of the valence band and includes corrections due to inhomogeneous strain. The same effects are accounted for (although not so explicitly) by Eq. (12).

Quantitative assessment

In the following, we compare the accuracy (with respect to the 8-band $\bm{k}\!\cdot\!\bm{p}$ results) of various approximations to the effective mass equation derived above. Our systematic approach allows us to discuss the possible approximate equations in the uniform framework, according to the level of approximation made to the $\mathcal{D}$ operator.

(1) As the very first step, obviously inaccurate but included for completeness, we apply the bulk approximation without strain, in which the operator $\mathcal{D}$ is taken diagonal and without the kinetic ($k^{2}$) terms,

[TABLE]

with local bad gaps calculated from unstrained bulk values, $E_{\mathrm{hh}}=E_{\mathrm{lh}}=E_{0}$, $E_{\mathrm{so}}=E_{0}+\Delta_{0}$, where the band edges are interpolated according to the local composition (see Tab. 1).

(2) The first reasonable approximation is to correct the above procedure by accounting for band-edge shifts due to local strain, which is done by including the diagonal strain-dependent matrix elements in $H_{\mathrm{v}}$ and $H_{\mathrm{c}}$. Thus, the $\mathcal{D}$ operator is still given by Eq. (21) but with $E_{\mathrm{hh}}=E_{0}+(a_{\mathrm{c}}-a_{\mathrm{v}})\operatorname{Tr}\eta+b_{\mathrm{v}}(2\eta_{zz}-\eta_{xx}-\eta_{yy})/2$, $E_{\mathrm{lh}}=E_{0}+(a_{\mathrm{c}}-a_{\mathrm{v}})\operatorname{Tr}\eta-b_{\mathrm{v}}(2\eta_{zz}-\eta_{xx}-\eta_{yy})/2$, $E_{\mathrm{so}}=E_{0}+\Delta_{0}+(a_{\mathrm{c}}-a_{\mathrm{v}})\operatorname{Tr}\eta$, where the local value of strain is used and parameters are interpolated according to the local composition. In this way, the major correction to the band gap (which is a crucial factor determining the effective mass and $g$ factor) is included in the model. This will be referred to as bulk approximation with strain.

(3) The value used for the band gap in the previous approach is still a crude approximation since the actual energy spacing between the states in the conduction and valence bands in a QD are affected by spatial confinement. A common approach, taken in many cases[43, 44, 45] to estimate the electron $g$-factor when interpreting experimental data, is to use the measured value of the fundamental transition energy $E_{\mathrm{g}}^{(\mathrm{exp})}$ as the effective band gap value for a given system. In our numerical study this corresponds to replacing the fundamental band edge offset in the operator $\mathcal{D}$ (together with its strain correction) by the splitting between the top vb and bottom cb states obtained from the full 8-band calculation, while using the bulk values of the spin-orbit splitting $\Delta_{0}$ and the Kane parameter $P$ interpolated according to composition. The $\mathcal{D}$ operator is still given by the diagonal form of Eq. (21) but now $E_{\mathrm{hh}}=E_{\mathrm{lh}}=E_{\mathrm{g}}^{(\mathrm{exp})}$, $E_{\mathrm{so}}=E_{\mathrm{g}}^{(\mathrm{exp})}+\Delta_{0}$, with a constant $E_{\mathrm{g}}^{(\mathrm{exp})}$ and a local composition-dependent $\Delta_{0}$. We will refer to this as semi-phenomenological approximation. For the $g$-factor, a further simplification along this line is achieved by introducing constant values of $\Delta_{0}$ obtained by averaging the position-dependent values weighted by the squared wave function, which yields an explicit effective Roth-Lax-Zwerdling formula for the $g$-factor of a nanostructure, that is, Eq. (20) with the effective values of the parameters.

(4) In the off-diagonal approximation, we include the full structure of the vb Hamiltonian $H_{\mathrm{v}}$, hence the full form of $\mathcal{D}$, but still neglect the $k$-dependent terms. In block notation,

[TABLE]

where $H_{\mathrm{8v8v}}^{(x)}$, $H_{\mathrm{8v7v}}^{(x)}$, and $H_{\mathrm{7v7v}}^{(x)}$ are the position-dependent parts of the corresponding vb Hamiltonian blocks, formally obtained by setting $k=0$ in Eq. (2b), Eq. (2f), and Eq. (2c), respectively, and $H_{\mathrm{7v8v}}^{(x)}=H_{\mathrm{8v7v}}^{(x)\dagger}$. This yields a cb Hamiltonian strictly quadratic in $k$ (and strictly equivalent to the original 8-band Hamiltonian up to the quadratic order), which corresponds to the usual notion of the effective mass equation in the parabolic approximation. Starting from this approximation, we go beyond the usual forms of the effective mass equations by introducing the full matrix structure of the vb-induced corrections.

(5) The next step consists in including the self-consistent averages of the $k$-dependent terms in $\mathcal{D}$, that is

[TABLE]

where $H_{\mathrm{8v8v}}^{(\mathrm{av})}$, $H_{\mathrm{8v7v}}^{(\mathrm{av})}$, and $H_{\mathrm{7v7v}}^{(\mathrm{av})}$ are obtained from Eq. (2b), Eq. (2f), and Eq. (2c), respectively, by self-consistently replacing each product $k_{i}k_{j}$ with its average (denoted by $\langle\ldots\rangle$) in the state of interest. In this way we go beyond the parabolic approximation. We will call this the off-diagonal$+\langle k^{2}\rangle$ approximation. In this approximation the implicit character of Eqs. (8a)–(8c) (the appearance of $\tilde{H}_{\mathrm{c}}$ and $\tilde{H}_{\mathrm{v}}$) is resolved to the leading order, that is, by replacing the two blocks of the Hamiltonian by the original $H_{\mathrm{c}}$ and $H_{\mathrm{v}}$.

(6) The final approach in our sequence of approximations consists in approximating $\tilde{H}_{\mathrm{c}}$ and $\tilde{H}_{\mathrm{v}}$ in Eqs.  (8a)–(8c) self-consistently by renormalizing the effective mass and Luttinger parameters according to Eqs. (9). The form of $\mathcal{D}$ is therefore the same as in Eq. (22) but with the renormalized parameters. We will refer to this as the self-consistent effective mass equation.

As a test of these approximations, we study the low-energy part of the spectrum of a self-assembled QD. We consider two models of a lens-shaped, self-assembled InAs/GaAs QD. In both models the QD has $24$ nm diameter and $4.2$ nm height but they differ in the composition profile: While the first model assumes a uniform composition with 100% InAs inside the QD, the second one has a more realistic trumpet-shape profile of InAs/GaAs composition with the InAs content defined by [46] $C(\bm{r})=C_{\mathrm{b}}+(C_{\mathrm{t}}-C_{\mathrm{b}})\exp[-\sqrt{x^{2}+y^{2}}\exp(-z/z_{\mathrm{p}})/r_{\mathrm{p}}]$, where we took $C_{\mathrm{b}}=0.4$ , $C_{\mathrm{t}}=0.8$, $r_{\mathrm{p}}=0.9$ nm and $z_{\mathrm{p}}=1.4$ nm. In both cases the QD is placed on a $0.6$ nm wetting layer (WL), which in the case of the first dot contains 100% and in the latter one 40% InAs. The two figures of merit that we will investigate here are the energy splitting $\Delta E_{sp}$ between the ground (“$s$-shell”) electron state and the lowest excited (“$p$-shell”) state at zero magnetic field (Fig. 1) and the ground state $g$-factor, extracted from the leading order (linear) term of the Zeeman splitting of the ground-state doublet at low magnetic fields (Fig. 2). In both cases we show in the plots the values obtained using only the first, usual term in Eq. (12) (points) and those from the full equation, including also the non-standard second term (crosses).

As shown by the results in Fig. 1, the effective mass methods typically fail to reproduce the $s$-$p$ shell splitting which, in fact, belongs to the most fundamental quantitative characteristics of a QD system. What matters here is the renormalization of the $\bm{k}\!\cdot\!\bm{p}$ parameters according to Eqs. (9), hence the self-consistent equation (method (6) from our series of approximations) is able to produce a result that is closest to the correct value. For the uniform QD, Fig. 1(a), the disagreement is still at the level of 30%, while for the more realistic model with the trumpet-shape composition profile Fig. 1(b), the effective mass result is much more exact. The 30% discrepancy is consistent with the estimate of the error induced by the linear truncation in Eq. (6). The much higher precicision in the case of the smooth trumpet-shape composition may be due to the fact that terms of higher order in $k$ in an envelope-function theory for an inhomoheneous system are not just “non-parabolic” but also induce higher order derivatives of system parameters, to which the uniform QD model with abrupt compositional boundaries is obviously much more sensitive. One can understand why approximation (6) is particularly suited for relatively correct modeling of this particular spectral feature by noting that the $s$-$p$ shell splitting is related to the in-plane excitation of the QD, for which the value of the effective mass plays a crucial role. The renormalization described by Eq. (9) is an important correction to this parameter. On the contrary, including the non-parabolicity correction, either by self-consistently adding the average $k^{2}$ terms in the $\mathcal{D}$ operator, as in the approach (5), or by using a phenomenological value of $E_{\mathrm{g}}$, as in the approach (3), mostly reflects the impact of the strong confinement in the growth direction, which shifts the lowest shells rigidly, hence does not affect the $s$-$p$ splitting considerably. It is found also that including the full matrix structure of $H_{\mathrm{v}}$ in $\mathcal{D}$ (approximation (4)) does not bring any clear advantage by itself. In addition, it turns out that correcting the effective mass equation by the second term in Eq. (12) has an opposite effect on the two models: it decreases the accuracy for the homogeneous QD, while it improves it (to a larger degree) for the trumpet-shape model. We were not able to find a plausible explanation of this fact. Interestingly, the correction stemming from this term is roughly the same for a given QD structure in each approximation.

The accuracy of various approximations to the effective mass equation is different when the $g$-factor is considered (Fig. 2). Here the approximations (4-6), where the full matrix structure of the vb is included, yield much more accurate results than the diagonal approximations (1) and (2). As one can see in (Fig. 2), the self-consistent model (6), with self-consistently renormalized parameters, again provides the most accurate results in the systematic series of approximations. The discrepancy is about 60% for the uniform QD model (Fig. 2a) while the result for the model with a trumpet-shape composition profile (Fig. 2b) is overestimated by a factor of 2. In the case of the $g$-factor, including the part of non-parabolicity that is captured by our self-consistent approach improves the accuracy to some extent but is not critical. Also the correction provided by the second term in Eq. (12) is of much less importance than it was for the spectrum, with a larger relative contribution in the trumpet-shape model. In both models, it turns out to lower the accuracy. The large error in the values of the $g$-factor apparently exceeds the estimate of the effect of truncating Eq. (6). In fact, however, the numerical values should be referred to the bulk value of $-15$, which is nearly entirely compensated by confinement and strain effects. Therefore, the mismatch is actually of only a few percent.

It seems interesting that accurate systematic modeling of the Zeeman splitting in an effective mass equation requires the full form of the operator $\mathcal{D}$. This means that the standard Eq. (19) is not very useful if the $g$-factor is calculated by Löwdin perturbation theory including only the diagonal contributions to the vb structure (in particular, the hh-lh splitting induced by the axial-strain-related $b_{\mathrm{v}}$ terms): for our two structures the approach (2) overestimates the Zeeman splitting by an order of magnitude. We were able to trace this back to the important role played by the off-diagonal terms proportional to the $b_{\mathrm{v}}$ deformation potential in the correct alignment of valence bands: including axial strain contributions only in the diagonal terms produces wrong vb alignment (in terms of band splitting and the sign of the band offsets between the QD and the barrier). As a consequence, completely neglecting these subtle details of the vb, as done in the semi-phenomenological approximation (3), leads to a result which is only slightly worse than the most exact one when compared to the full 8-band model. Here the lh-hh splitting and band offsets are completely washed out by setting a single constant value for the fundamental band gap and only the spin-orbit ($\Gamma_{7}$) band splitting is position dependent. The result becomes even more accurate if the latter is also replaced by a single number obtained from spatial averaging of the standard bulk values (black triangles in (Fig. 2)). Thus, our numerical results support the commonly used and very convenient Roth-Lax-Zwerdling formula for estimating the electron $g$-factor.

Conclusions

In this paper we have proposed a systematic derivation of a cb effective mass equation from the 8-band $\bm{k}\!\cdot\!\bm{p}$ Hamiltonian. We have shown that this derivation develops into a systematic series of approximations that differ in the way the vb is represented in the final equation. Possible approximations range from using a set of fixed band gap parameters to a full matrix structure with self-consistent non-parabolicity corrections and parameter renormalization. We have assessed the accuracy of the approximations in calculating selected spectral and spin-related characteristics for a self-assembled QD system within two models of the composition profile.

We have shown that a quantitatively correct description of the lowest sector of the electron spectrum, which involves intraband dynamics and therefore relies primarily on the accurate modeling of the effective mass, requires a self-consistent renormalization of the Hamiltonian parameters that goes beyond the second order Löwdin perturbation. The accuracy is also improved by accounting for cb non-parabolicity by self-consistently including terms of higher order in the electron momentum. When studying the Zeeman splitting of the ground electron level we found that the most accurate value of the $g$ factor is obtained within the systematic scheme only after including the full structure of the vb Hamiltonian in the equation. Here, again, including non-parabolicity corrections improves the accuracy of the result. Surprisingly, the values obtained from the resulting rather complicated equation can be reproduced by a simple diagonal model with a fixed value of the fundamental band gap, which supports the use of the effective Roth-Lax-Zwerdling formula for a nanostructure. The accuracy of the latter is remarkable, given its simplicity. The effective mass equations reproduce 8-band $\bm{k}\!\cdot\!\bm{p}$ value of the electron $g$-factor for the QD ground state within a factor of 2, which may seem disappointing. One should take into account, however, that the value of the $g$-factor for bulk InAs is about $-15$, which is nearly entirely compensated by confinement and strain effects. The absolute mismatch, which is on the order of 0.1, is a tiny fraction of the original value (that is, the compensation is quantitatively reproduced) but yields a large relative error when compared to the very small final value.

In general, we have shown that the effective mass equation offers a limited accuracy when modeling the lowest-energy sector of the electron spectrum in a self-assembled QD, unless it is extended to a rather complicated, inconvenient and computationally expensive form. In particular, an equation that is rigorously derived from the $\bm{k}\!\cdot\!\bm{p}$ theory strictly up to order $k^{2}$ (which may correspond to the most usual notion of an effective mass theory) quantitatively fails in all respects. The more accurate self-consistent equation proposed in this paper is not particularly transparent and does not even allow a separation into kinetic and Pauli terms. On the other hand, we have shown that the electron $g$-factor can be estimated by a very simple version of the effective mass theory, thus justifying the phenomenological Roth-Lax-Zwerdling formula for a nanostructure. Our derivation highlights the correspondence between the BenDaniel-Duke and Gora-Williams-Bastard orderings of the kinetic term, both emerging as different approximations within a general scheme, and yields a generalized Rashba coupling including inhomogeneity and strain effects. This confirms once more that the effective mass equation, even though not perfectly accurate, can be very useful from the conceptual point of view and often provides valuable physical insight.

Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Acknowledgments

This work was supported by the Polish National Science Centre (Grant No. 2014/13/B/ST3/04603). Calculations have been partly carried out in Wroclaw Centre for Networking and Supercomputing (http://www.wcss.wroc.pl), Grant No. 203. A. M.-P. acknowledges support from the subsidy granted by the Polish MNiSW to the Faculty of Fundamental Problems of Technology for the research activity contributing to the development of young researchers and PhD students.

Author contributions statement

All authors performed analytical calculations. K.G. developed most of the code, with contributions from A.M.-P. Simulations were performed by A.M.-P. and K.G. All authors contributed equally to analyzing and interpreting the results. All authors contributed to writing and reviewing the manuscript.

Additional information

The authors declare that they have no competing interests.