k·p Method for Strained Wurtzite GaN Band Structure

In the tree, select .

The geometry can be a simple square for setting up the Floquet-Bloch boundary conditions. The size of the domain should be smaller than the shortest wavelength, to prevent spurious modes (a related mechanical model is thin_film_baw_resonator_dispersion_diagram in the MEMS Module, where a narrow domain is used).

Geometry 1

The Model Wizard ended at the pane for the node in the Model Builder tree structure. We can use this opportunity to select a convenient length unit.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Square 1 (sq1)

In the toolbar, click

In the window for , click

In this 2D model we will look at the valence band dispersion along the x and z directions. Change the default second axis name from y to z for clarity.

Component 1 (comp1)

In the window, click .

In the window for , locate the section.

Find the subsection. In the table, enter the following settings:


First	Second	Third
x	z	y

Find the subsection. In the table, enter the following settings:


First	Second	Third
X	Z	Y

Enter the GaN material parameters as listed in Table III of Chuang and Chang’s paper. Only the valence band parameters are needed for this model. First enter the energy parameters Δ1, Δso, Δ2, Δ3, and Δ.

Global Definitions

Parameters 1 - GaN

In the window, under click .

In the window for , type Parameters 1 - GaN in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
Del1	16[meV]	2.5635E-21 J	Energy parameters
Delso	12[meV]	1.9226E-21 J
Del2	Delso/3	6.4087E-22 J
Del3	4[meV]	6.4087E-22 J
Del	sqrt(2)*Del3	9.0633E-22 J

Next enter the valence band effective-mass parameters A1 ~ A6.


Name	Expression	Value	Description
A1	-6.56	-6.56	Valence band effective mass parameters
A2	-0.91	-0.91
A3	5.65	5.65
A4	-2.83	-2.83
A5	-3.13	-3.13
A6	-4.86	-4.86

Next enter the deformation potentials D1 ~ D4.


Name	Expression	Value	Description
D1	0.7[eV]	1.1215E-19 J	Deformation potentials
D2	2.1[eV]	3.3646E-19 J
D3	1.4[eV]	2.243E-19 J
D4	-0.7[eV]	-1.1215E-19 J

Next enter the elastic stiffness constants C13 and C33.

Finally to prevent too long a list of parameters when setting up sweeps, exclude these parameters from the selection list.


Name	Expression	Value	Description
C13	15.8e11[dyn/cm^2]	1.58E11 N/m²	Elastic stiffness constants
C33	26.7e11[dyn/cm^2]	2.67E11 N/m²

Click to expand the section. Clear the check box.

Now set up the swept parameters.

Parameters 2 - Sweeps

In the toolbar, click

and choose .

In the window for , type Parameters 2 - Sweeps in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
epsxx	-0.01	-0.01	Strain,xx
epsyy	epsxx	-0.01	Strain,yy
epszz	-2C13/C33epsxx	0.011835	Strain,zz

To conveniently create a plot to compare with Fig. 5 in the paper, set up a swept parameter kp such that its positive axis represents the kx axis and the negative axis represents the kz axis.

Parameters 3 - Analytic formula


Name	Expression	Value	Description
kp	0[rad/nm]	0 rad/m	Positive axis: kx, negative axis: kz
kx	if(kp>0,kp,0[rad/nm])	0 rad/m
ky	0[rad/nm]	0 rad/m
kz	if(kp<0,-kp,0[rad/nm])	0 rad/m

For comparing with analytic solutions, enter the formulas from Eq. (34) and Eq. (42) in the paper. These parameters also are not needed to be in the selection list. In these formulas, the wave vectors kx and kz are numerical parameters. Therefore a 3X3 Hermitian matrix can be constructed and diagonalized to provide the eigenvalues for analytic comparison. Later on when setting up the Schrödinger equation, the wave vectors kx and kz will be replaced by differential operators, not just numbers. This distinction should be kept in mind. First enter the diagonal parameters F, G, λ, λε, θ, and θε in Eq. (34).

In the toolbar, click

and choose .

In the window for , type Parameters 3 - Analytic formula in the text field.

