The importance of second order deformation potentials in modeling of InAs/GaAs nanostructures

TL;DR

This paper introduces a non-linear 8-band k.p modeling approach with strain terms fitted to tight-binding results, improving accuracy in simulating InAs/GaAs nanostructures while maintaining computational efficiency.

Contribution

It develops and validates a non-linear 8-band k.p method with strain terms fitted to atomistic results for better modeling of InAs/GaAs nanostructures.

Findings

Good agreement with tight-binding results for electron and hole states.

Provides non-linear k.p parameters for InAs and GaAs.

Proposes an improved effective mass parameterization.

Abstract

Accurate modeling of electronic properties of nanostructures is a challenging theoretical problem. Methods making use of continuous media approximation, such as k.p, sometimes struggle to reproduce results obtained with more accurate atomistic approaches. On the contrary, atomistic schemes generally come with a substantially larger cost of computation. Here, we bridge between these two approaches by taking 8-band k.p method augmented with non-linear strain terms fit to reproduce sp3d5s* tight-binding results. We illustrate this method on the example of electron and hole states confined in quantum wells and quantum dots of photonics applications relevant InAs/GaAs material system, and demonstrate a good agreement of a non-linear k.p scheme with empirical tight-binding method. We discuss limits of our procedure as well as provide non-linear 8-band k.p parameter sets for InAs and GaAs. Finally, we propose a parameterization for effective term used to improve the accuracy of the standard effective mass method.

Tables

Table 1. Table 1: Deformation potentials (in eV), and parameters of the valence force field model, used in the calculations.

	GaAs	InAs
$a_c^{\left(1\right)}$	-6.79f(-7.17a)	-4.78f(-5.08a)
$a_v^{\left(1\right)}$	1.84f(1.16a)	1.24f(1.00a)
$b_v^{\left(1\right)}$	-1.85f(-2.0a)	-1.77f(-1.8a)
$d_v^{\left(1\right)}$	-4.8a	-3.6a
$C_2$	2.5 – 5.5 (3.3)	6.0f
$a_c^{\left(2a\right)}$	-2.71f	-3.40f
$a_c^{\left(2b\right)}$	24.6f	18.1f
$a_v^{\left(2a\right)}$	-3.56f	-1.56f
$a_v^{\left(2b\right)}$	4.81f	1.07f
$b_v^{\left(2a\right)}$	-6.38f	-5.95f
$b_v^{\left(2b\right)}$	-9.23f	-6.64f
${\tilde{a}}_c^{\left(2\right)}$	-3.8f	-3.2f
$\alpha$ $({10}^3$ dyne $)$	$41.49$b	$35.18$b
$\beta$ $({10}^3$ dyne $)$	$8.94$b	$5.49$b
fValues fitted to the tight-binding band structure.
aValues taken from Ref.vurgaftman01 .
bValues taken from Ref.pryor98b .

Table 2. Table 2: Band structure parameters used in the calculations. We show the values of reduced Luttinger parameters (where the contribution from the Γ 6 ​ c subscript Γ 6 c \Gamma_{\mathrm{6c}} is substracted), which directly enter the Hamiltonian Winkler2003 .

	GaAs	InAs
$E_v$	$0.0$ eV	$0.21$ eVa
$E_g$	1.519 eVa	0.418 eVf
$E_p$	$21.5$ eVf	$19.5$ eVf
$P$	calculated from $P=\sqrt{E_p{\hslash }^2/\left(2m_0\right)}$
${\gamma }_1^{\prime }$	2.26a	2.15f
${\gamma }_2^{\prime }$	-0.299a	-0.325f
${\gamma }_3^{\prime }$	0.571a	0.542f
$C_k$	-0.0034 $\mathrm{eV}\AA$b	-0.0112 $\mathrm{eV}\AA$b
${\Delta }_0$	0.341 eVa	0.38 eVf
fValues fitted to the tight-binding band structure.
aValues taken from Ref.vurgaftman01 .
bValues taken from Ref.Winkler2003 .

Equations

$$t=t_0{\left(\frac{d_0}{d}\right)}^{\eta },$$

$$E_{xy}=E_d+2b_d\text{quantity}[{ϵ}_{zz}-({ϵ}_{xx}+{ϵ}_{yy})/2],$$

$$H_{k,6c6c}$$

$$H_{k,8v8v}$$

$$\displaystyle\phantom{=}-4\gamma^{\prime}_{3}\left[\{J_{x},J_{y}\}\{k_{x},k_{y}\}+{\mathrm{c.p.}}\right]\bigg{\}}$$

$$-i\frac{{\hslash }^2}{m_0}\text{quantity}({\kappa }^{\prime }[k_x,k_y]J_z+q^{\prime }[k_x,k_y]J_z^3+c.p.)$$

$$+\frac{2}{\sqrt{3}}{C}_k\left[\{J_x,J_y^2-J_z^2\}{k}_x+c.p.\right],$$

$$H_{k,7v7v}$$

$$-i\frac{{\hslash }^2}{m_0}\text{quantity}({\kappa }^{\prime }[k_x,k_y]{\sigma }_z+c.p.),$$

$$H_{k,6c8v}$$

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

$$-\frac{\sqrt{3}}{2}{B}_{8v}^{-}{T}_{zz}(k_x^2-k_y^2),$$

$$H_{k,6c7v}$$

$$H_{k,8v7v}$$

$$\displaystyle\phantom{=}-12\gamma^{\prime}_{3}(U_{xy}\{k_{x},k_{y}\}+{\mathrm{c.p.}})\big{]}$$

