Department of physics, NTNU TFY4340 Mesoscopic Physics Spring 2010

Solution to Exercise 5

a) We have

E(α) = \frac{∫ Φ_{0}^{∗}HΦ_{0}dz}{∫ Φ_{0}^{∗}Φ_{0}dz} ,

with

Φ_{0}(z, α) = z e^{−αz/2}
H = −\frac{h̄^{2}}{2m^{∗}} \frac{d^{2}}{dz^{2}} + F z

and we must integrate from z = 0 to z = ∞. Relevant integrals are of the type

∫_{0}^{∞} z^{n} e^{−αz} dz = (−1)^{n} \frac{d^{n}}{dα^{n}} ∫_{0}^{∞} e^{−αz} dz = \frac{n!}{α^{n+1}} ,

with n between 1 and 3. Hence,

E(α) = \frac{∫_{0}^{∞} dzze^{−αz/2} [\frac{h̄^{2}}{2m^{∗}} (α − α^{2}z/4) + F z^{2}] e^{−αz/2}}{∫_{0}^{∞} dzz^{2}e^{−αz}}

= \frac{∫_{0}^{∞} dze^{−αz} [\frac{h̄^{2}α}{2m^{∗}} z − \frac{h̄^{2}α^{2}}{8m^{∗}} z^{2} + F z^{3}]}{∫_{0}^{∞} dzz^{2}e^{−αz}}
= \frac{\frac{h̄^{2}α}{2m^{∗}} · \frac{1}{α^{2}} − \frac{h̄^{2}α^{2}}{8m^{∗}} · \frac{2}{α^{3}} + \frac{6F}{α^{4}}}{\frac{2}{α^{3}}}
= \frac{h̄^{2}α^{2}}{8m^{∗}} + \frac{3F}{α}

and the variational estimate based on this trial wave function is found by setting dE/dα = 0, i.e.,

\frac{h̄^{2}α}{4m^{∗}} − \frac{3F}{α^{2}} = 0,

or

α = (\frac{12m^{∗}F}{h̄^{2}} )^{1/3} .

Inserting this value of α into E yields, after some algebra,

E_{min} = (\frac{h̄^{2}}{2m^{∗}} )^{1/3} · F^{2/3} · \frac{9}{2 · 6^{1/3}} ≃ (\frac{h̄^{2}}{2m^{∗}} )^{1/3} · F^{2/3} · 2.48,

which is slightly higher than the exact result, where the numerical factor was 2.34. b) The numerical solution of the Schrödinger equation in the matlab m-script given in the exercise, is based on a straightforward discretization of the second derivative that appears in the kinetic energy part of the Hamiltonian:

\frac{d^{2}Φ(z)}{dz^{2}} ≃ \frac{Φ(z − s) − 2Φ(z) + Φ(z + s)}{s^{2}},

where s is the step size for the z coordinate. Then, with Φ_{j} = Φ(z_{j}) = Φ((j − N − 1)s), the discretized version of the Schrödinger equation is

−\frac{h̄^{2}}{2m^{∗}s^{2}} (Φ_{j+1} − 2Φ_{j} + Φ_{j−1}) + V_{j}Φ_{j} = EΦ_{j}.

Here, V_{j} = V (z_{j}) = V ((j − N − 1)s), i.e., the value of the potential in position z_{j}. If we choose our potential profile in the range 1 ≤ j ≤ 2N + 1, the potential is infinite in j ≤ 0 and j ≥ 2N + 2, i.e., the wavefunction is zero in z ≤ −(N + 1)s and in z ≥ (N + 1)s. In matrix form, we have

[ \begin{matrix}\frac{h̄^{2}}{m^{∗}s^{2}} + V_{1} & −\frac{h̄^{2}}{2m^{∗}s^{2}} \\ −\frac{h̄^{2}}{2m^{∗}s^{2}} & \frac{h̄^{2}}{m^{∗}s^{2}} + V_{2} & −\frac{h̄^{2}}{2m^{∗}s^{2}} \\ & −\frac{h̄^{2}}{2m^{∗}s^{2}} & \frac{h̄^{2}}{m^{∗}s^{2}} + V_{3} & −\frac{h̄^{2}}{2m^{∗}s^{2}} \\ & & ... & ... & ... \end{matrix} ] [ Φ_{1} Φ_{2} Φ_{3} ⋮ ] = E [ Φ_{1} Φ_{2} Φ_{3} ⋮ ] .

This tridiagonal Hamiltonian matrix H may be diagonalized with a matrix S:

S^{−1}HS = D.

The elements of the diagonal matrix D are now the energy eigenvalues E_{n} (n = 1, ⋯, 2N + 1) of H, whereas the columns of the transformation matrix S represent the corresponding eigenfunctions Φ^{(n)}. In Matlab, the command

[S,D] = eig(H);

will return the desired matrices S and D.

In the matlab program

triangular_well.m

a smooth potential is used between Z_{min} = −50 nm and the interface at z = 0, with V = 270 meV for large negative z and V = 300 meV at the interface. For positive z, the potential rises linearly, with ”force factor” F = 10 meV/nm. Outside the interval −50 < z < 50 nm, we have V = ∞ and Φ = 0. The resulting 5 lowest eigenvalues are 75, 141, 195, 243, and 273 meV (see figure next page), so the distance between the ground state and the first excited state is 66 meV. In the lectures, we found that a 2D carrier density of 4 · 10^{11} cm^{−2} corresponds to a Fermi level 14 meV above the lowest 2D subband edge (at 75 meV). Hence, the Fermi level is well below the first excited 2D subband edge (at 141 meV).

The corresponding wavefunctions Φ^{(1)}, ⋯, Φ^{(5)} look like this: In both figures, the potential profile is also included. In the latter figure, the potential has been scaled, in order to fit in with plots of the wavefunctions.