Locate the section. Clear the check box.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
F	Del1+Del2+lm+th	2.7003E-21 J	Eq. (34)
G	Del1-Del2+lm+th	1.4186E-21 J
lm	hbar_const^2/2/me_const(A1kz^2+A2*(kx^2+ky^2))+lmeps	-5.4018E-21 J
lmeps	D1epszz+D2(epsxx+epsyy)	-5.4018E-21 J
th	hbar_const^2/2/me_const(A3kz^2+A4*(kx^2+ky^2))+theps	4.8977E-21 J
theps	D3epszz+D4(epsxx+epsyy)	4.8977E-21 J

Then enter the off-diagonal terms Kt and Ht in Eq. (42).

Global ODEs and DAEs (ge)


Name	Expression	Value	Description
Kt	hbar_const^2/2/me_constA5kt^2	-0 J	Eq. (42)
Ht	hbar_const^2/2/me_constA6kt*kz	-0 J
kt	sqrt(kx^2+ky^2)	0 rad/m

Set up a global equation to diagonalize the 3X3 Hermitian matrix - the upper-left Hamiltonian in Eq. (45) in the paper. Use the reserved name lambda for the eigenvalue in the global equation. Scale the Hamiltonian with the energy scale of 1 meV so that the source term is about order of unity, and the eigenenergy is the eigenvalue in the unit of meV. Add some blank spaces in the expressions to align the columns of the matrix.

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Global Equations 1

In the window, under click .

In the window for , locate the section.

Schrödinger Equation (schr)


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (u_t0) (1/s)	Description
u1	(F u1+Kt u2-iHt u3)/1[meV]-lambda*u1	0	0
u2	(Kt u1+G u2+(Del-iHt)u3)/1[meV]-lambda*u2	0	0
u3	(iHtu1+(Del+iHt)u2+lm u3)/1[meV]-lambdau3	0	0

The text is yellow colored because at this point the eigenvalue variable lambda is not defined yet.

Set up the physics. Unlike the analytic formulas where everything is a number, here in the Schrödinger equation the wave vectors kx and kz will be replaced by differential operators. Thus some terms contain just numbers while others contain operators, and they will be entered in different physics features. Some patience and attention to details will help a lot to avoid mistakes. (To compare with Fig. 5 in the paper and for simplicity, let ky stay at zero.)

For example, consider the (1,1) element of the matrix F =Δ 1 +Δ 2 +λ +θ. It contains terms that are just numbers: Δ1 +Δ 2 +λ ε + θε. It also has terms containing the wave vectors kx and kz, which will be replaced by the differential operators id/dx and id/dz, respectively (the differentiations here and next are partial derivatives). Note that there is no minus sign in front of the imaginary unit i because of the engineering sign convention adopted by all COMSOL physics interfaces: a plane wave is exp(-ikx + iωt), not exp(ikx − iωt). Thus the terms for F that contains differential operators are (hbar2/2me)(id/dz)(A1 + A3)(id/dz) and (hbar2/2me)(id/dx)(A2 + A4)(id/dx).

Since these two terms have second order differentiations, they are entered using the feature. Each term occupies one row in the . The position (1,1) in the matrix and the two differentiation operators are specified by drop-down menus. The A parameters are entered in the input field. A factor of (hbar2/2me) is already included in the feature.

To conveniently create a plot to compare with Fig. 5 in the paper, set the eigenvalue scale to the same energy unit as the vertical axis (meV), so that the eigenvalue takes on the numerical value of the eigenenergy in the same unit (meV).

In the window, under click .

In the window for , locate the section.

Find the subsection. In the λscale text field, type 1[meV].

Second Order Hamiltonian 1: Diagonal F, G, lambda

In the toolbar, click

and choose .

In the window for , type Second Order Hamiltonian 1: Diagonal F, G, lambda in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
1	1	i ?/?z	A1+A3	i ?/?z	F: lambda+theta