$$-i\frac{3{\hslash }^2}{2m_0}\text{quantity}({\kappa }^{\prime }{U}_z[k_x,k_y]+c.p.)$$

$$-i\sqrt{3}{C}_k\left(U_{yz}{k}_x+c.p.\right),$$

$$H_{\mathrm{str},6c6c}^{(1)}$$

$$H_{\mathrm{str},6c8v}^{(1)}$$

$$H_{\mathrm{str},6c7v}^{(1)}$$

$$H_{\mathrm{str},8v8v}^{(1)}$$

$$-\frac{d_{\mathrm{v}}^{(1)}}{\sqrt{3}}\left[2\{J_x,J_y\}{ϵ}_{xy}+c.p.\right],$$

$$H_{\mathrm{str},8v7v}^{(1)}$$

$$-\sqrt{3}{d}_{\mathrm{v}}^{(1)}\left(2U_{xy}{ϵ}_{xy}+c.p.\right),$$

$$H_{\mathrm{str},7v7v}^{(1)}$$

$$H_{\mathrm{str},6c6c}^{(2)}$$

$$H_{\mathrm{str},8v8v}^{(2)}$$

$$-b_{\mathrm{v}}^{(2a)}\left[\left(J_x^2-\frac{1}{3}{J}^2\right){ϵ}_{xx}^2+c.p.\right]$$

$$-b_{\mathrm{v}}^{(2b)}\left[\left(J_x^2-\frac{1}{3}{J}^2\right){ϵ}_{yy}{ϵ}_{zz}+c.p.\right],$$

$$H_{\mathrm{str},8v7v}^{(2)}$$

$$-3b_{\mathrm{v}}^{(2b)}\left(U_{xx}{ϵ}_{yy}{ϵ}_{zz}+c.p.\right),$$

$$H_{\mathrm{str},7v7v}^{(2)}$$

$${\mathbf{R}}_i=(\text{mathds}1+\widehat{ϵ}){{\mathbf{R}}^0}_i\pm \zeta \frac{a}{4}({ϵ}_{yz},{ϵ}_{zx},{ϵ}_{xy}),$$

$$\zeta =\frac{\alpha -\beta }{\alpha +\beta },$$

$${\tilde{H}}_{\mathrm{eff}}$$

$$+{\tilde{a}}_{\mathrm{c}}^{(2)}[({ϵ}_{xx}-{ϵ}_{yy}{)}^2+({ϵ}_{yy}-{ϵ}_{zz}{)}^2+({ϵ}_{zz}-{ϵ}_{xx}{)}^2],$$

$$U$$

$$+\frac{3}{8}\sum _i\sum _j^{\mathrm{NN}(i)}\sum _{k>j+1}^{\mathrm{NN}(i)}{B}_{ijk}{\left({\mathbf{r}}_{ij}{\mathbf{r}}_{ik}-d_{ij}{d}_{ik}\cos {\theta }_0\right)}^2,$$

$$H_{\mathbf{k}}$$

$$+\frac{\hslash }{4m_0^2{c}^2}(\mathbf{p}+\hslash \mathbf{k})\cdot \mathbf{\sigma }\times (\text{gradient}{V}_0(\mathbf{r})),$$

$$⟨\lambda \nu |H_0^{\prime }|\mu {\lambda }^{\prime }⟩$$

$$+\sum _{ijkl}{F}_{ijkl}^{\nu \mu }{ϵ}_{ij}{ϵ}_{kl}{\delta }_{\lambda {\lambda }^{\prime }},$$

$${\mathcal{D}}_S={\mathcal{D}}_A\times {\mathcal{D}}_B^{\ast }.$$

$$H_{\mathrm{str},6c6c}^{(2)}$$

$$+a_{\mathrm{c}}^{(2c)}({ϵ}_{xy}^2+c.p.),$$

$$H_{\mathrm{str},8v8v}^{(2)}$$

$$+a_{\mathrm{v}}^{(2c)}({ϵ}_{xy}^2+c.p.)$$

Full text

I Introduction

The realm of nanostructures is a very broad family of systems varying from colloidal nanocrystals, delerue-book self-assembled bimberg-book ; michler ; michler2017quantum ; tartakovskii_2012 and nanowire quantum dots bjork ; borgstrom-nwqd and quantum dashes dery2004nature ; Khan2014 , nanowires, vls-wagner ; dalacu-selective quantum wells, single dopants in silicon, Zwanenburg_RMP13 ; Usman2016 and more. One of the key points necessary for understanding of physical properties of nanostructures involves accurate and computationally efficient theoretical studies of their electronic structure. Since calculations for nanostructures may involve multi-million computational boxes, nanostructures are usually beyond the reach of current ab-initio methods. szabo-ostlund ; siesta Therefore, various semi-empirical approaches are typically employed. Moreover, usually researchers opt for one of two apparently opposite schemes. One group of methods is based on the assumption of underlying continuous media. This set involves in particular one-band effective mass approximation, as well as different flavours of multi-band (e.g. 4-, 8- or 14-) $\bm{k}\!\dotproduct\!\bm{p}$ method. stier-grundmann-bimberg ; kp-schliwa ; Schliwa111 ; Schulz111 ; kp14 These approaches combine very high-computational efficiency with unambiguous parameter sets, that can be obtained directly from bulk properties. Moreover, the $\bm{k}\!\dotproduct\!\bm{p}$ approaches turned out to be highly successful in understanding of main spectral pryor-zunger ; wojs ; weidong-localization ; sheng-cheng-prb2005 ; ardelt16 ; climente08 ; leger ; marynski2013electronic ; gawelczyk2017exciton and spin-relatedandlauer09 ; Gawarecki2018a ; Gawelczyk2018a features for a broad group of nanostructures.

