Efficient quantum dot \mathbf{k}\cdot\mathbf{p} in wurtzite systems including spatially varying elastic and dielectric constants and smooth alloy profile
EEfficient quantum dot k · p in wurtzite systems including spatially varying elastic anddielectric constants and smooth alloy profile Luc Robichaud and Jacob J. Krich Department of Physics, University of Ottawa, Ottawa, Canada K1N 6N5 and School of Electrical Engineering and Computer Science,University of Ottawa, Ottawa, Canada K1N 6N5 (Dated: 16 December 2020)We present Fourier-space based methods to calculate the electronic structure of wurtzite quantumdot systems with continuous alloy profiles. We incorporate spatially varying elastic and dielectricconstants in strain and piezoelectric potential calculations. A method to incorporate smooth alloyprofiles in all aspects of the calculations is presented. We demonstrate our methodology for the caseof a 1D InGaN quantum dot array and show the importance of including these spatially varyingparameters in the modeling of devices. We demonstrate that convergence of the lowest bound stateenergies is to good approximation determined by the largest wave vector used in constructing thestates. We also present a novel approach of coupling strain into the k · p Hamiltonian, greatlyreducing the computational cost of generating the Hamiltonian.
I. INTRODUCTION
Given their large range of bandgaps, from 0.78 eV to 3.51 eV, InGaN materials have attracted attention fromvarious applications such as LEDs, single-photon emitters, water splitting and solar cells [1–5]. For any application,device performance depends on having an electronic structure well tuned to its target application. Given that theelectronic structure of quantum dots can be drastically changed by varying their size and composition, they can bequite attractive for applications. The main problem in modeling complex structures such as quantum dots is includingall the necessary effects for the model to be accurate while also keeping computational cost down.Tight binding and k · p theory are standard approaches for calculating single-particle electronic structures forbulk materials and nanostructures [6]. The k · p method gives a good balance between accuracy and computationalrequirements, especially when considering large dots that contain large number of atoms where the tight bindingmethod becomes costly. k · p theory has been developed in both real space and Fourier space [7, 8]. Following theFourier-space method, symmetry adapted basis approaches have been developed to reduce the required size of theHamiltonian, which block diagonalize the Hamiltonian, reducing the computational cost of calculating the system’seigenstates [9–11].InGaN materials are strongly piezoelectric, having both spontaneous and strain-induced contributions to the piezo-electric polarization. Strain calculations have been performed using valence force field and Green’s function basedmethods [8, 12]. The latter method has the advantage that it respects the symmetry of the crystal lattice. Fromthe strain, the piezoelectric potential can be calculated from Maxwell’s equations [8]. References [8, 10, 13] use theGreen’s function method for calculating strain and calculated the piezoelectric potential from Maxwell’s equations.These works assume uniform elastic and dielectric constants, which was justified for their respective InAs/GaAs andGaN/AlN systems. However, in the case of InGaN, these constants vary more significantly between dot and host.Additionally, InGaN devices frequently do not have sharp interfaces between dot and barrier, with indium diffusingover several nanometers. This smooth alloy profile gives a spatial profile to every material parameter of the system,effectively changing the confining potential seen by the electrons.In this paper, we show the importance of including spatially varying elastic and dielectric constants in strain andpiezoelectric potential calculations in the case InGaN systems. For strain calculations, we implement a formalismpreviously presented for including spatially varying elastic constants [8]. We present a new Fourier-space formalism forthe calculation of piezoelectric potentials with spatially varying dielectric constants. We also present an approach toinclude smooth indium profiles in the strain, piezoelectric potential and electronic structure calculations, modeling thesmooth alloy profiles found in experimental devices. Considering smooth indium profiles both increase the accuracyof the simulations and decrease their computational cost by decreasing the number of plane waves required forconvergence.Strain plays an important role in the electronic structure properties of quantum dots. In quantum dot k · p , a singlereal space unit cell is typically used when working in a Fourier-space approach. However, strain decays more slowlythan bound state wavefunctions. When studying isolated dots, the difference in decay lengths makes it computationallyexpensive to fully capture both the strain and electronic structure using a single unit cell. Reference [11] presents anapproach that implements two different unit cells; one for the electronic structure and on for strain. This methodallows for the modeling of the electronic structure and strain, but introduces some complexity in calculating the a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Hamiltonian, which requires the calculation of multiple composed convolutions on different Fourier-space meshes.These convolutions can be computationally costly depending on the sizes of meshes needed for convergence. By fixingthe strain unit cell to be commensurate with the electronic unit cell, we present an approach that reduces the numberof needed convolutions, significantly reducing the computational cost.We demonstrate our methodology by calculating the electronic structure for a 1D array of InGaN quantum dots,modeling devices grown as LEDs and for water splitting [1, 3]. In this example, we show the importance of the inclusionof spatially varying elastic and dielectric constants and smooth indium profiles for accurate electronic structures. Wealso show that the most important criterion for convergence of the lowest quantum dot electron and hole energies isthe maximum wave vector included in the Fourier-space sampling, which can be increased with low computationalcost by using a small unit cell.Section II contains strain and piezoelectric potential calculations using spatially varying elastic and dielectric pa-rameters. Section III presents the k · p model used for electronic structure calculations and our novel approach toefficiently include strain through choices of unit cells. Section IV introduces a method to use smooth indium profilesin all aspects of our calculations. Section V demonstrates our entire methodology for the case of a 1D quantum dotarray, such as quantum dots grown inside of nanowires [1]. II. SPATIALLY VARYING ELASTIC AND PIEZOELECTRIC CONSTANTS CORRECTIONS
We begin by considering quantum dot heterostructures with abrupt changes in alloy fraction. Alloying the hostmaterial changes the local lattice constants, leading to a lattice mismatch at the host and dot material boundary.This lattice mismatch is a source of strain throughout the QD system, affecting the electronic states of the system.For example, InN has a larger lattice constant than GaN, so alloying GaN with indium to form quantum dots induceschange in the lattice constant. Additionally, strain can generate strong piezoelectric potentials in materials such asIII-nitrides. The piezoelectric potential in III-nitrides is particularly important along the c-axis and can be strongenough to spatially separate electron and hole states through the quantum-confined Stark effect [14].In prior work, elastic and dielectric constants are largely assumed to be spatially uniform in Fourier-based calcu-lations of strain and the piezoelectric potential. In fact, these material properties are different in the dot and hostmaterials, which can cause significant errors when determining electronic structures. Here, we calculate the strainand piezoelectric potential of a quantum dot superlattice with elastic and dielectric constants that vary with alloyfraction, while focusing on the changes brought on by spatially changing parameters. In the case of the spatiallyvarying elastic constants, we use a method outlined in Ref. [8]. We present a version with typos in Eqs. A3, A7 andA8 of Ref. [8] corrected in Section II B. For the piezoelectric potential, we use a procedure similar to Ref. [8], but weconstruct a theory to include spatially varying dielectric constants. The strain field and piezoelectric potential arecoupled into a k · p model, presented in Section III, for electronic structure calculations. A. Quantum dot system
We consider a superlattice of cylindrical wurtzite quantum dots embedded in a bulk host material, as shown inFig. 1. InGaN quantum dots such as those described in Ref. [1] have a lens-like shape and do not have a sharplydefined boundary. We approximate these quantum dots as being cylindrical. This choice of dot geometry simplifiescalculations, as described in Sec. II B, and preserves the C v symmetry of the material, which we take advantage ofin Section III for electronic structure calculations. Hexagonal periodic boundary conditions are used to also preservethe material’s C v symmetry. For single-dot calculations, the superlattice unit cell must be large enough that thechoice of cell size does not affect results. For actual quantum dot arrays, we consider only hexagonal superlattices inthe plane.In this periodic system, the real space quantum dot superlattice is defined by the set of lattice vectors L i , as shownin Fig. 1. We denote the real space unit cell by Ω e , its volume V e , and the reciprocal-space unit cell by Ω − e . The index“e” indicates that these quantities relate to the electronic cell, as opposed to the strain unit cell, which is introducedin Section III C. Imposing periodic conditions in real space implies a discrete reciprocal space with wave vectors q = i b + i b + i b i , i , i ∈ Z (1) FIG. 1. Quantum dot superlattice and its basis vectors. White regions are the host material and grey regions are thequantum dots. Dashed lines show the unit cell boundaries of the quantum dot superlattice. Due to the symmetry, we have L ≡ | L | = | L | . with the reciprocal basis vectors b = 2 πL [1 , − √ , b = 2 πL [1 , √ , b = 2 πL [0 , , (2)Due to the symmetry of the system, we have L = L , which we define as L . In our reciprocal-space calculations,we sample on sets of the wave vectors q ∈ Ω − e . We define m and m such that i , i = {− m , · · · , , · · · , m } and i = {− m , · · · , , · · · , m } . This sampling produces a hexagonal mesh of size N = N N N where N i = 2 m i + 1 . Toobtain a C symmetric mesh, we remove points such that | q x | > m
12 2 πL , leaving a mesh whose size we denote by N e .By choosing the unit cell dimensions L i large enough, it is possible to remove electronic coupling between neighboringdots. This flexibility allows us to model 3D, 2D and 1D arrays of coupled dots. The isolated dot case can also beobtained by choosing both L and L sufficiently large. Section II B presents a method that also uncouples dots interms of strain, which is based on calculating strain and the electronic structure using different unit cells.We illustrate the methods presented in this manuscript by modeling a quantum dot system inspired by Ref. [1].That system consists of InGaN dots grown in GaN nanowires. We approximate this system as a 1D quantum dotarray, by choosing L to match the measured dot-dot spacing and L large enough to avoid dot-dot coupling. Wefix the dot indium alloy fraction, radius and height based on the experimental device. System parameters are listedin Table I and material parameters are in Appendix A. B. Strain
In this section, we present how we calculate strain with elastic constants that depend on alloy fraction for 3D, 2Dand 1D quantum dot superlattices and isolated dots. Our method follows from Refs. [8, 11].Materials such as InGaN have elastic constants that vary based on the alloy fraction. Therefore, the spatial variationof the elastic constants throughout the superlattice unit cell must be included for accurate calculations of strain. Wepresent a method, originally derived in Ref. [8], to include spatially varying elastic constants in strain calculations.We calculate the strain produced by a single isolated dot and construct the quantum dot superlattice strain by linearsuperposition.The calculated strain is to be coupled into the electronic structure calculations. However, strain decays considerablyslower than bound electronic wavefunctions. In the case of isolated dots, the unit cell must be large enough toaccomodate the strain decay. Choosing a unit cell large enough to capture the strain decay reduces the maximumwave vector attainable when using a fixed number of plane waves. As we demonstrate in Section V, accurately
TABLE I. Quantum dot superlattice system parameters used for calculations unless specified otherwise.Parameter Value X h ˚ AR ˚ AL ˚ AL ˚ Am m n n δ [1 . , . , .
5] ˚ A describing the electronic states requires using sufficiently large wave vectors, and thus a large unit cell requires a largenumber of plane waves. Following Ref. [11], we consider that the electronic model and strain model each have their ownreal space unit cells. This additional degree of freedom allows accurate and computationally efficient determination ofboth electronic structure of rapidly decaying confined quantum dot states and longer-range strain effects in isolateddots. In the case of a quantum dot superlattice, different real-space electronic and strain unit cells are not required.
1. Isolated quantum dot strain
In prior work, lattice-mismatch-driven strain has been calculated for a single dot using a continuum theory with aGreen’s function approach while assuming spatially uniform elastic constants [8, 13, 15]. Here, we present the methodoutlined in Appendix A of [8] to include spatially varying elastic constants. In this section, we show how the spatiallyvarying elastic constants modify strain and how this modified strain changes the piezoelectric potential in SectionII C. We show that the elastic constant correction is necessary to obtain accurate strain and piezoelectric potentials.Consider a single InGaN QD in bulk GaN with spatially varying elastic constants λ ijmn ( r ) that depend on thelocal alloy fraction, λ ijmn ( r ) = λ d ijmn + λ h ijmn [1 − χ d ( r )] , (3)where λ h ijmn and λ d ijmn are the host and dot’s elastic constants, respectively. Assuming spatially varying elasticconstants, the Green’s tensor ˜ G in for the displacement field in an infinite anisotropic elastic medium must satisfy ∂∂x k λ iklm ( r ) ∂∂x m G ( r , r (cid:48) ) = − δ ( r − r (cid:48) ) δ in (4)Taking the Fourier transform of Eq. 4, we obtain λ h iklm q k q m ˜ G ln ( q , r (cid:48) )+ ∆ λ iklm (cid:88) q (cid:48) ˜ χ d ( q − q (cid:48) ) q k q (cid:48) m ˜ G ln ( q (cid:48) , r (cid:48) )= 1(2 π ) e i q · r (cid:48) δ in . The system strain is given by the superposition ˜ (cid:15) lm ( q ) = e T lm ˜ χ d ( q ) + ˜ (cid:15) c lm ( q ) where e T lm is the stress-free strain dueto the initial lattice mismatch and ˜ (cid:15) c lm is the interface-driven strain [8, 13]. Reference [8] showed that the Green’stensor can be related to the strain ˜ (cid:15) c lm ( q ) to obtain λ h iklm q k ˜ (cid:15) clm ( q )+ ∆ λ iklm q k (cid:88) q (cid:48) ˜ χ d (cid:16) q − q (cid:48) (cid:17) ˜ (cid:15) clm (cid:16) q (cid:48) (cid:17) = − λ d ikpr ˜ (cid:15) Tpr q k ˜ χ d ( q ) , (5)with typos fixed, where χ d is the characteristic function of the dot, which is unity inside the dot and zero outside (seeAppendix B for its Fourier transform in our case of cylindrical dots), and e T pr is e T pr = ε a δ pr + ε ca δ p δ r with ε a = (cid:0) a h − a d (cid:1) /a d , ε c = (cid:0) c h − c d (cid:1) /c d and ε ca = ε c − ε a . Here, a h , c h are the lattice constants of the host materialand, a d and c d are of dot material. More specifically, a is the xy-plane lattice constant and c is the lattice constantalong the z-axis. A solution for ˜ (cid:15) c lm ( q ) can be found by expanding ˜ (cid:15) c lm ( q ) in a power series, ˜ (cid:15) c lm ( q ) = ˜ (cid:15) (0) lm ( q ) + ˜ (cid:15) (1) lm ( q ) + ˜ (cid:15) (2) lm ( q ) + · · · , (6)where ˜ (cid:15) ( N ) lm ( q ) ∝ (cid:0) ∆ λλ (cid:1) N , ∆ λ ijmn = λ h ijmn − λ d ijmn and the condition ∆ λλ (cid:28) ensures convergence of the series. Theleading term ˜ (cid:15) (0) lm corresponds to uniform elastic constants of the dot with each subsequent term being a correction toinclude spatial variations due to the alloy profile. Using the Einstein summation convention, each term has the form ˜ (cid:15) ( N ) lm ( q ) = (2 π ) (cid:104) F ( N ) p ( q ) q l ˜ G h mp ( q )+ F ( N ) p ( q ) q m ˜ G h lp ( q ) (cid:105) (7)where F (0) i ( q ) = − λ d ikpr ˜ (cid:15) Tpr q k ˜ χ d ( q ) (8) F ( N ) i ( q ) = − ∆ λ iklm q k (2 π ) V (cid:88) q (cid:48) ˜ χ d (cid:16) q − q (cid:48) (cid:17) ˜ (cid:15) ( N − lm (cid:16) q (cid:48) (cid:17) (9)where Eqs. 7 and 9 are corrected from Ref. [8]. Here, ˜ G h in is the Green’s tensor for the host material and is fullywritten out in Appendix D.It has been shown, when assuming uniform elastic constants, that using the parameters for the host material givesmore accurate results. We compare the strain corrected at various orders according to Eq. 6 to the usually consideredcase of uniform elastic constants of the host material. Figure 2 shows the convergence of the strain corrections for the1D quantum dot array system described in Section II A. We quantify convergence with the following metric for thenorm of the strain: | ˜ (cid:15) | = (cid:115)(cid:88) m (cid:62) l V (2 π ) (cid:90) d q | ˜ (cid:15) lm ( q ) | (10)where m (cid:62) l indicates the sum of the unique elements of the strain tensor ( ˜ (cid:15) , ˜ (cid:15) , ˜ (cid:15) , ˜ (cid:15) , ˜ (cid:15) , ˜ (cid:15) ). The green linein Fig. 2 compares the norm of the corrected strain ˜ (cid:15) , calculated from Eq. 7, to the norm of the strain ˜ (cid:15) GaN , whichis calculated assuming spatially uniform elastic constants of GaN. Blue line shows the self-convergence of the powerseries in Eq. 6. From these results, we conclude that a 2nd order correction is sufficient to have strain convergedwithin 1% in self-convergence and that this converged strain differs from the uniform case by about 6%, indicatingthat the elastic constant corrections are important for accurate strain fields in InGaN systems. In Sec. II C, we showthat the calculated piezoelectric potential remains essentially unchanged from 3rd order corrections and up. Giventhat including these corrections are not computationally costly, we have included 3rd order corrections in all of ourcalculations unless stated otherwise. Figure 3 shows the hydrostatic strain along a cut through the axis of the dot,showing relaxation of strain inside the dot with each additional correction.
2. Quantum dot superlattice strain
Because the strain is linear in stress, the strain produced by the QD superlattice can be obtained from linearsuperposition of the single-dot strain. However, we want the ability to study dots that are completely uncoupled,both electronically and from strains of the periodic array. Reference [11] proposed a method to allow simultaneoustreatment of large unit cell for the strain problem and small unit cell for the electronic problem, which togetherallow isolated dots to be considered in a computationally tractable manner. In this case of two independent cells,
FIG. 2. Convergence of the strain with spatially varying elastic constants. Green line is the relative difference between thecorrected strain and the case with uniform λ of GaN for a 1D quantum dot array as described in Table I. The zeroth term isthe case of uniform elastic constants of InGaN with alloy fraction of the dot. Converged strain rests at a 6% relative differencefrom the case of uniform elastic constants of GaN, indicating that the corrections are necessary for accurate strain fields. Blueline is self-convergence of the power series in Eq. 6 in respect to the zeroth term. Correction magnitudes are found to be lessthan 1% starting from 2nd order. We conclude that corrections up to and including order 3 are sufficient for the calculation ofaccurate strain fields.FIG. 3. Hydrostatic strain for increasing correction orders in the elastic constants. Black dashed line is strain assuming uniformelastic constants of the host material λ GaN while solid lines are for spatially varying elastic constants with correction order λ ( N ) .As the correction order λ ( N ) increases, the strain moves toward the uniform case λ GaN . Note that the lines for λ (2) and λ (3) areoverlapping. However, the uniform case λ GaN underestimates the strain in the dot. the strain is calculated in its own real space unit cell Ω s with volume V s . We denote the strain reciprocal unit cellas Ω − s such that it contains the wave vectors Q , which are defined similarly to Eq. 1 for the electronic cell. Giventhat strain relaxes more slowly than bound electronic wavefunctions, we only consider V s ≥ V e . In this two-unit-cellsapproximation, the Fourier transform of the strain produced by the quantum dot array is ˜ (cid:15) a ij ( q ) = 1 V s (cid:88) Q ∈ Ω − s ˜ (cid:15) ij ( Q ) ˜ χ e ( q − Q )= 1 V s (˜ (cid:15) ij ∗ ˜ χ e ) s ( q ) , (11)where χ e is the characteristic function of the electronic unit cell Ω e in Ω s , which is given for our case in AppendixB. Superscript “ a ” indicates array. We follow the notation that q ∈ Ω − e and Q ∈ Ω − s . (˜ (cid:15) ij ∗ ˜ χ e ) s ( q ) denotes aconvolution where the subscript “ s ” indicates that the convolution is over the wave vectors Q ∈ Ω − s , see AppendixC for Fourier transform and convolution definitions. We show in Sec. III C that choosing the linear dimensions of Ω s to be integer multiples of the linear dimensions of Ω e ensures that all vectors q ∈ Ω − e are also in Ω − s . This choiceallows Eq. 11 to be evaluated efficiently. C. Piezoelectric potential
III-nitride materials are strongly piezoelectric, having both spontaneous and strain-driven polarizations [16, 17].Calculation of the polarization from an electric field requires knowledge of the static dielectric constant ε of thematerial. In prior work, all Fourier-space based approaches assumed a uniform dielectric constant. We present amethod to obtain the Fourier transform of the scalar potential ˜ ϕ ( q ) assuming ε ( r ) changes with the local alloyfraction. We find that correcting for the spatial dependence of the dielectric function leads to important changes inthe piezoelectric potential. We also show in Sec. V that this change in piezoelectric potential significantly shifts thelowest quantum dot energy levels. We do not discuss metallic screening, which can be important in highly dopedmaterials [18–20].Generally, we can write the displacement field D ( r ) as D ( r ) = ε E ( r ) + P tot ( r ) , where E ( r ) is the electric field, ε is the vacuum permittivity, and P tot is the total polarization. In the strainedmaterial, there are three sources of polarization: bound charge, strain and spontaneous polarization, P tot ( r ) = P bnd ( r ) + P st ( r ) + P sp ( r ) . Here, we assume no free charge screening and so an intrinsic material. Assuming P bnd to be linear with the electricfield and incorporated into ε ( r ) as usual, D ( r ) = ε ( r ) E ( r ) + P st ( r ) + P sp ( r ) , (12)where P st ( r ) + P sp ( r ) = P ( r ) is the residual polarization after electric-field induced bound charge has been includedin ε ( r ) .We take ε ( r ) to be ε h in the host material and ε d in the dot material, so ε ( r ) = ε h + (cid:0) ε d − ε h (cid:1) χ d ( r ) . (13)We obtain ε d by linear interpolation of the binary compounds’ bulk dielectric constants. Taking the divergence of Eq.12, using ∇ · D = 0 , taking the Fourier transform and solving for the electric field gives E m ( r ) = − ε ( r ) F − (cid:26) q n q m ˜ P n ( q ) (cid:27) (14)where F − represents the inverse Fourier transform. Using E m = − ∂ m ϕ , where ∂ m ≡ ∂∂x m and ϕ ( r ) is the scalarpotential, ˜ ϕ ( q ) = − iq m F (cid:26) ε ( r ) F − (cid:26) q n q m ˜ P n ( q ) (cid:27) ( r ) (cid:27) ( q ) . (15)For the case of sharp alloy interfaces, χ d ( r ) is either 1 or 0 and Eq. 13 gives ε ( r ) = 1 ε h + (cid:18) ε d − ε h (cid:19) χ d ( r ) . (16)We treat the case of smoothly varying alloy fraction in Sec. IV. Putting this result in Eq. 15 gives ˜ ϕ ( q ) = ˜ ϕ huni ( q ) + ∆ ˜ ϕ ( q ) (17)with ˜ ϕ huni ( q ) = − iq m ε h q n q m ˜ P n ( q ) (18) ∆ ˜ ϕ ( q ) = − iq m (cid:18) ε d − ε h (cid:19) F (cid:26) χ d ( r ) F − (cid:26) q n q m ˜ P n ( q ) (cid:27)(cid:27) (19)Here, ˜ ϕ huni is the contribution to ϕ with ε r ( r ) = ε h , and ∆ ˜ ϕ is the change in ˜ ϕ due to the dot material having adifferent dielectric constant. The polarization fields ˜ P n ( q ) for the wurtzite crystal structure are given in terms ofstrain in Appendix E.We now show the piezoelectric potentials that result from this formulation, for our model system described in Sec.II A. Figure 4 shows ϕ ( z ) along the central axis of the quantum dot calculated with constant ε of the dot and hostand with Eq. 17. The calculation with spatially varying ε ( r ) agrees with ϕ duni inside the dot and also agrees with ϕ huni outside the dot, with a transition near the boundary that is captured by neither of the uniform cases.We showed in Fig. 3 how spatially varying elastic constants change strain profiles. Figure 5 shows how ϕ changesdue to the elastic constants’ correction propagates into the piezoelectric potential. We find that the changes inpiezoelectric potential, a peak correction of 8 mV, are significant if looking to converge the energy levels within a fewmeV’s. FIG. 4. Piezoelectric potential ϕ ( z ) along the central axis, beginning in the center of the dot, with parameters as in Table I.Blue and red dashed lines show ϕ calculated with uniform ε ( r ) = ε h and ε ( r ) = ε d , respectively. Black line shows the case withspatially varying ε ( r ) . Note that the ϕ is antisymmetric in z . These results show that the case of uniform ε cannot capture ϕ throughout the system. Strain calculations include third order corrections for the nonuniform elastic constants.FIG. 5. Piezoelectric potential difference along the central axis of the dot for various elastic constant corrections. Potentialdifference is with respect to ϕ h ( r ) , which is calculated with a uniform λ ( r ) = λ GaN . Calculations implement an alloy smoothingof δ = [1 . , . , .
5] ˚ A , which is described in Sec. IV, and the vertical dashed line indicates the nominal material interfacewithout smoothing. These results are for the same system as Fig. 4, calculated using various orders of correction for the effectsof nonuniform elastic constants, described in Sec. II B 1. Spatially varying dielectric constants are included. The zeroth ordercase (blue line) shows λ ( r ) = λ d . III. SYMMETRY ADAPTED BASIS k · p FOR WURTZITE QUANTUM DOTS
Here, we present the quantum dot k · p model we use for electronic structure calculations. We first present a theoryfor bulk materials and use it to construct a theory for quantum dots. This quantum dot Hamiltonian is written in asymmetry adapted basis, which reduces the computational cost for calculating and diagonalizing the Hamiltonian. Inthis symmetry adapted basis, we show how the strain produced by the quantum dots contributes to the Hamiltonian.We also introduce strain effects using a different unit cell than that of the electronic cell defined in Fig. 1. In thissection, our goal is to show our method of efficiently including strain in the quantum dot k · p model, which we doby choosing the strain unit cell’s dimensions to be integer multiples of the unit cell used for the electronic structurecalculations. A. Bulk k · p model To describe the electronic structure of bulk wurtzite systems, we use an 8-band k · p model, which includes spin-orbitcoupling, crystal field splitting and strain. An 8-band k · p model for bulk wurtzite material has been presented byRef. [7] in the basis of Γ -point Bloch functions. References [7, 21] presented a 6-band model using eigenfunctions ofthe angular momentum operator ˆ J z . Since choosing ˆ J z eigenfunctions aids in the construction of a symmetry adaptedbasis, which is presented in Sec. III B, we have used these two references to construct an 8-band Hamiltonian in the ˆ J z eigenfunctions basis. More precisely, we have constructed the Hamiltonian using S , X , Y and Z Γ -point Blochfunctions as a basis and then performed a basis transformation to obtain the ˆ J z eigenfunctions basis. While k · p parameters are usually obtained in S , X , Y and Z basis, recent work has obtained k · p parameters directly in thesymmetry adapted basis using ab initio calculations [22].We consider the time-independent Schrödinger equation for a single electron ˆ H | ψ (cid:105) = E | ψ (cid:105) (20)where ˆ H = ˆ K + ˆ V + ˆ H so + ˆ H cr + ˆ H st (21)Here, ˆ K is the kinetic term of the electrons, ˆ V the potential from the electron-ion interaction, ˆ H so is spin-orbitcoupling, ˆ H cr is crystal field splitting and ˆ H st is strain coupling. We expand | ψ (cid:105) in terms of the ˆ J z eigenfunctions | u i (cid:105)| ψ (cid:105) = e i k · r (cid:88) α =1 C α | u i (cid:105) (22)where | u (cid:105) = | iS, ↑(cid:105) | u (cid:105) = |− iS, ↓(cid:105)| u (cid:105) = |− X + iY √ , ↑(cid:105) | u (cid:105) = | X − iY √ , ↓(cid:105)| u (cid:105) = | X − iY √ , ↑(cid:105) | u (cid:105) = |− X + iY √ , ↓(cid:105)| u (cid:105) = | Z, ↑(cid:105) | u (cid:105) = | Z, ↓(cid:105) (23)Here, S , X , Y and Z are Γ -point Bloch functions with arrows indicating spin. The eigenvalues of the ˆ J z eigenfunctionsare J z = (cid:26) , , − , , − , − , , − (cid:27) , respectively. Inserting Eq. 22 into 20, the eigenvalue problem can be written as H α (cid:48) α C α = EC α (cid:48) (24)Keeping terms only up to order k , the 8x8 k · p Hamiltonian is H = (cid:34) g ( k ) γ − γ ∗ g ∗ ( k ) (cid:35) (25)where g ( k ) = g ( k ) + g ( k ) + g cr + g so + g st g ( k ) = E (cid:48) c − P √ k + P √ k − P k z − P √ k − E (cid:48) v P √ k + E (cid:48) v P k z E (cid:48) v g ( k ) = A (cid:48) (cid:0) k x + k y (cid:1) + A (cid:48) k z (cid:16) L (cid:48) + M (cid:17) (cid:0) k x + k y (cid:1) + M k z − N (cid:48) k − − √ N (cid:48) k − k z − N (cid:48) k (cid:16) L (cid:48) + M (cid:17) (cid:0) k x + k y (cid:1) + M k z √ N (cid:48) k + k z − √ N (cid:48) k + k z √ N (cid:48) k − k z M (cid:0) k x + k y (cid:1) + L (cid:48) k z g st = a ( (cid:15) xx + (cid:15) yy ) + a (cid:15) zz ( l + m ) ( (cid:15) xx + (cid:15) yy ) + m (cid:15) zz − ( l − m ) ( (cid:15) xx − (cid:15) yy ) + in (cid:15) xy − n ( (cid:15) xz − i(cid:15) yz ) √ − ( l − m ) ( (cid:15) xx − (cid:15) yy ) − in (cid:15) xy ( l + m ) ( (cid:15) xx + (cid:15) yy ) + m (cid:15) zz n ( (cid:15) xz + i(cid:15) yz ) √ − n ( (cid:15) xz + i(cid:15) yz ) √ n ( (cid:15) xz − i(cid:15) yz ) √ m ( (cid:15) xx + (cid:15) yy ) + l (cid:15) zz (26) g cr = ∆ cr g so = ∆ so − γ = √ so − Here, ∆ cr and ∆ so are the crystal field splitting and spin-orbit coupling, respectively. The band edges are E (cid:48) c = E v + E g + ∆ cr + ∆ so + ϕ and E (cid:48) v = E v + ϕ where ϕ is any additional scalar potential such as the piezoelectric potential.The A (cid:48) i parameters are related to the Kane parameters P i and L (cid:48) i , M i , N (cid:48) i to the Luttinger-like parameters A i , all ofwhich are shown in Appendix A. g st is the contribution to the Hamiltonian due to strain (cid:15) ij . The parameters a i , l i . m i and n i for the strain contribution are given in Appendix A in terms of deformation potentials.In example calculations, alloy parameters have been obtained by linearly interpolating between bulk GaN and InNparameters, which are given in Appendix A, except the band gap, which has bowing included. B. Quantum dot k · p For the quantum dot system, we construct the Hamiltonian from the bulk system described in Section III A. We useslowly varying envelope functions and apply a spatial dependence to the bulk Hamiltonian. The problem is expressedin a symmetry adapted basis to obtain a block diagonal Hamiltonian from which we calculate the eigenstates of thequantum dot.We start from Eq. 20, but expand | ψ (cid:105) in terms of envelope functions F α ( r ) that are slowly varying compared tothe lattice constant [8, 9, 11, 23], | ψ (cid:105) = (cid:88) α =1 | F, α (cid:105)(cid:104) r | F, α (cid:105) = F α ( r ) u α ( r ) (27)1where the u α ( r ) are defined by Eq. 23 and are periodic with the crystal lattice. Analogous to Eq. 24, this envelopefunction expansion leads to (cid:88) α =1 H α (cid:48) α F β ( r ) = EF α (cid:48) ( r ) (28)where H α (cid:48) α are the bulk Hamiltonian matrix elements from Eq. 25. Due to the broken translation symmetry in thequantum dot system, we apply the substitution k j → − i ∂∂x j (29)to the bulk Hamiltonian in Eq. 25. Each parameter in the bulk Hamiltonian also has a spatial dependence based onthe alloy distribution, f ( r ) = f d χ d ( r ) + f h [1 − χ d ( r )] (30)Here, f stands for any of the parameters in the bulk Hamiltonian that are material dependent. f h and f d are theparameter values of the host and alloyed dot material, respectively. Applying the substitution in Eq. 29 to Eq. 25, theHamiltonian consists of terms of the form f ( r ) , f ( r ) ∂∂x j and f ( r ) ∂ ∂x i ∂x j . To preserve Hermiticity, we symmetrizethe derivatives [9, 23, 24]: f ( r ) ∂∂x j → (cid:18) f ( r ) ∂∂x j + ∂∂x j f ( r ) (cid:19) (31) f ( r ) ∂ ∂x i ∂x j → (cid:18) ∂∂x i f ( r ) ∂∂x j + ∂∂x j f ( r ) ∂∂x i (cid:19) (32)The envelope functions F α ( r ) are periodic with the superlattice and can be expanded in Fourier domain using thesuperlattice reciprocal wave vectors q defined in Eq. 1. Writing the envelope functions F α ( r ) in terms of plane wavesleads to a non-sparse Hamiltonian [8, 23]. For computational efficiency, we use a symmetry adapted basis, which takesadvantage of the C symmetry of the wurtzite crystal structure by block diagonalizing the Hamiltonian. Symmetryadapted bases have been fully described for both zincblende and wurtzite systems [9, 10]. We use a symmetry-adaptedbasis with elements | m f , α, q (cid:105) where q = ( q x , q y , q z ) are chosen within a single sextant, so ≤ q y ≤ tan (cid:0) π (cid:1) q x , and m f = {− / , − / , − / , / , / , / } can be interpreted as a total quasi angular momentum [10, 11]. This basisconsists of the basis functions of the irreducible representations of the double group ¯ C . Using this basis reduces theFourier space sampling to a single sextant of the full space and block diagonalizes the Hamiltonian into 6 blocks, whichare labeled by m f . This block diagonalization greatly reduces the computational cost to diagonalize the Hamiltonian.Figure 6 shows an example of the Fourier space sampling used in the symmetry adapted basis. Written out, the basisstates are | m f , α, q (cid:105) = Λ ( m f , α, q , r ) | u α (cid:105) (33) Λ ( m f , α, q , r ) = √ (cid:80) l =0 e i (cid:16) ←→ R l q (cid:17) · r e il π [ m f − J z ( α )] q x (cid:54) = 0 or q y (cid:54) = 0 e i q · r q x = q y = 0 J z ( α ) = m f (34)where ←→ R l is the l π rotation around the z-axis. Equation 34 distinguishes wave vectors that are purely along thez-axis from those that have an xy-component, which we denote by q z and q , respectively. These two cases differbecause a z-axis rotation leaves q z invariant while sending q to a new wave vector ←→ R l q . The case of | m f , α, q z (cid:105) with J z ( α ) (cid:54) = m f does not exist in the basis set. Using the symmetry adapted basis, the eigenstates can be written | ψ i,m f (cid:105) = (cid:88) α =1 (cid:88) q A αim f ( q ) | m f , α, q (cid:105) (35)where the q summation is restricted to the sextant, shown in Fig. 6.2Writing the envelope functions in the symmetry adapted basis, the eigenvalue problem in Eq. 28 can then be written (cid:88) α =1 (cid:88) q H m f α (cid:48) α ( q (cid:48) , q ) A αim f ( q ) = E i A α (cid:48) im f ( q (cid:48) ) (36)with H m f α (cid:48) α ( q (cid:48) , q ) ≡ (cid:104) m f , α (cid:48) , q (cid:48) | ˆ H| m f , α, q (cid:105) = 1 V e (cid:90) V e d r Λ ∗ ( r ) H α (cid:48) α Λ ( r ) (37)where H α (cid:48) α are the bulk Hamiltonian matrix elements presented in Sec. III A. Expressions for H m f α (cid:48) α ( q (cid:48) , q ) are fullywritten out in Appendix F in terms of the bulk Hamiltonian matrix elements and quantum dot characteristic function. FIG. 6. Fourier space sampling used for the symmetry adapted basis. Red circles are the Fourier space points q used in thesymmetry adapted basis. Blue dots are the full Fourier space sampled by rotations ←→ R l q . Dashed lines highlight the 6 sextants. C. Including strain and piezoelectric effects
Deformation potentials and piezoelectric effects, which are both strain-driven, are important for accurate calcula-tions of electronic structure in III-N materials. However, including deformation potentials can be computationallycostly for the case of isolated dots. The two-unit cell approach presented in Sec. II B 2 allows for the study of isolateddots, but at the cost of computationally expensive convolutions. Additionally, another layer of convolutions appearsin the Hamiltonian matrix elements, leading to composed convolutions. Here, we present the matrix elements due tostrain and show our computationally efficient approach of dealing with these composed convolutions by choosing thelinear dimensions of the real-space strain cell Ω s to be integer multiples of those of the electronic cell Ω e .The bulk strain Hamiltonian matrix elements in Eqs. 25 and 26 can be written as H α (cid:48) α = (cid:88) ij f ijα (cid:48) α (cid:15) ij ( r ) where the f ijα (cid:48) α consist of k · p parameters ( a i , l i , m i and n i ). Using the prescription of Sec. III B, the straincontributions to the quantum dot Hamiltonian are3 H ij, st m f α (cid:48) α ( q (cid:48) , q ) = 16 (cid:88) l (cid:48) =0 5 (cid:88) l =0 e i π { l [ m f − J z ( α )] − l (cid:48) [ m f − J z ( α (cid:48) )] } h ij, st α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , ←→ R l q (cid:17) H ij, st m f α (cid:48) α ( q (cid:48) , q z ) = 1 √ (cid:88) l (cid:48) =0 e − il (cid:48) φ [ m f − J z ( α (cid:48) )] h ij, st α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , q z (cid:17) H ij, st m f α (cid:48) α ( q (cid:48) z , q z ) = h ij, st α (cid:48) α ( q (cid:48) z , q z ) where h ij, st α (cid:48) α ( q (cid:48) , q ) = (2 π ) f ij, h α (cid:48) α V e ˜ (cid:15) a ij ( q (cid:48) − q ) + (2 π ) (cid:16) f ij, d α (cid:48) α − f ij, h α (cid:48) α (cid:17) V e (cid:0) ˜ χ d ∗ ˜ (cid:15) a ij (cid:1) e ( q (cid:48) − q ) . (38)Here, f ij, h α (cid:48) α and f ij, d α (cid:48) α are the k · p parameters for bulk host and dot materials, respectively. ˜ (cid:15) a ij is the strain producedby the quantum dot array calculated in Sec. II B. The subscript “e” in (cid:0) ˜ χ d ∗ ˜ (cid:15) a ij (cid:1) e indicates that the convolution is overthe wave vectors q ∈ Ω − e . Inserting the superlattice strain from Eq. 11 into the k · p strain matrix elements from Eq.38 leads to composed convolutions, (cid:0) ˜ χ d ∗ ˜ (cid:15) a ij (cid:1) e ( q ) = (cid:88) q (cid:48) ∈ Ω − e ˜ χ d ( q (cid:48) ) ˜ (cid:15) a ij ( q − q (cid:48) ) (39) = 1 V s (cid:88) q (cid:48) ∈ Ω − e ˜ χ d ( q (cid:48) ) (cid:88) Q ∈ Ω − s ˜ (cid:15) ij ( Q ) ˜ χ e ( q − q (cid:48) − Q ) (40)which can be computationally demanding depending on the number of wave vectors used. The original proposal ofusing a large strain cell with a smaller electronic cell imposed no relationship between their sizes [11]. Equation 40then requires evaluating ˜ χ e at points q − q (cid:48) − Q , which are contained on neither the electronic nor strain meshes,requiring a unique convolution be calculated for every q (cid:48) . It is well known that using the convolution theorem tocompute a convolution between two vectors of length N has a computational cost that scales as N log ( N ) . Similarly,the computational cost for a convolution on a 3D N × N × N mesh scales as N log ( N ) . Computing the composedconvolutions in Eq. 40 would then scale as N e log ( N e ) N s log ( N s ) since a convolution in Q has to be calculated foreach individual q (cid:48) . Note that the convolutions from Eqs. 39-40 are linear convolutions, which implies that the arraysof function values must be padded with zeros before using the convolution theorem as detailed in Appendix C. Thiszero padding increases both N e and N s . We show that choosing a strain unit cell to be a supercell of the electronicunit cell reduces the number of convolutions to compute, leading to an improved scaling of N e log ( N e ) + N s log ( N s ) .Choosing the strain unit cell linear dimensions to be multiples of the electronic cell, we have L s i = n i L e i i = 12 , (41)where the n i take positive integer values. This choice of real-space unit cells leads to the electronic Fourier-spacemesh being contained in the strain mesh Ω − e ⊂ Ω − s . The wave vectors Q then have a spacing that is a fraction ofthe spacing of the electronic wave vectors q , ∆ Q i = ∆ q i n i (42)Note that from Eq. 39, ˜ (cid:15) a ij is only sampled at points ∆ q = q − q (cid:48) , which belong to the electronic mesh. Ourprocedure starts with using the convolution theorem (see Appendix C) to efficiently calculate the inner convolution (˜ (cid:15) ij ∗ ˜ χ e ) s ( Q ) on the strain mesh to obtain ˜ (cid:15) a ij ( Q ) . Since the wave vectors Q also contain the wave vectors q ,we can then extract the points that lie on the electronic mesh to obtain ˜ (cid:15) a ij ( q ) . Lastly, we perform the secondconvolution (cid:0) ˜ χ d ∗ ˜ (cid:15) a ij (cid:1) e ( q ) , again utilizing the convolution theorem. This workflow is shown in Fig. 8(a). In ourmethod, we compute only two 3D convolutions and so get a complexity scaling of N e log ( N e ) + N s log ( N s ) , which isa considerable improvement compared to the non-overlapping case. Note that N s is generally much larger than N e toobtain appropriate convergence, so the computational cost is dominated by the convolutions on Ω − s .The piezoelectric potential brings no additional complexity, and the workflow for calculating the piezoelectricpotential is shown in Fig. 8(b). The potential is initially calculated on the strain mesh, and the electronic meshportion is extracted to calculate the piezoelectric potential contributions to the Hamiltonian, which are written outin Appendix F.4 FIG. 7. Example of Fourier space meshes where n = 2 , leading to ∆ Q i = ∆ q i , which implies the vectors Q contain all thevectors q . (a)Extractto e-mesh ˜ (cid:15) a ij ( Q ) (cid:124) (cid:123)(cid:122) (cid:125) −→ ˜ (cid:15) a ij ( q ) −→ (cid:0) ˜ χ d ∗ ˜ (cid:15) a ij (cid:1) e ( q ) (cid:124) (cid:123)(cid:122) (cid:125) Array strain Calculate on e-meshon s-mesh (b) ˜ (cid:15) arr ( Q ) −→ ˜ ϕ ( Q ) (cid:124) (cid:123)(cid:122) (cid:125) −→ ˜ ϕ arr ( Q ) = ( ˜ ϕ ∗ ˜ χ e ) s ( Q ) (cid:124) (cid:123)(cid:122) (cid:125) −→ ˜ ϕ arr ( q ) (cid:124) (cid:123)(cid:122) (cid:125) Strain gives potential Truncate to Extractin strain box electronic unit cell on e-meshFIG. 8. (a) Workflow for calculation of composed convolutions of the form shown in Eq. 39. (b) Workflow used to obtainpiezoelectric potential on the electronic mesh from the strain on the strain mesh. e-mesh and s-mesh signify the Fourier spaceelectronic and strain meshes, respectively.
IV. SMOOTH ALLOY PROFILE
When InGaN devices are grown by molecular beam epitaxy (MBE), indium diffuses between layers [1]. While moststudies of MBE-grown materials simulate abrupt junctions, this diffusion leads to smoothing of the material interfaces,producing a continuously varying alloy fraction, which changes the local band properties and lattice constant, which inturn change strain and polarization fields. This smooth alloy profile must be included for accurate modeling. Smoothindium profiles also provide a computational benefit, since sharp features of the confining potentials are removed, sofewer wave vectors are required to attain convergence. In this section, we present a method to include alloy diffusioneffects by effectively smoothing the characteristic function of the dot. We focus on indium alloying here for theexamples, but the methods are general for all k · p calculations of alloy structures.
1. Smoothing method
In the case of a sharp material interface, the local alloy fraction X ( r ) can be defined by the characteristic functionof the dot X ( r ) = X χ d ( r ) , where the characteristic function χ d ( r ) defines the geometry of the dot with indium fraction X . By convolving witha Gaussian G ( r , δ ) = √ πδ x δ y δ z e − (cid:18) x δ x + y δ y + z δ z (cid:19) or other kernel, we can obtain a smooth version of the characteristic5 FIG. 9. (a) Band gap in a quantum dot system obtained by bowed interpolation using Eq. 46 along the z and y axes throughthe center of the dot. Band gap obtained from linear interpolation using Eq. 47 is visually indistinguishable. (b) Differenceof the bowed and linearly interpolated band gaps along the z and x axes. Quantum dot parameters are listed in Table I.However, for computational simplicity and smooth curves, these results were obtained using a rectangular real-space unit cellwith dimensions L x = L y = 500 ˚ A and L z = 70 ˚ A and a smoothing of δ = [3 , ,
5] ˚ A . function X sm ( r ) = ( X χ d ∗ G ) ( r )= X χ sm ( r ) , where δ = [ δ x , δ y , δ z ] controls the radius of smoothing and needs to be chosen to model the desired alloy diffusion. G ( r , δ ) is normalized to preserve the total amount of alloying element, and χ sm ( r ) is a smoothed characteristicfunction. Using the convolution theorem, the smoothed characteristic function satisfies ˜ χ sm ( q ) = ˜ χ d ( q ) e − ( δ xq x + δ yq y + δ zq z ) . Note that χ sm ( r ) is no longer strictly a characteristic function, as it takes values continuously between 0 and 1. Wenow show that it can be inserted in place of the characteristic function in the previous sections to give k · p parameters,strain and piezoelectric fields accurately with a smooth alloy profile.
2. Material parameters
We now focus on the case of InGaN to illustrate the interpolation of material parameters. In the case of sharpmaterial interfaces, the host and dot regions each consist of uniform material. The host material is a binary materialand has well-defined parameters. The dot region consists of alloyed InGaN, and its parameters are obtained by eitherlinear or bowed interpolation of bulk GaN and InN parameters, which are listed in Appendix A.In the case of a smooth alloy profile, the dot and host regions are no longer uniform, giving the material parameters asmooth spatial dependence. Parameters that were linearly interpolated in the sharp interface case can still be obtainedfrom a simple linear interpolation based on the local alloy fraction X ( r ) . The band gap E g is nonlinear in the alloyfraction due to a bowing factor. This nonlinearity prevents us from using the convolution theorem in calculating the6 FIG. 10. (a) Inverse dielectric constant from Eq. 48 along the z and x axes through the center of the dot. Linearly interpolatedinverse from Eq. 49 is visually indistinguishable. (b) Relative difference of the linear and nonlinear dielectric constants alongthe z and x axes. Quantum dot parameters are listed in Table I. However, for computational simplicity and smooth curves,these results were obtained using a rectangular real-space unit cell with dimensions L x = L y = 500 ˚ A and L z = 70 ˚ A and asmoothing of δ = [3 , ,
5] ˚ A . Hamiltonian matrix elements. However, we show that neglecting the bowing parameters in the alloy-smoothing regioncan still give computationally efficient and accurate smoothed profiles when the alloy fraction is not too large.The local value for any of the linearly interpolated material parameters depends on the local alloy fraction f ( r ) = f B X ( r ) + [1 − X ( r , X )] f A (43)where f can be a parameter such as lattice constant, and subscripts A and B stand for the two binary materials, GaNand InN for example. For this case of linearly interpolated quantities, smoothed parameters can be written: f ( r ) = f d ( X ) χ sm ( r ) + [1 − χ sm ( r )] f A (44)where f d ( X ) is the linearly interpolated material parameter at the nominal alloy fraction X of the quantum dot.Band gaps do not vary linearly with alloy fraction and are generally well described with a bowing term, as E g ( r ) = E B g X ( r ) + [1 − X ( r )] E A g − X ( r ) [1 − X ( r , X )] C (45)where C is a bowing constant. Following the same procedure as in Eq. 44, a smoothed version can be written: E g ( r ) = X χ sm ( r ) E Bg + [1 − X χ sm ( r )] E Ag − CX χ sm ( r ) [1 − X χ sm ( r )]= E Ag + (cid:2) E Bg − E Ag (cid:3) E Ag X χ sm ( r ) − CX χ sm ( r ) [1 − X χ sm ( r )] (46)where the first two terms are the linear interpolation and the last term is the bowing. This bowing term bringsadditional complexity when performing k · p calculations due to the nonlinearlity in χ sm ( r ) . We approximate theband gap by a linear interpolation between the host and dot band gaps, E g ( r ) ≈ E d g ( X ) χ sm ( r ) + [1 − χ sm ( r )] E h g . (47)Here, E d g is the bulk band gap at an alloy fraction of X and E h g is the bulk band gap of the host material. Thislinear interpolation gives a good approximation for the band gap for most regions and as well for moderate indiumfractions, as shown in Fig. 9. The regions with largest deviation are in the same locations where E g changes over 1.5eV, so we expect the slight shift of position where each band gap value occurs to have minimal effect. The neglect ofthe χ sm term allows the theory to stay linear and therefore efficiently calculated with the convolution theorem.7
3. Strain and the piezoelectric potential
Here we show how smoothing is included in the strain and piezoelectric potential calculations. Once calculated,those strains and piezoelectric potentials can be included in the k · p model exactly as shown in Sec. III C.Following the derivations from Refs. [8, 13], it is not obvious how smoothing is to be implemented in straincalculations since they begin from the stress of the sharp interface dot/barrier interface. However, Ref. [15] presents analternative derivation for the same strain calculation indicating that χ ( r ) in the strain expressions can be exchangedfor the smoothed version χ sm ( r ) without any further changes.For the piezoelectric potential, Eq. 16 for the spatially varying inverse dielectric constant assumed sharp interfaces.In the case of a smooth indium profile, we use Eq. 43 to write ε ( r ) = ε B X ( r ) + [1 − X ( r )] ε A (48)In the scenario where X ( r ) is spatially varying, Eq. 16 can no longer be applied, because the inverse of the dielectricconstant is not a linear function of indium. However, similar to the band gap, we find that ε ( r ) ≈ ε d ( X ) χ sm ( r ) + [1 − χ sm ( r )] 1 ε A (49)still gives an accurate representation of ε − ( r ) . Figure 10 shows a disagreement of less than 1 percent between theinverse dielectric from Eq. 48 and the linear interpolation in Eq. 49. The form of Eq. 49 allows us to use Eqs. 17-19for the piezoelectric potential with a simple substitution of χ ( r ) by χ sm ( r ) . FIG. 11. Convergence of the fundamental gap of the quantum dot E for a 1D array of quantum dots with the maximummagnitude of q included in the k · p calculations. L e and n vary while (a) has R = 200 ˚ A with strain box size held constantat L s = 10 R and (b) has R = 40 ˚ A with L s = 16 R . Plane wave sampling m from 3 to 11 are shown for each choice of L e ; the number next to each point indicates m . Different values of m , L e that produce the same q max can be seen toproduce approximately the same E , showing that q max is a useful metric for convergence of these states. Since q max = m π/L e and computational cost scales with m , smaller L e allows easier access to large q max . In both panels, the black curves have L e = 2 R , so the dots touch each other at the 6 edges of the hexagonal unit cell. For larger dot, there is no visible deviationof E from the trend with separated dots. With the smaller dot, tunneling of the wavefunctions into neighboring dots causes E to have a significant change, labeled ∆ . FIG. 12. Electronic structure of the QD superlattice system along the central axis of the dot for δ = [1 . , . , .
5] ˚ A , corre-sponding to s = 1 in Fig. 13(a). Other system parameters are in Table I. Lowest bound electron and hole states energies areshown by the horizontal dashed lines. Thick solid lines are the bulk band edges under the influence of the piezoelectric fieldand strain. Thin solid lines are z-axis projections of the probability distributions obtained from the envelope functions. Dashedvertical lines are the nominal material interfaces before smoothing. These calculations include spatially varying elastic anddielectric constants.TABLE II. Energy shifts due to spatially varying λ ( r ) and ε ( r ) relative to the case with uniform constants of the host material λ GaN and ε GaN , respectively. System parameters in Table I.Energy shifts λ ( r ) & ε GaN λ GaN & ε ( r ) λ ( r ) & ε ( r )∆ E c (meV) 16.7 46.7 64.7 ∆ E v (meV) -5.1 -30.6 -37.4 ∆ E (meV) 21.7 77.4 102.1 V. IMPACTS OF CORRECTIONS
In this section, we apply our methodology to study the case of a 1D array of quantum dots, such as described inRef. [1], though we do not consider the boundaries of the nanowire. We achieve this 1D array by taking n = 1 tofully couple the dots in the z-direction and n > to avoid strain effects from neighboring dots in the xy-plane. Inthis section, we investigate convergence of the lowest electron and hole state energies E c and E v , which define thefundamental gap of the dot E = E c − E v . More specifically, we show that the largest wave vector sampled plays adominant role in convergence. We also show the energy shifts experienced by these two states when using uniform orspatially varying material parameters and when including alloy smoothing.We model an infinite 1D quantum dot array with parameters listed in Table I. The 1D dot array has an experi-mentally well characterized dot-dot spacing in z, which fixes L e = L s = L , leaving L and n to be fixed. Thesequantum dots have a rather large radius, so the smallest spatial feature that we need to resolve is the decay of thebound wavefunctions into the classically forbidden region. Given that bound wavefunctions decay faster than strain,we need wave vectors that are relatively large to be able to resolve the wavefunctions. Increasing m increases themaximum wave vector contained in the mesh, but we can also sample at larger wave vectors by using a smaller L e .However, if the electronic cell is chosen too small, then there can be electronic wavefunction overlap between statesof neighboring dots. We must then choose L e as small as possible while also avoiding dot-dot interactions. As forstrain, in order to study a 1D array, we must choose n sufficiently large to have L s = n L e large enough that thestrain of the quantum dot superlattice does not extend across neighboring strain unit cells in the xy-plane.With this intuition, we turn to the convergence of E in terms of m , L e and n . Figure 11(a) shows theimportance of the largest q in the electronic mesh, q max , for convergence of E . In this study, L e and n are chosento keep a constant L s = L e n = 2000 ˚ A . We observe that E is to good approximation a function of only q max andnot of m and L e , converging towards the same value for all choices of L e . We also observe that the smallest L e with highest m gives the most converged E , since q max = m π/L e . The black line in Fig. 11(a) represents thecase of dots touching in the xy plane and, interestingly, does not break the convergence trend. However, we do finda break in the convergence trend for smaller dots in Fig. 11(b). This difference in convergence is due to the largerquantum dots having better confined states compared to the smaller dots. Smaller dots have wavefunctions thatextend further outside the dot region, which makes them more able to tunnel to a neighboring dot. Consequently,9 FIG. 13. Effects of indium diffusion, given by δ = s [1 . , . , .
5] ˚ A , for (a) E c and E v and (b) bulk conduction band edge.Remaining parameters are as listed in Table I. Increase in indium diffusion leads to less confinement, which pushes the twostates apart in energy, widening the electronic gap E . Electronic structure for s = 1 is shown in Fig. 12. care has to be taken in choosing the unit cell dimensions for small quantum dots.The lowest quantum dot confined electron and hole energies, E c and E v , have respectively been converged within meV by choosing m , m and n sufficiently large, see Table I. Material parameters are listed in Appendix A. Bandedges and lowest-energy confined states are shown in Fig. 12. The thick black solid lines represent the bulk bandedges modified by the piezoelectric potential and strain. To include strain effects in the bulk band edges, we haveused the (1 , matrix element from Eq. 26 to modify the conduction band edge and a third of the trace of the × valence band block for the valence band edge.The modifications in both the strain and piezoelectric potential due to spatially varying elastic and dielectricconstants also have effects on the electronic structure. Table II shows how much E v and E c shift due to the corrections.We find that both corrections push the states apart, leading to an energy gap 100 meV larger than from simplercalculations with uniform ε and λ , a significant change that shows the importance of accurate modeling of dielectricand elastic parameters.Figure 13(a) shows that indium diffusion pushes the lowest electron and hole states apart, which is due to changesin the confining potentials. From Fig. 13(b), we see that indium diffusion reduces the depth of the confining potential.We have observed similar behavior for the hole state confining potential, leading to E v being pushed down in energy.Consequently, the gap E increases in energy as indium diffusion is increased.In the case of sharp material interfaces, large wave vectors are needed to resolve the discontinuous parameterprofiles. Smoothing removes the sharp interfaces and yields smoothly varying material parameters. Consequently,the required q max for the same degree of convergence is smaller, which means that smaller m and m and thereforereduced computational cost are needed with increasing δ . In Figs. 13 (a) and (b), we converged for the case of s = 1 ,guaranteeing convergence for the rest of the sweep. VI. CONCLUSIONS
We have demonstrated techniques for and results of four modifications of standard quantum dot k · p theory. Wehave included spatially varying elastic and dielectric constants as alloy fraction changes in strain and piezoelectricpotential calculations. The effects of the spatially varying parameters are non-negligible on the strain and piezoelectricpotential and also produce important shifts of the lowest electron and hole states, significantly changing the calculatedgap of the quantum dot. We have also presented a method to include smoothly varying alloy profiles in Fourier-based0strain, piezoelectric potential and k · p calculations. This smoothing has to be chosen to represent the device of interest,such as for indium diffusion in InGaN systems. For the case of k · p theory for isolated dots, we have presented a newmethodology of overlapping electronic and strain meshes to facilitate the coupling of strain into the k · p Hamiltonian,greatly reducing the computational cost of calculating the Hamiltonian matrix elements. Lastly, we have shown thatthe maximum wave vector contained in the electronic sampling mesh is the most important criterion for determiningconvergence of quantum dot levels.
ACKNOWLEDGMENTS
We acknowledge useful conversation with Stanko Tomić about Fourier-space k · p methods. We acknowledge fundingfrom the Ontario Early Researcher Award and NSERC CREATE TOP-SET program, Award number 497981. Appendix A: Bulk k · p parameters From Ref. [7], the k · p parameters A (cid:48) i are related to the Kane parameters P i as A (cid:48) = (cid:126) m (cid:107) e − P E g A (cid:48) = (cid:126) m ⊥ e − P E g where P = (cid:126) m (cid:18) m m (cid:107) e − (cid:19) E g (∆ so + E g ) + ∆ cr (2∆ so + 3 E g )2∆ so + 3 E g P = (cid:126) m (cid:18) m m ⊥ e − (cid:19) E g [3 E g (∆ so + E g )] + ∆ cr (2∆ so + 3 E g )∆ cr ∆ so + 3∆ cr E g + 2∆ so E g + E g . Here, m (cid:107) e and m ⊥ e are the electron effective masses along the z-axis and in the xy-plane, respectively. ∆ cr and ∆ so are the crystal field splitting and spin-orbit coupling, respectively. The Luttinger-like parameters L (cid:48) i , M i and N (cid:48) i arerelated to the A i parameters by L (cid:48) = (cid:126) m ( A + A + A ) + P E g L (cid:48) = (cid:126) m A + P E g M = (cid:126) m ( A + A − A ) M = (cid:126) m ( A + A ) M = (cid:126) m A N (cid:48) = (cid:126) m A + P E g N (cid:48) = (cid:126) m √ A + P P E g Note there is an error in the relations for L (cid:48) , L (cid:48) and N (cid:48) in Ref. [7], which we have corrected in agreement withAppendix B of Ref. [25]. The parameters A i , P i and E g used in our numerical study of InGaN systems are given inTable III.1Similarly to the k · p parameters L (cid:48) i , M i and N (cid:48) i , the strain parameters are l = D + D + D l = D m = D + D − D m = D + D m = D + D m = D n = 2 D n = √ D where the deformation potentials D i are listed in Table III. TABLE III. Material parameters used to model the InGaN system. Parameters were taken from Ref. [7].Parameters GaN InN a ( ˚A ) 3.189 3.545 c ( ˚A ) 5.185 5.703 C (GPa) 390 223 C (GPa) 145 115 C (GPa) 106 92 C (GPa) 398 224 C (GPa) 105 48 e ( C / m ) 0.326 0.264 e ( C / m ) -0.527 -0.484 e ( C / m ) 0.895 1.06 P sp ( C / m ) -0.034 -0.042 ε r E g (eV) 3.51 0.78 E v (eV) 0 0.5 ∆ cr (eV) 0.010 0.040 ∆ so (eV) 0.017 0.005 m (cid:107) /m m ⊥ /m A -7.21 -8.21 A -0.44 -0.68 A A -3.46 -5.23 A -3.40 -5.11 A -4.90 -5.96 a (eV) -4.9 -3.5 a (eV) -11.3 -3.5 D (eV) -3.7 -3.7 D (eV) 4.5 4.5 D (eV) 8.2 8.2 D (eV) -4.1 -4.1 D (eV) -4.0 -4.0 D (eV) -5.5 -5.5 Appendix B: Characteristic functions
The characteristic function of a single dot is unity inside the dot and zero outside, χ d ( r ) = (cid:40) r ∈ Ω d otherwisewhere Ω d is the space inside the dot. For a cylindrical dot centered on the origin with radius R and height h alongthe z-axis, the Fourier transform of χ d is ˜ χ d ( q ) = 1(2 π ) πRq (cid:113) q x + q y sin (cid:18) h q (cid:19) J (cid:16) R (cid:113) q x + q y (cid:17) . The characteristic function of the electronic cell, the hexagonal prism shown in Fig. 1, is defined as χ e ( r ) = (cid:40) r ∈ Ω e otherwiseand its Fourier transform is ˜ χ e ( q ) = L sinc (cid:18) q L (cid:19) q cos (cid:16) L q (cid:17) + q cos (cid:16) L q (cid:17) − ( q + q ) cos (cid:0) L q + q (cid:1) √ q q ( q + q ) . Appendix C: Fourier and Convolution Conventions
We define the Fourier forward and inverse transforms of a function g ( r ) as F { g ( r ) } ( q ) ≡ ˜ g ( q ) = 1(2 π ) ∞ (cid:90) −∞ d r g ( r ) e − i q · r F − { ˜ g ( q ) } ( r ) = g ( r ) = ∞ (cid:90) −∞ d q ˜ g ( q ) e i q · r A convolution is denoted by ( f ∗ g ) ( r ) = ∞ (cid:90) −∞ f ( r (cid:48) ) g ( r − r (cid:48) ) d r (cid:48) The convolution theorem states (cid:16) ˜ f ∗ ˜ g (cid:17) ( q ) = F { f ( r ) g ( r ) } ( q ) . For a system that is periodic in real space with a unit cell Ω of volume V , we define ˜ g ( q ) = 1 V (cid:90) Ω d r g ( r ) e − i q · r g ( r ) = (cid:88) q ∈ Ω − ˜ g ( q ) e i q · r where Ω − is the reciprocal space to the unit cell Ω . Defining the Fourier space convolution as (cid:16) ˜ f ∗ ˜ g (cid:17) ( q ) = (cid:88) q (cid:48) ∈ Ω − ˜ f ( q (cid:48) ) ˜ g ( q − q (cid:48) ) , (C1)3the convolution theorem is (cid:16) ˜ f ∗ ˜ g (cid:17) ( q ) = F { f ( r ) g ( r ) } ( q ) . Since we consider different superlattice unit cells for the electronic and strain properties, (cid:16) ˜ f ∗ ˜ g (cid:17) e indicates a convo-lution on the electronic space Ω − e and (cid:16) ˜ f ∗ ˜ g (cid:17) s on the strain space Ω − s .The reciprocal space Ω − contains a discrete infinity of wave vectors on which to evaluate ˜ f and ˜ g . The convolutionin Eq. C1 then sums over the infinite number of wave vectors, with ˜ f and ˜ g being functions that decay at large wavevectors. For the calculations in this manuscript, we choose a finite number of wave vectors. By choosing this mesh tocontain wave vectors sufficiently large to capture the decay of ˜ f and ˜ g , we can calculate the linear convolution in Eq.C1 to good approximation by padding the ˜ f and ˜ g arrays with zeros and performing a circular convolution, definedbelow.We denote the finite Fourier-space mesh by q i i i with − m j < i j < m j such that N j = 2 m j + 1 is the dimension ofthe mesh in each direction. For simplicity, we use the mapping p j = i j + m j + 1 to start indexing from 1. Evaluating afunction ˜ f on the mesh q p p p gives the array ˜ f p p p . The Fourier-space array ˜ f p p p and its real-space counterpart f v v v are then related through the discrete Fourier transform and its inverse, ˜ f p p p = F { f } p p p = 1 N N ,N ,N (cid:88) v ,v ,v =1 f v v v e − i π ( p v /N + p v /N + p v /N ) (C2) f v v v = F − (cid:110) ˜ f (cid:111) v v v = N ,N ,N (cid:88) p ,p ,p =1 ˜ f p p p e i π ( p v /N + p v /N + p v /N ) . (C3)where N = N N N . The circular convolution is defined as (cid:16) ˜ f • ˜ g (cid:17) p p p = N ,N ,N (cid:88) p (cid:48) ,p (cid:48) ,p (cid:48) =1 ˜ f p (cid:48) p (cid:48) p (cid:48) ˜ g ( p − p (cid:48) ) , ( p − p (cid:48) ) , ( p − p (cid:48) ) where ˜ f ( p + N ) , ( p + N ) , ( p + N ) = ˜ f p p p and likewise for ˜ g . The convolution theorem is then (cid:16) ˜ f • ˜ g (cid:17) p p p = F { f g } p p p . To perform a linear convolution, we pad the arrays with zeros, which increases the dimensions of the mesh to N p j = N j + 2 m j and yields the padded array ˜ f p p p p containing N p = N p N p N p elements. The linear convolution isthen (cid:16) ˜ f ∗ ˜ g (cid:17) p p p = (cid:16) ˜ f p • ˜ g p (cid:17) p p p (C4)This result is independent of the basis used to generate the mesh and is valid in the case of hexagonal meshes. Appendix D: Displacement field Green’s tensor
The Green’s function for the displacement field for spatially varying elastic constants must satisfy [8] ∂∂x k λ iklm ( r ) ∂∂x m G ln ( r − r (cid:48) ) = − δ ( r − r (cid:48) ) δ in From Ref. [8], the Green’s tensor when λ jklm is spatially invariant is ˜ G h in = 1(2 π ) βq + ρq (cid:26) δ in − (cid:2)(cid:0) ρq + γq (cid:1) δ i + ( κ + ρ ) q i q (cid:3) F δ n − IF P − IQ q n − [( α + β ) q i + ( κ + ρ ) δ i q ] P − Qδ n F P − IQ q n (cid:27) , F ( q ) = ( C + 2 C − C ) q + C q Q ( q ) = ( C − C − C + C ) q + ( C + 2 C − C ) q P ( q ) = C q + ( C − C − C ) q I ( q ) = ( C + C ) q α = C β = 12 ( C − C ) γ = C − C − C + C κ = C − C ρ = C + C − C with C ij being the elastic constants. Appendix E: Polarization fields
There are two contributions to the piezoelectric polarization fields: strain driven and spontaneous. The strain-drivenpolarization can be written [10] P st P st P st = e (cid:15) e (cid:15) e ( (cid:15) + (cid:15) ) + e (cid:15) , where e ij are the piezoelectric constants and (cid:15) ij the strain fields. We write the piezoelectric constants using thecharacteristic function of the quantum dot e ij ( r ) = e d ij χ d ( r ) + e h ij [1 − χ d ( r )] . The Fourier transform of the polarization fields is then ˜ P st ˜ P st ˜ P st = e h ˜ (cid:15) ( q ) + π ) ( e d − e h ) V ( ˜ χ d ∗ ˜ (cid:15) ) e ( q )2 e h ˜ (cid:15) ( q ) + π ) ( e d − e h ) V ( ˜ χ d ∗ ˜ (cid:15) ) e ( q ) e h [˜ (cid:15) ( q ) + ˜ (cid:15) ( q )] + (2 π ) ( e d − e h ) V [ ˜ χ d ∗ (˜ (cid:15) + ˜ (cid:15) )] e ( q ) + e h ˜ (cid:15) ( q ) + (2 π ) ( e d − e h ) V ( ˜ χ d ∗ ˜ (cid:15) ) e ( q ) In a bulk wurtzite material, the spontaneous polarization is along the c-axis and uniform throughout the material P sp ( r ) = P sp Then ˜ P sp ( q ) = P sp , h δ q , + (cid:0) P sp , d − P sp , h (cid:1) ˜ χ ( q ) where P sp , d and P sp , h are the spontaneous polarizations in the dot and host materials, respectively. Appendix F: QD k · p Hamiltonian
For the quantum dot superlattice system, the k · p Hamiltonian matrix elements of Eq. 37 in the symmetry adaptedbasis are given in terms of the parameters f d α (cid:48) α and f h α (cid:48) α that make up the bulk k · p Hamiltonian matrix elements H α (cid:48) α , presented in Sec. III. We take the convention where superscript “ (0) ” indicates a bulk Hamiltonian matrix5element containing no wave vector, “ ( i ) ” indicates a single wave vector k i and “ ( i, j ) ” two wave vectors k i and k j . Bydefining φ = π and S ll (cid:48) m f α (cid:48) ,α = e iφ { l [ m f − J z ( α )] − l (cid:48) [ m f − J z ( α (cid:48) )] } S l (cid:48) m f α (cid:48) = e − il (cid:48) φ [ m f − J z ( α (cid:48) )] , the quantum dot Hamiltonian matrix elements are H (0) m f α (cid:48) α ( q (cid:48) , q ) = 16 (cid:88) l (cid:48) =0 5 (cid:88) l =0 S ll (cid:48) m f α (cid:48) ,α h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , ←→ R l q (cid:17) H ( i ) m f α (cid:48) α ( q (cid:48) , q ) = 16 (cid:88) l (cid:48) =0 5 (cid:88) l =0 S ll (cid:48) m f α (cid:48) ,α (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) i + (cid:16) ←→ R l q (cid:17) i h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , ←→ R l q (cid:17) H ( i,j ) m f α (cid:48) α ( q (cid:48) , q ) = 16 (cid:88) l (cid:48) =0 5 (cid:88) l =0 S ll (cid:48) m f α (cid:48) ,α (cid:16) ←→ R l q (cid:17) j (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) i + (cid:16) ←→ R l q (cid:17) i (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) j h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , ←→ R l q (cid:17) H (0) m f α (cid:48) α ( q (cid:48) , q z ) = 1 √ (cid:88) l (cid:48) =0 S l (cid:48) m f α (cid:48) h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , q z (cid:17) H ( i ) m f α (cid:48) α ( q (cid:48) , q z ) = 1 √ (cid:88) l (cid:48) =0 S l (cid:48) m f α (cid:48) (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) i + ( q z ) i h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , q z (cid:17) H ( i,j ) m f α (cid:48) α ( q (cid:48) , q z ) = 1 √ (cid:88) l (cid:48) =0 S l (cid:48) m f α (cid:48) ( q z ) j (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) i + ( q z ) i (cid:16) ←→ R l (cid:48) q (cid:48) (cid:17) j h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , q z (cid:17) H (0) m f α (cid:48) α ( q (cid:48) z , q z ) = h α (cid:48) α ( q (cid:48) z , q z ) H ( i ) m f α (cid:48) α ( q (cid:48) z , q z ) = ( q (cid:48) z ) i + ( q z ) i h α (cid:48) α ( q (cid:48) z , q z ) H ( i,j ) m f α (cid:48) α ( q (cid:48) z , q z ) = ( q z ) j ( q (cid:48) z ) i + ( q z ) i (0 , , q (cid:48) z ) j h α (cid:48) α ( q (cid:48) z , q z ) where h α (cid:48) α (cid:16) ←→ R l (cid:48) q (cid:48) , ←→ R l q (cid:17) = f h α (cid:48) α δ l,l (cid:48) δ q,q (cid:48) + (2 π ) (cid:0) f d α (cid:48) α − f h α (cid:48) α (cid:1) V ˜ χ d (cid:16) ←→ R l (cid:48) q (cid:48) − ←→ R l q (cid:17) The contributions of the piezoelectric potential to the Hamiltonian are H pz m f α (cid:48) α ( q (cid:48) , q ) = − δ α (cid:48) ,α (2 π ) e c V (cid:88) l (cid:48) =0 5 (cid:88) l =0 e iφ ( l − l (cid:48) ) [ m f − J z ( α )] ˜ ϕ (cid:16) ←→ R l (cid:48) q (cid:48) − ←→ R l q (cid:17) H pz m f α (cid:48) α ( q (cid:48) , q z ) = − δ α (cid:48) ,α (2 π ) e c V √ (cid:88) l (cid:48) =0 e − il (cid:48) φ [ m f − J z ( α (cid:48) )] ˜ ϕ (cid:16) ←→ R l (cid:48) q (cid:48) − q z (cid:17) H pz m f α (cid:48) α ( q (cid:48) z , q z ) = − δ α (cid:48) ,α (2 π ) e c V ˜ ϕ ( q (cid:48) z − q z ) ˜ ϕ is described in Secs. II C and III C and e c is the electric charge. [1] H. P. T. Nguyen, S. Zhang, K. Cui, X. Han, S. Fathololoumi, M. Couillard, G. A. Botton, and Z. Mi, “p-Type ModulationDoped InGaN/GaN Dot-in-a-Wire White-Light-Emitting Diodes Monolithically Grown on Si(111),” Nano Lett. , 1919–1924 (2011).[2] Tim J. Puchtler, Tong Wang, Christopher X. Ren, Fengzai Tang, Rachel A. Oliver, Robert A. Taylor, and Tongtong Zhu,“Ultrafast, Polarized, Single-Photon Emission from m-Plane InGaN Quantum Dots on GaN Nanowires,” Nano Lett. ,7779–7785 (2016).[3] Md G. Kibria, Hieu P.T. Nguyen, Kai Cui, Songrui Zhao, Dongping Liu, Hong Guo, Michel L. Trudeau, Suzanne Paradis,Abou Rachid Hakima, and Zetian Mi, “One-step overall water splitting under visible light using multiband InGaN/GaNnanowire heterostructures,” ACS Nano , 7886–7893 (2013).[4] Liwen Sang, Meiyong Liao, Qifeng Liang, Masaki Takeguchi, Benjamin Dierre, Bo Shen, Takashi Sekiguchi, YasuoKoide, and Masatomo Sumiya, “A Multilevel Intermediate-Band Solar Cell by InGaN/GaN Quantum Dots with a Strain-Modulated Structure,” Adv. Mater. , 1414–1420 (2014).[5] Ross Cheriton, Sharif M. Sadaf, Luc Robichaud, Jacob J. Krich, Zetian Mi, and Karin Hinzer, “Two-photon photocurrentin InGaN/GaN nanowire intermediate band solar cells,” Communications Materials , 63 (2020).[6] T. Saito and Y. Arakawa, “Electronic structure of piezoelectric In . Ga . N quantum dots in GaN calculated using atight-binding method,” Physica E , 169–181 (2002).[7] Momme Winkelnkemper, Andrei Schliwa, and Dieter Bimberg, “Interrelation of structural and electronic properties inIn x Ga − x N/GaN quantum dots using an eight-band k · p model,” Phys. Rev. B , 155322 (2006).[8] A. D. Andreev and E. P. O’Reilly, “Theory of the electronic structure of GaN/AlN hexagonal quantum dots,” Phys. Rev.B , 15851–15870 (2000).[9] Nenad Vukmirović, Dragan Indjin, Vladimir D. Jovanović, Zoran Ikonić, and Paul Harrison, “Symmetry of k.p Hamiltonianin pyramidal InAs/GaAs quantum dots: Application to the calculation of electronic structure,” Phys. Rev. B , 075356(2005).[10] Nenad Vukmirović, Zoran Ikonić, Dragan Indjin, and Paul Harrison, “Symmetry-based calculation of single-particle statesand intraband absorption in hexagonal GaN/AlN quantum dot superlattices,” J. Phys. Condens. Matter , 6249–6262(2006).[11] Nenad Vukmirovć and Stanko Tomić, “Plane wave methodology for single quantum dot electronic structure calculations,”J. Appl. Phys. , 103718 (2008).[12] O. Stier, M. Grundmann, and D. Bimberg, “Electronic and optical properties of strained quantum dots modeled by 8-bandk.p theory,” Phys. Rev. B , 5688–5701 (1999).[13] A. D. Andreev, J. R. Downes, D. A. Faux, and E. P. O’Reilly, “Strain distributions in quantum dots of arbitrary shape,”J. Appl. Phys. , 297–305 (1999).[14] J. Renard, R. Songmuang, G. Tourbot, C. Bougerol, B. Daudin, and B. Gayral, “Evidence for quantum-confined starkeffect in gan/aln quantum dots in nanowires,” Phys. Rev. B , 121305 (2009).[15] A. V. Nenashev, A. A. Koshkarev, and A. V. Dvurechenskii, “Approximate analytical description of the elastic strain fielddue to an inclusion in a continuous medium with cubic anisotropy,” Journal of Applied Physics , 105104 (2018).[16] Fabio Bernardini, Vincenzo Fiorentini, and David Vanderbilt, “Spontaneous polarization and piezoelectric constants ofIII-V nitrides,” Phys. Rev. B , R10024–R10027 (1997).[17] Agostino Zoroddu, Fabio Bernardini, Paolo Ruggerone, and Vincenzo Fiorentini, “First-principles prediction of struc-ture, energetics, formation enthalpy, elastic constants, polarization, and piezoelectric constants of AlN, GaN, and InN:Comparison of local and gradient-corrected density-functional theory,” Phys. Rev. B , 045208 (2001).[18] S. F. Chichibu, A. C. Abare, M. S. Minsky, S. Keller, S. B. Fleischer, J. E. Bowers, E. Hu, U. K. Mishra, L. A. Coldren,S. P. DenBaars, and T. Sota, “Effective band gap inhomogeneity and piezoelectric field in ingan/gan multiquantum wellstructures,” Applied Physics Letters , 2006–2008 (1998).[19] J. P. Ibbetson, P. T. Fini, K. D. Ness, S. P. DenBaars, J. S. Speck, and U. K. Mishra, “Polarization effects, surface states,and the source of electrons in algan/gan heterostructure field effect transistors,” Applied Physics Letters , 250–252(2000).[20] Hwa-mok Kim, Yong-hoon Cho, Hosang Lee, Suk Il Kim, Sung Ryong Ryu, Deuk Young Kim, Tae Won Kang, andKwan Soo Chung, “High-Brightness Light Emitting Diodes Using Dislocation-Free Indium Gallium Nitride/Gallium NitrideMultiquantum-Well Nanorod Arrays,” Nano Letters , 1059–1062 (2004).[21] S. L. Chuang and C. S. Chang, “k · p method for strained wurtzite semiconductors,” Phys. Rev. B , 2491–2504 (1996).[22] Milan Jocić Jocić and Nenad Vukmirović, “Ab initio construction of symmetry-adapted k.p hamiltonians for the electronicstructure of semiconductors,” Phys. Rev. B , 085121 (2020).[23] Stanko Tomić, Andrew G. Sunderland, and Ian J. Bush, “Parallel multi-band k·p code for electronic structure of zincblend semiconductor quantum dots,” J. Mater. Chem. , 1963–1972 (2006).[24] Richard A. Morrow and Kenneth R. Brownstein, “Model effective-mass Hamiltonians for abrupt heterojunctions and theassociated wave-function-matching conditions,” Phys. Rev. B , 678–680 (1984).[25] E. Berkowicz, D. Gershoni, G. Bahir, E. Lakin, D. Shilo, E. Zolotoyabko, A. C. Abare, S. P. Denbaars, and L. A. Coldren, “Measured and calculated radiative lifetime and optical absorption of in x ga − x N / GaN quantum structures,” Phys. Rev. B61