Note in the table how one differential operator is located on the left side of the A parameter and the other one is on the right. In general, such as in heterostructures or quantum dots, the A parameter can be spatially varying. In those cases, the differential operator on the left side acts on both the A parameter and the wave function. This is shown in the Equation section as well. Here for bulk material the A parameter is just a constant. Now add a row for the other term (hbar2/2me)(id/dx)(A2 + A4)(id/dx).

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
1	1	i ?/?x	A2+A4	i ?/?x	F: lambda+theta

This finishes the two terms with differential operators for the (1,1) element of the matrix F. Before entering the terms that are just numbers for F, while we are still in this feature, enter the terms with differential operators for the (2,2) and (3,3) elements G and λ.

The following instructions add the contributions from G and λ one row at a time. Alternatively, it is possible to save some time by copying and pasting. For example, since the contribution to the second order Hamiltonian from F is the same as the one from G, except the position in the matrix is different - F is at (1,1) and G is at (2,2), we can select the two rows for F in the table by mouse drag, right-click to copy, click on the button to add a new row, and then right click the first cell of the new row to paste the two rows that have been copied. After that, change the position of the two pasted rows from (1,1) to (2,2) and update the description from F to G, leaving the differential operators and the A parameters the same as the ones for F.

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
2	2	i ?/?z	A1+A3	i ?/?z	G: lambda+theta

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
2	2	i ?/?x	A2+A4	i ?/?x	G: lambda+theta

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
3	3	i ?/?z	A1	i ?/?z	lambda

Electron Potential Energy 1: Diagonal F, G, lambda

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
3	3	i ?/?x	A2	i ?/?x	lambda

The feature that we just finished entering for the diagonal elements corresponds to the kinetic energy terms normally added by the default feature. So we should disable it to remove the unwanted default contribution to the Hamiltonian.

Effective Mass 1

In the window, right-click and choose .

Now enter the terms that are just numbers for the 3 diagonal elements F, G and λ. They are Δ1 +Δ 2 +λ ε + θε for F, Δ1 −Δ 2 +λ ε + θε for G, and simply λε for λ. For the diagonal elements, it is easiest to use the default feature.

In the window, click .

In the window for , type Electron Potential Energy 1: Diagonal F, G, lambda in the text field.

Locate the section. In the Ve,11 text field, type Del1+Del2+lmeps+theps.

In the Ve,22 text field, type Del1-Del2+lmeps+theps.

In the Ve,33 text field, type lmeps.

We have finished entering the 3 diagonal elements F, G and λ. Now enter the off-diagonal elements. Again some are just numbers (Δ), and others have partial derivatives (Kt and Ht). They will be entered into different features. First the ones with partial derivatives.

Second Order Hamiltonian 2: Off-diagonal Kt, Ht

In the toolbar, click

and choose .

In the window for , type Second Order Hamiltonian 2: Off-diagonal Kt, Ht in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
1	2	i ?/?x	A5	i ?/?x	Kt

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
2	1	i ?/?x	A5	i ?/?x	Kt

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
1	3	i ?/?x	-i*A6	i ?/?z	Ht

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
3	1	i ?/?x	i*A6	i ?/?z	Ht

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
2	3	i ?/?x	-i*A6	i ?/?z	Ht

Zeroth Order Hamiltonian 1: Off-diagonal Delta

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Left differential operator (i)	A parameter (1)	Right differential operator (j)	Description
3	2	i ?/?x	i*A6	i ?/?z	Ht

Finally enter the off-diagonal elements that are just numbers (Δ), using the feature. Remember to compensate for the factor of (hbar2/2me) that is inherent in all of this type of Hamiltonian features.

In the toolbar, click

and choose .

In the window for , type Zeroth Order Hamiltonian 1: Off-diagonal Delta in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Zeroth-order Hamiltonian parameter (1/m^2)	Description
2	3	Del2me_const/hbar_const^2	Delta

Select Boundaries 1 and 4 only.

In the table, enter the following settings:


Hamiltonian row index (m)	Hamiltonian column index (n)	Zeroth-order Hamiltonian parameter (1/m^2)	Description
3	2	Del2me_const/hbar_const^2	Delta

Next set up the Floquet-Bloch periodic boundary conditions. Use the parameters kx and kz for the wave vector. Select opposing boundary pairs.