There is however a second group of widely acknowledged, so-called ‘‘atomistic‘‘ methods based on explicit accounting for low, atomistic symmetry of nanostructures, which are defined in calculations atom-by-atom. Notably, these 111We note that there are some intermediate schemes in particular the EBOM method, which however does not have atomistic spatial resolution. include empirical pseudopotential method, Zunger-EPM2 ; bester-zunger-nair ; bester-zunger-p-shell ; singh-bester-ordering ; singh-bester-lower and empirical tight-binding method. slater-koster ; vogl ; jancu ; santoprete03 ; saito-arakawa ; schultz12 ; schulz-schum-czycholl ; alloyed-qd ; boykin-fit ; voisin-fit ; jaskolski06 ; jaskolski04 ; zielinski12 ; Zielinski2013 Atomistic approaches often provide more accurate results (such as e.g. the magnitude of the bright exciton splitting bayer-eh ; takagahara ; Swiderski2017 ; Zielinski2013 ) and are methods of choice when one deals with effect of alloying, rand-mlinar-zunger ; singh-bester-ordering ; singh-bester-lower ; Zielinski-natural ; Swiderski2017 ; zielinski12 ; dash-mrowinski low-shape symmetry and faceting, bryant-prl ; bryant-prb ; Zielinski.PRB.2015 ; wojnar2014strain monolayer-thin sizes Diaz.PhysRevB.75.245433 or sections of semiconductor, zhang ; nika-cpqd ; our-cpqd and atomistically sharp interfaces. Zwanenburg_RMP13 ; maaike-strain ; nw-sharp Atomistic methods are however usually much more computationally demanding, alloyed-qd ; zielinski-multiscale ; rozanski-zielinski ; rozanski-zielinski-cpc with often complicated and non-trivial schemes for semi-empirical parameters fitting. alloyed-qd ; boykin-fit ; voisin-fit ; tan-fit

The properties of semiconductor nanostructures are strongly modified by the presence of strain, which is inevitable for a system composed of lattice mismatched materials. In atomistic methods, strain is accounted for by displacements of the atomic positions which alters the bond lengths and bond angles. On the other hand, in continuous media approaches, strain is represented by macroscopic strain tensor field. In a standard way, strain enters the $\bm{k}\!\dotproduct\!\bm{p}$ model via well-established Bir-Pikus Hamiltonian, containing strain-tensor elements in linear orderBir1974 ; Trebin1979 ; bahder90 . This approach accurately describes the band structure, if strain is low or moderate. To describe nonlinear effects (visible in DFT simulations for InAs and GaAs materialsKadantsev2012 ) which are relevant at stronger strain, higher order strain terms would be needed. Despite the second order scheme was proposedSuzuki1974 no reliable parameterization is available.

In this work we aim to bridge between these two seemingly excluding ways of calculations and improve $\bm{k}\!\dotproduct\!\bm{p}$ results considerably by accounting for second-order strain effects. We thus emphasize and study the importance of nonlinear strain effects on the spectral properties of nanostructures. To this end, we go beyond the standard Bir-Pikus Hamiltonian and implement quadratic strain terms with the second order deformation potentialsSuzuki1974 . By comparison to the $sp^{3}d^{5}s^{*}$ tight binding method, we find the relevant parametrization for InAs and GaAs bulk materials. Then, using $8$-band $\bm{k}\!\dotproduct\!\bm{p}$ Hamiltonian (supplemented by the 2nd order strain terms), we calculate the electron and hole energy levels in classes of quantum well (QW) and quantum dot (QD) structures. We show, that single particle energies are strongly affected by nonlinear terms related to the biaxial strain. We also demonstrate that, in the case of electron in the QW and QD, an excellent agreement between the $8$-band $\bm{k}\!\dotproduct\!\bm{p}$ and the $sp^{3}d^{5}s^{*}$ tight-binding model can be achieved if second order strain terms are taken into account. These findings are not only important from the fundamental point of view, but also to properly predict on the the optoelectronic or nanophotonic device characteristics, especially when calculating the excitonic properties in single quantum dots of this material system broadly considered as efficient sources of single or entangled photon states also at the telecommunication spectral range of 1.3 - 1.55 $\mu$m.

The paper is organized as follows. In Sec. II, we describe the models, calculate the band structure, and fit the material parameters used in further simulations. In Sec. III, we calculate the electron and hole states in the QW and the QD structures. Finally, Sec. IV contains concluding remarks. Furthermore, the derivation of strain Hamiltonian is presented in the Appendix.

II Bulk material

The reference band structures for InAs and GaAs are obtained within $sp^{3}d^{5}s^{*}$ tight-binding modelSlater1954 , where we took the material parameters from Ref.Jancu1998 . To account for the strain, we manipulate the atomic positions in the primitive cell. The bond length variations enter the Hamiltonian via the orbital-dependent exponents rescaling the two-center integrals (a generalization of the Harrison law).Jancu1998 ; harrison-big-book Inter-atomic (hopping) matrix elements $t$ are thus modified as:

[TABLE]

where $t_{0}$ ($t$) is unstrained (strained) hopping matrix element, $d_{0}$ ($d$) is unperturbed (perturbed) bond-length,$\eta$ is scaling matrix element whose magnitude is fitted to reproduce bulk deformation potentials. Jancu1998 ; boykin-fit While the hydrostatic strain alters only the bond lengths; the uniaxial, biaxial and shear strain change the bond angles as well. Furthermore, the uniaxial and biaxial strain lead to the energy splitting in the atomic $d$ shell. The $d_{xy}$, $d_{yz}$ and $d_{zx}$ orbital on-site matrix elements are then given byJancu1998

[TABLE]

where $E_{yz}$, $E_{zx}$ are given by cyclic permutations of indices; $E_{d}$ represents the unperturbed $d$-shell energy, $b_{d}$ is a material-dependent parameter, and $\epsilon_{ij}$ are the strain tensor components. We implemented the matrix elements using the form described in Ref.zielinski12 , where the strain tensor components are expressed in terms of the bond lengths and direction cosines. We note here that strain related power-dependence of hopping matrix elements, as accounted by Harrison law, is inherently non-linear and may lead to non-linear bulk bands evolution under external strain, in particular for larger strain as present in InGaAs nanostructures.

The second approach, which we utilize in the paper is the $8$-band $\bm{k}\!\dotproduct\!\bm{p}$ model, where the Hamiltonian explicitly contains $\Gamma_{\mathrm{6c}},\Gamma_{\mathrm{8v}}$ and $\Gamma_{\mathrm{7v}}$ blocks corresponding to the irreducible representations of the $T_{d}$ symmetry point groupmayer91 ; Winkler2003 . According to the invariant expansion schemeLuttinger1955 , the Hamiltonian can be expressed in terms of the matrices representing crystal symmetry invariants. Then, the kinetic part of the Hamiltonian is writtenWinkler2003

[TABLE]

where $E_{g}$ is the energy gap, $\Delta_{0}$ describes the spin-orbit coupling, $P$ is a parameter proportional to the interband momentum matrix element, $m_{0}$ is the free electron mass, $A^{\prime}_{\mathrm{c}}$ accounts for the remote bands contributions to the electron effective mass, $\gamma^{\prime}_{i}$ are modified Luttiner parameters, $\kappa^{\prime}=-\frac{1}{3}(\gamma^{\prime}_{1}-2\gamma^{\prime}_{2}-3\gamma^{\prime}_{3}+2)$, $q^{\prime}$ is an anisotropy parameter, $B_{\mathrm{8v}}^{\pm},B_{\mathrm{7v}}$ are parameters related to the Dresselhaus spin-orbit coupling, $\{A,B\}=\frac{1}{2}(AB+BA)$, ,,${\mathrm{c.p.}}$” denotes cyclic permutations, $\sigma_{i}$ are the Pauli matrices, matrices $J_{i}$ are related to the $j=3/2$ representation of angular momentum, $T_{i}$ are matrices connecting the $j=1/2$ representation to the $j=3/2$, $T_{ij}=T_{i}J_{j}+T_{j}J_{i}$, $U_{i}=T^{\dagger}_{i}$, and $U_{ij}=T^{\dagger}_{ij}$. The explicit definitions of the matrices are given in Refs.mayer91, ; Winkler2003, ; Trebin1979, .

Within the standard Bir-Pikus model, strain tensor elements enters the Hamiltonian in linear orderWinkler2003 ; bahder90 . At $\bm{k}=0$, this is given by

[TABLE]

where $a^{\mathrm{(1)}}_{\mathrm{c}}$, $a^{\mathrm{(1)}}_{\mathrm{v}}$, $b^{\mathrm{(1)}}_{\mathrm{v}}$ and $d^{\mathrm{(1)}}_{\mathrm{v}}$ are the deformation potentialsbir61 ; bahder90 . The term proportional to the $C_{2}$ parameter results from the inversion asymmetryTrebin1979 , and provides a channel of spin-orbit coupling.

We go beyond the standard framework and include the terms which are quadratic in strain tensor elementsSuzuki1974 and correspond to the hydrostatic and biaxial strain, while second order terms with shear strain components are neglected. The relevant part of the Hamiltonian in the invariant expansion is then given by

[TABLE]

where six additional deformation potentials are introduced. A detailed derivation is presented in the Appendix.

We calculated band edges of the InAs and GaAs bulk materials using the tight-binding method. Then, the deformation potentials of $8$-band $\bm{k}\!\dotproduct\!\bm{p}$ model are fitted to the TB data (see Table 1 for values). Fig. 1(a) presents the results of band edges in InAs as a function of hydrostatic strain. In this case, the effect of nonlinearity is weak, and the $H^{\mathrm{(1)}}_{\mathrm{str}}$ already provides a good approximation. The results for uniaxial- and biaxial strain are shown in Figs. 1(b,c)]. For the biaxial strain we use the Poisson ratio $\epsilon_{zz}=-(2C_{12}/C_{11})\epsilon_{\parallel}$, where $C_{12}$ and $C_{11}$ are the elastic constants. The mismatch between the Bir-Pikus model and tight-binding results becomes significant, however a good agreement can be achieved, if the second order strain terms ($H^{\mathrm{(2)}}_{\mathrm{str}}$) are included.

In the case of hydrostatic and uni-/biaxial strain (represented by macroscopic tensor $\hat{\epsilon}$), the displacements of atoms in the unit cell are well defined. On the other hand, if the shear strain is taken into account, the situation becomes more complicated. In such case, the atomic positions are also affected by microscopic relative displacements between the two crystal sublattices. The atomic positions ($\bm{R}_{i}$) are then given byKleinman1962 ; Yu2010 ; Niquet2009