Periodic Condition 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Specify the kF vector as


kx	X
kz	Z

Periodic Condition 2

Right-click and choose .

In the window for , locate the section.

Select Boundaries 2 and 3 only.

Now set up the mesh. Since we have created a domain that is much smaller than the wavelength, there is no need for many mesh elements. In addition, the periodic boundary conditions require the meshes on the pair of source and destination boundaries to match. In this case, a simple Mapped Mesh works well.

Mesh 1

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 2.

In the toolbar, click

Set up studies to compute the valence band structure for an unstrained and a 1% compressively strained GaN wurtzite crystal, to be compared to Fig. 5 in the paper. Sweep the parameter epsxx for the strain. Sweep the parameter kp prepared earlier such that on its negative axis, kx=0 and kz=-kp, while on its positive axis, kz=0 and kx=kp. Use one study for the Schrödinger equation and a second study for the global equation (analytic formula).

Study 1: Fig. 5. Schrödinger Eq.

In the window, click .

In the window for , type Study 1: Fig. 5. Schrödinger Eq. in the text field.

Locate the section. Clear the check box.

In the window, under click .

In the window for , locate the section.

In the text field, type 10.

Locate the section. In the table, clear the check box for .

Click to expand the section. Select the check box.

From the list, choose .


Parameter name	Parameter value list	Parameter unit
epsxx (Strain,xx)	0 -0.01

Create plots to compare to Fig. 5 in the paper. First the unstrained case. Make one global plot for each band by selecting the appropriate eigenvalue from the drop-down menu for the dataset.


Parameter name	Parameter value list	Parameter unit
kp (Positive axis: kx, negative axis: kz)	range(-0.12,0.005,0.12)	1/angstrom

In the toolbar, click

Global 1 - Schrödinger Eq. HH

Fig. 5. Unstrained

In the toolbar, click

and choose .

In the window for , type Fig. 5. Unstrained in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

Locate the section. Select the check box.

In the associated text field, type Negative axis: kz Positive axis: kx (1/angstrom).

Select the check box.

In the associated text field, type Energy (meV).

Locate the section. Select the check box.

In the text field, type -0.12.

In the text field, type 0.12.

In the text field, type -100.

In the text field, type 50.

Right-click and choose .

In the window for , type Global 1 - Schrödinger Eq. HH in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. Click

Global 1 - Schrödinger Eq. LH


Expression	Unit	Description
lambda		HH

Locate the section. From the list, choose .

From the list, choose .

In the text field, type kp/1[angstrom^-1].

Right-click and choose .

In the window for , type Global 1 - Schrödinger Eq. LH in the text field.

Locate the section. From the list, choose .

In the text field, type 2.

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
lambda		LH

Global 1 - Schrödinger Eq. CH

Right-click and choose .

In the window for , type Global 1 - Schrödinger Eq. CH in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
lambda		CH

Duplicate the plot for the strained case. Change the plot label and the selection of the strain parameter epsxx.

Fig. 5. -1% Strained

In the window, right-click and choose .

In the window for , type Fig. 5. -1% Strained in the text field.

Global 1 - Schrödinger Eq. HH

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Global 1 - Schrödinger Eq. LH

In the window, click .

In the window for , locate the section.

From the list, choose .

Global 1 - Schrödinger Eq. CH

In the window, click .

In the window for , locate the section.

From the list, choose .

Now set up a study to solve the analytic 3X3 matrix equation configured with the global equation. Use the All eigenvalue option to obtain all 3 eigenvalues of the 3X3 matrix equation.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 1: Fig. 5. Schrödinger Eq.

Study 2: Fig. 5. Analytic

In the window, under right-click and choose .

In the window, click .

In the window for , type Study 2: Fig. 5. Analytic in the text field.

Locate the section. Clear the check box.

Right-click and choose .

Add the analytic solution to the 2 plots to compare with the numerical solution.

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Physics interface	Solve for	Equation form
Schrödinger Equation (schr)		Automatic (Eigenvalue)
Global ODEs and DAEs (ge)	√	Not applicable