[TABLE]

where ($\bm{R^{0}}_{i}$) is the initial (in unstrained crystal) position of the $i$-th atom, $\zeta$ is the Kleinmann parameter, $a$ is the lattice constant, and the sign $\pm$ depends on the sublattice (anion or cation). In general, the value of Kleinman parameter depends on hydrostatic and shear strain in the systemGarg2009 . Due to uncertainty related to $\zeta$ and the number of free parameters (see the Appendix A), we skipped fitting of the second order shear strain terms. We approximate the value of Kleinman parameter by a constant value with the formulaSteiger2011

[TABLE]

where $\alpha$, $\beta$ are parameters of the Keating VFF modelpryor98b . The comparison between the methods is given in Fig.1(d). As the fitting of $C_{2}$ for GaAs requires high shear strain, its value is not estimated satisfactory (see Table 1). Instead, following Ref.Mielnik-Pyszczorski2018b we take its value as $3.3$ eV, which was extracted from the experimental data of spin relaxation timedyakonov86a . Moreover we note that $sp^{3}d^{5}s^{*}$ parametrization by Jancu Jancu1998 has limited accuracy for representing shear strains, Jancu2007 and a such may not be the best target for $\bm{k}\!\dotproduct\!\bm{p}$ fitting.

In the next step, we fit the band structure parameters for unstrained crystal [see Table 2] for values]. While the TB parametrization Jancu1998 for GaAs is in very good agreement with Ref.vurgaftman01 , the comparison for InAs requires some rescaling of the Luttinger parameters. Moreover, following Ref.birner11 , to avoid spurious solutions in further calculations in nanostructure, we set $A^{\prime}_{c}=1$ and then rescale $E_{\mathrm{p}}$.

In the case of the (single band) effective mass model, the most important second order strain terms enter the Hamiltonian with a single effective parameter. If we neglect relatively small nonlinearity of the hydrostatic strain, the effective mass Hamiltonian for the electron can be written as

[TABLE]

where $m^{\prime}$ is the electron effective mass, $E_{\mathrm{c}}$ is the conduction band (cb.) edge, $a_{\mathrm{c}}$ is the standard cb. deformation potential, $\widetilde{a}^{(2)}_{\mathrm{c}}$ is the effective parameter related to the biaxial strain in the second order. The values of $\widetilde{a}^{(2)}_{\mathrm{c}}$ for InAs and GaAs, (where the fitting procedure was optimized for the negative biaxial strain) are given in Table 1.

III Nanostructures

In this Section, we utilize the models described previously for the calculations of carrier states in QWs and QDs. We find the strain distribution in the system within standard Keating VFF modelpryor98b . The elastic energy is given by

[TABLE]

where $\mathrm{NN}(i)$ denotes nearest neighbors of $i$-th atom, $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ and $d_{ij}$ are actual and idealized (unrelaxed) distances between $i$-th and $j$-th atom, $A_{ij}=\alpha_{ij}/d_{ij}^{2}$ with bond-streching constant $\alpha_{ij}$, $B_{ijk}=(\beta_{ij}+\beta_{ik})/(2d_{ij}d_{ik})$, here $\beta_{ij}$ is a constant which represents bond bending. We use PETSC TAOpetsc library to minimize the elastic energy of the system via relaxation of the atomic positions. The strain tensor elements are obtained from the direction cosines and the bond lengths, in a way described in Ref. zielinski12 . We find the electron and hole states within multimilion atom simulations using the $sp^{3}d^{5}s^{*}$ tight-binding modelSlater1954 ; zielinski12 ; Jancu1998 . The passivation of dangling-bonds at the limit of computational domain, has been performed following Ref.Lee2004 .

For the 8-band $\bm{k}\!\dotproduct\!\bm{p}$ simulations, we move from the atomic lattice to the cartesian grid by taking $h_{x}=h_{y}=a_{\mathrm{GaAs}}$, $h_{z}=a_{\mathrm{GaAs}}/2$ mesh size, and averaging the strain and composition over the two cations in each mesh cell. To calculate electron and hole states in a nanostructure within the 8-band $\bm{k}\!\dotproduct\!\bm{p}$ model, we perform the standard substitution $k_{i}=-i\hbar\partialderivative{x_{i}}$ in the bulk Hamiltonian. Since all of the material parameters are position dependent, $k_{i}$ does not commute with them and the operator ordering becomes important. The details of the implementation are described in the Appendix of Ref.Gawarecki2018a .

III.1 Quantum Wells

In the previous section, we found that the non-linear, 8-band $\bm{k}\!\dotproduct\!\bm{p}$ model is capable to reproduce accurately bulk band edges evolution under strain, as given by the atomistic tight-binding model. In the following, we will verify how both methods compare for nanostructures, where the quantum confinement plays a significant role. We start by calculating the ground state energy of the electron [Fig.2(a)] and hole [Fig.2(b)], in a quasi-two dimensional system, namely the InAs quantum well embedded in bulk GaAs. We study electron and hole ground state energies as a function of quantum well width changing from $5a$ to $25a$, where $a$ is the lattice constant. This corresponds to quantum well thickness from approximately $3$ to $15$ nm. The computational domain size for the strain calculation is taken $220a\times 220a\times 220a$, while the single particle states are calculated in a $80a\times 80a\times 60a$ box. zielinski-multiscale Fig. 2 shows results obtained by the tight binding calculation, the standard Bir-Pikus model, and the Bir-Pikus model with non-linear terms. Apparently (and similarly to the strained bulk results) the linear Bir-Pikus model significantly overestimates the single particle electron energy and as well it underestimates the ground hole state energy (please note reversed ordering of hole states as compared to the electron). This is however expected since the quantum well is under high biaxial strain (the lattice mismatch between GaAs and InAs is about $7$%), which according to Fig. 1(c) has strongly non-linear character and is not well reproduced by linear $\bm{k}\!\dotproduct\!\bm{p}$ approach already at the bulk level. In particular, in case of linear $\bm{k}\!\dotproduct\!\bm{p}$ the electron ground state energy $e_{1}$ is overestimated by approximately $100$ meV with respect to the tight-binding approach. Additionally, approximately $30$ meV difference is also present in the energy of the hole ground state $h_{1}$. Overall, the linear $\bm{k}\!\dotproduct\!\bm{p}$ systematically overestimates the single particle energy gap $e_{1}-h_{1}$ by approximately $130$ meV. On the other hand, we obtain an excellent agreement between the 8-band $\bm{k}\!\dotproduct\!\bm{p}$ and the tight-binding method, if the terms second order in strain are taken into account. For all considered quantum well thicknesses we obtain at most $\approx 10$ meV difference between both methods, with only some variations due to different quantum well heights.

III.2 Quantum Dots

Disk-shaped quantum dots

Next, we have calculated several lowest electron and hole states for a series of disk-shaped InAs/GaAs quantum dots. Here, we have assumed fixed quantum dot basis diameter equal to approximately $25$ nm and systematically varied quantum dot height from approximately $3$ to $5.4$ nm ($5$ to $9$ $a$). The dot is placed on a lattice constant thick ($\approx 0.6$ nm) wetting layer. Similarly to results for quantum wells, for disk-shaped quantum dots, if the second order strain terms are taken into account, we obtain an excellent agreement for the electron ground state [Fig.3(a)]. Additionally, we study here the energy difference [Fig.3(c)] between the first excited electron and the ground electron states (i.e. energy difference between electron $p-$ and $s-$ shells, i.e. $e_{2}-e_{1}$). This energy difference increases with the quantum dot height, and the magnitude of this spacing is also much better (within few meV error) reproduced by the augmented, non-linear $\bm{k}\!\dotproduct\!\bm{p}$ approach as compared to the straightforward linear $\bm{k}\!\dotproduct\!\bm{p}$, with predictions systematically smaller by approximately $10$ meV.

In the case of the hole ground state energy [Figs.3(b)], the relative discrepancies between 8-band $\bm{k}\!\dotproduct\!\bm{p}$ and tight-binding model become somewhat larger. This could be expected since Jancu‘s Jancu1998 tight-binding model is a 20-band, atomistic model with d-orbitals, and hole states in quantum dots have complicated multi-band character. However, interestingly all three models predict quite similar energy of the ground hole state ($\approx 0.3$ eV), with non-linear results only several meV‘s closer to the tight-binding results for flat (small height) quantum dots, than those given by linear $\bm{k}\!\dotproduct\!\bm{p}$. Moreover, contrarily to a quantum well, the ground hole state energy of a disc-shape quantum dot rather weakly depends on quantum dot height. Actually, it even shows an opposite trend with respect to the height, due to strain distribution changing its character with an increasing quantum dot height. zielinski-tallqd Therefore, size-dependent strain distribution in finite size quantum dots (rather than in quasi-2D quantum well) seems to play here a dominant role. Should strain effects be artificially neglected (see the Appendix B) the ground hole state energy would simply increase with dot‘s height due to reduced confinement. In fact, as shown in the Appendix, $\bm{k}\!\dotproduct\!\bm{p}$ and tight-binding results for disk-shaped quantum dots with strain artificially neglected are quite similar, again demonstrating the key role of strain in modeling of quantum dots, and in particular the role of non-linear strain treatment in $\bm{k}\!\dotproduct\!\bm{p}$.

To summarize, the single particle gap ($e_{1}-h_{1}$) in disk-shaped InAs/GaAs quantum dots is synthetically overestimated by $\approx 100$ meV by the linear $\bm{k}\!\dotproduct\!\bm{p}$ approach (mostly due to large overestimation of $e_{1}$), and only by $\approx 15$ meV by the non-linear $\bm{k}\!\dotproduct\!\bm{p}$ method (mostly due to small underestimation of $h_{1}$).

The energy difference between the ground and the first excited hole states ($s-$ and $p-$ shells spacing, i.e., $h_{1}-h_{2}$,222Note reversed ordering of hole states with respect to the electron) increases with the quantum dot height as in the case of electron [Figs.3(d)], again with non-linear $\bm{k}\!\dotproduct\!\bm{p}$ predictions closer to the tight-binding, than the linear approach results. The second order $\bm{k}\!\dotproduct\!\bm{p}$ overestimates the magnitude of this spacing by $\approx 2$ meV, whereas the linear $\bm{k}\!\dotproduct\!\bm{p}$ underestimates this difference by $\approx 4$ meV. Therefore, both for electron and hole $s-p$ shell spacings, the agreement between the non-linear $\bm{k}\!\dotproduct\!\bm{p}$ and the tight-binding could be described as satisfactory considering how inherently different (continuous media vs atomistic) $\bm{k}\!\dotproduct\!\bm{p}$ and tight-binding models are.

Lens-shaped quantum dots