Locate the section. In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
kp (Positive axis: kx, negative axis: kz)	range(-0.12,0.01,0.12)	1/angstrom

In the toolbar, click

Add a new study to compute the 2D band structure in the kx-kz plane, first using the Schrödinger equation.

Global 1 - Analytic

In the window, right-click and choose .

In the window for , type Global 1 - Analytic in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
lambda		Analytic

Click to expand the section. Find the subsection. From the list, choose .

Find the subsection. From the list, choose .

In the toolbar, click

Right-click and choose .

Global 1 - Analytic

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 3: 2D Dispersion Schrödinger Eq.

In the window for , type Study 3: 2D Dispersion Schrödinger Eq. in the text field.

Locate the section. Clear the check box.

Study 1: Fig. 5. Schrödinger Eq.

Study 3: 2D Dispersion Schrödinger Eq.

In the window, under right-click and choose .

In the window, right-click and choose .

In the window for , locate the section.


Parameter name	Parameter value list	Parameter unit
kx	range(0,0.1,1)*0.07	1/angstrom
kz	range(0,0.1,1)*0.07	1/angstrom

In the toolbar, click

Create one table for each band to plot the 2D band structure using the plot type for the 3 hole bands. First the heavy hole band (HH).

Global Evaluation 1: Schrödinger eq.

Evaluation Group 1: HH

In the toolbar, click

In the window for , type Evaluation Group 1: HH in the text field.

Locate the section. From the list, choose .

From the list, choose .

Right-click and choose .

In the window for , type Global Evaluation 1: Schrödinger eq. in the text field.

Locate the section. Click


Expression	Unit	Description
lambda		HH (Schr)

In the toolbar, click

Table

Go to the window.

Click in the window toolbar.

Table Surface 1: HH (Schr)

In the window, under click .

In the window for , type Table Surface 1: HH (Schr) in the text field.

Locate the section. From the list, choose .

Height Expression 1

Right-click and choose .

2D Band Structure, -1% Strained

In the window for , type 2D Band Structure, -1% Strained in the text field.

Click to expand the section. From the list, choose .

Then the light hole band (LH).

Evaluation Group 2: LH

In the window, right-click and choose .

In the window for , type Evaluation Group 2: LH in the text field.

Locate the section. From the list, choose .

In the text field, type 2.

Global Evaluation 1: Schrödinger eq.

In the window, expand the node, then click .

In the window for , locate the section.

Table Surface 2: LH (Schr)


Expression	Unit	Description
lambda		LH (Schr)

In the toolbar, click

In the window, right-click and choose .

In the window for , type Table Surface 2: LH (Schr) in the text field.

Locate the section. From the list, choose .

Click to expand the section. From the list, choose .

Then the crystal-field split-off hole band (CH).

Evaluation Group 3: CH

In the window, right-click and choose .

In the window for , type Evaluation Group 3: CH in the text field.

Locate the section. From the list, choose .

Global Evaluation 1: Schrödinger eq.

In the window, expand the node, then click .

In the window for , locate the section.

Table Surface 3: CH (Schr)


Expression	Unit	Description
lambda		CH (Schr)

In the toolbar, click

In the window, right-click and choose .

In the window for , type Table Surface 3: CH (Schr) in the text field.

Locate the section. From the list, choose .

In the toolbar, click

Finally set up a study to compute the 2D band structure in the kx-kz plane using the global equation of the analytic formula.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 4: 2D Dispersion Analytic

In the window for , type Study 4: 2D Dispersion Analytic in the text field.

Locate the section. Clear the check box.

Study 3: 2D Dispersion Schrödinger Eq.

Study 4: 2D Dispersion Analytic

In the window, under right-click and choose .

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Physics interface	Solve for	Equation form
Schrödinger Equation (schr)		Automatic (Eigenvalue)
Global ODEs and DAEs (ge)	√	Not applicable

In the toolbar, click

Similarly create one table for each band to plot the analytic 2D band structure using the plot type for the 3 hole bands. Use the option for the analytic result to compare to the numerical result by overlaying the plotted surfaces. First the heavy hole band (HH).