It is instructive to verify, how non-linear $\bm{k}\!\dotproduct\!\bm{p}$ performs for other geometries, rather than flat quantum wells and disk-shaped quantum dots. Therefore, Fig.4 shows results obtained, for lens-shaped quantum dot, with its height varying from approximately $3$ to $\approx 6.6$ nm ($5$ to $11$ $a$), while keeping the base diameter fixed to about $24$ nm. Similarly to the previous case, the lens-shaped quantum dot is placed on a lattice constant thick ($\approx 0.6$ nm) wetting layer. Again, for the ground electron state we obtain an excellent agreement [Fig.4(a)] between the non-linear $\bm{k}\!\dotproduct\!\bm{p}$ and the tight-binding.

Again, a very good agreement [Fig.4(b)] is found for $p-s$ electron energy spacing (i.e., $e_{2}-e_{1}$), with the non-linear method differing by at most $1$ meV from the tight-binding. Notably, for the lens-shaped quantum dot (and differently from the disk-shaped quantum dots) the electron $p-s$ energy difference is reduced with the quantum dot height.

Regarding, the ground hole state [Fig.4(b)], similarly to the disk-shape quantum dots, both $\bm{k}\!\dotproduct\!\bm{p}$ variants underestimate this energy by about $20$ to $30$ meV. Thus, the single particle band gap ($e_{1}-h_{1}$) is overestimated by $\approx 70$ meV by the linear approach, and only by at most $\approx 30$ meV by the non-linear approach. Again, should the strain effect be artificially neglected (see the Appendix B), the difference between 8-band $\bm{k}\!\dotproduct\!\bm{p}$ and the tight-binding predicted ground hole energies would be withing several meV‘s.

A more notable difference between various methods is observed for the $s-p$ energy spacing of hole states ($h_{2}-h_{1}$) which is systematically overestimated ($\approx 4$ meV) by both $\bm{k}\!\dotproduct\!\bm{p}$ approaches. The difference between atomistic and continuous media approximation is thus somewhat larger for hole states in lens-shaped quantum dot than in disk-shaped quantum dots. This difference is present for both unstrained (the Appendix) and strained quantum dots, yet strain further increases its magnitude.

Part of this discrepancy is probably related to shear strains which are more notable in lens-shaped quantum dots due to their curved shape. Deformation potentials due to shear strains are in fact reproduced with limited accuracy by Jancu‘s tight-binding parameterization. Jancu1998 ; Jancu2007 Having said the above, we should note that is not our goal to claim that the continuous media approximation (with complicated atoms to grid strain transfer) in all considered cases is able reproduce the atomistic results.

IV Conclusions

We have investigated the influence of strain nonlinearity on the electron and hole energy levels in semiconductor nanostructures. We have used Bir-Pikus Hamiltonian extended by the terms which are quadratic in strain. We have obtained the set of parameters for InAs and GaAs bulk semiconductors by fitting to the results obtained from the $sp^{3}d^{5}s^{*}$ tight-binding method. Next, we have provided the ready-to-use parameterization for the effective term used to improve the accuracy of the standard effective mass method. Then, we have calculated the electron and hole states for quantum wells, and various quantum dot systems. We have shown that, while the standard (linear in strain) $8$-band $\bm{k}\!\dotproduct\!\bm{p}$ model overestimates the electron and underestimates the hole energy, a very good agreement to the tight-binding results can be achieved if the terms second order in strain are taken into account. This agreement is particularly good for electron states confined in quantum wells and flat quantum dot systems. More discrepancies are found for hole states, especially in curved, lens-shaped quantum dots. We are convinced this differences origin from the different treatment of strains (in particular shear strain) by $\bm{k}\!\dotproduct\!\bm{p}$ and tight-binding, with a further research needed to tract these divergences.

Acknowledgements.

K.G. acknowledges the support from the National Centre for Research and Development in Poland within the FI-SEQUR Project of the 2nd Poland-Berlin Photonics Programme, grant No. 2/POLBER-2/2016 (project value 2 089 498 PLN) and from the Polish National Agency for Academic Exchange. M.Z. acknowledges the support from the Polish National Science Centre based on Decision No. 2018/31/B/ST3/01415. Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), Grant No. 203.

Appendix A Second order strain Hamiltonian

According to the standard $\bm{k}\!\dotproduct\!\bm{p}$ model, the Hamiltonian can be written byWinkler2003

[TABLE]

where $m_{0}$ is the free electron mass, $V_{0}(\bm{r})$ describes the crystal potential, $c$ is the speed of light, and $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ denotes vector of Pauli matrices.

The strain changes crystal lattice and lowers symmetry of the system. The Hamiltonian of the deformed crystal $H^{\prime}_{\bm{k}^{\prime}}$ can be described in terms of the transformed coordinates $\bm{r}^{\prime}$ Bir1974 . Then, the Hamiltonian is written in basis of $\bm{k}=0$ Bloch functions $\ket{\mu\lambda}$, where $\mu$ describes its orbital part and $\lambda$ is a spinWinkler2003 . The matrix elements, expressed up to the second order in strain tensor elements, takes the form

[TABLE]

where $E_{\nu}$ describes the unstrained band edges, $D^{\nu\mu}_{ij}$ and $F^{\nu\mu}_{ijkl}$ are first- and second-order deformation potentials respectively.

With the group theory, one can predict the number of independent tensor compontents ($N_{0}$) for a given symmetry point group. If tensor $\bm{S}$ connects $\bm{A}$ and $\bm{B}$ (eg. $A_{i}=\sum_{j}S_{ij}B_{j}$), and $\bm{A}$, $\bm{B}$ belong to the $\mathcal{D}_{A}$, $\mathcal{D}_{B}$ representations respectively, the number of independent $\bm{S}$ components is given by the number of trivial representations in Bir1974

[TABLE]

Let us consider the tensor describing deformation potentials ($\hat{D}$) as $\bm{S}$, a given block of the Hamiltonian as $\bm{A}$, and the strain tensor as $\bm{B}$. The strain tensor $\hat{\epsilon}$ belongs to the representation $\Gamma_{5}\otimes\Gamma_{5}=\Gamma_{1}\oplus\Gamma_{3}\oplus\Gamma_{4}\oplus\Gamma_{5}$, and since it is symmetric, the part related to $\Gamma_{4}$ vanishes.Bir1974 ; Yu2010 If we neglect the spin-orbit coupling, the conduction band (cb.) is related to the $\Gamma_{1}\otimes\Gamma_{1}$, and the valence band (vb.) to the $\Gamma_{5}\otimes\Gamma_{5}$ representations. Hence, in the case of the $D^{\nu\mu}_{ij}$, the number of the independent components is: one for the cb., and three for the vb. block. In consequence, only four deformation potentials are needed: $a_{\mathrm{c}}$, $a_{\mathrm{v}}$, $b_{\mathrm{v}}$, and $d_{\mathrm{v}}$ (we do not consider here, the block connecting cb. with vb.).

In the case the strain quadratic terms, the product $\hat{\epsilon}\otimes\hat{\epsilon}$ belongs to $(\Gamma_{5}\otimes\Gamma_{5})\otimes(\Gamma_{5}\otimes\Gamma_{5})=3\Gamma_{1}\oplus 3\Gamma_{3}\oplus 3\Gamma_{5}$, which leads to $3$ independent terms in the cb. block ($a^{\mathrm{(2a)}}_{\mathrm{c}}$, $a^{\mathrm{(2b)}}_{\mathrm{c}}$, $a^{\mathrm{(2c)}}_{\mathrm{c}}$) and $9$ components in the vb. block (which we denote as $a^{\mathrm{(2a)}}_{\mathrm{v}},a^{\mathrm{(2b)}}_{\mathrm{v}},a^{\mathrm{(2c)}}_{\mathrm{v}},b^{\mathrm{(2a)}}_{\mathrm{v}},b^{\mathrm{(2b)}}_{\mathrm{v}},b^{\mathrm{(2c)}}_{\mathrm{v}},d^{\mathrm{(2a)}}_{\mathrm{v}},d^{\mathrm{(2b)}}_{\mathrm{v}},d^{\mathrm{(2c)}}_{\mathrm{v}}$). In presence of spin-orbit interaction, the Hamiltonian blocks are described in terms of the double group representation, but it does not change the number of independent tensor components related to strain.

We write the Hamiltonian using the invariant expansion methodLuttinger1955 . We note, that strain tensor elements transform like symmetrized products $\{k_{i},k_{j}\}$. Then,we utilize the list of irreducible tensor components for the point group $T_{d}$, given (up to the 4th order in $k$) in the Appendix C of Ref.Winkler2003 . This leads to the form

[TABLE]

We neglected quadratic strain in off-diagonal blocks between the cb. and valence bands (namely $H^{\mathrm{(2)}}_{\mathrm{str},\mathrm{6c8v}}$ and $H^{\mathrm{(2)}}_{\mathrm{str},\mathrm{6c7v}}$). Due to large number of free parameters and uncertainty in the value of Kleinman parameter ($\zeta$), we skipped the fitting of parameters describing the shear strain ($a^{\mathrm{(2c)}}_{\mathrm{c}}$, $a^{\mathrm{(2c)}}_{\mathrm{v}}$, $b^{\mathrm{(2c)}}_{\mathrm{v}}$, $d^{\mathrm{(2a)}}_{\mathrm{v}}$, $d^{\mathrm{(2b)}}_{\mathrm{v}}$, and $d^{\mathrm{(2c)}}_{\mathrm{v}}$). We also neglected terms containing products of wave vector and strain tensor elements.

Appendix B Quantum dots with strain effects neglected

Here, for completeness, we show results obtained for disk-shaped [Fig. 5] and lens-shaped [Fig. 6] quantum dots of the same shapes, dimensions, and compositions, as in the main text, yet with strain effects artificially neglected (i.e. assuming identical bond lengths of InAs and GaAs, as well as, perfect tetrahedral bond angles).

Both figures show the same trends of decreasing (increasing) ground electron (hole) state energy with quantum dot height, consistent with the quantum confinement effect. For both geometries the agreement between $\bm{k}\!\dotproduct\!\bm{p}$ and the tight-binding description of carriers ground state energies is very good. Somewhat smaller agreement is found for the energy difference between the ground and the first excited states. Here, the $\bm{k}\!\dotproduct\!\bm{p}$ results systematically overestimate the magnitude of this splitting, with respect to the tight-binding, for both disk- and lens-shaped quantum dots, and both for the electron and the hole.

In the unstrained case, the hole $s-p$ shells spacing, ($h_{1}-h_{2}$ energy spacing) [Fig. 6 (c)] is reduced with growing quantum dot height for both disk- and lens-shaped systems. However, interestingly, the electron $s-p$ spacing either increases with the quantum dot height for disk-shaped nanostructures or show a rather flat (and non-monotonous) trend for lens-shaped systems. In all considered cases, contrary to a situation where the strain is accounted for, $\bm{k}\!\dotproduct\!\bm{p}$ and tight-binding results are both qualitatively and quantitatively similar.