Code covered by the BSD License

# Chebfun V4

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

Model of a quantum dot array for solar energy

# Model of a quantum dot array for solar energy

Toby Driscoll, May 7, 2011

(Chebfun example ode-eig/SolarQDA.m)

Researchers at the University of Delaware are exploring the use of quantum dot arrays (QDA) to capture solar energy. One simplified model of a QDA is as a 1-D Schrödinger eigenvalue problem. The governing equation is

  -(hbar^2/2m) psi'' + U(x) psi = E psi,

where psi(x) is the wavefunction, hbar is Planck's constant divided by 2 pi, m is an effective mass, E is an allowed energy, and U(x) is a potential with wells representing the quantum dots. We will take psi=0 far from the wells to supply boundary conditions.

Here are some dimensional parameters corresponding to physical experiments. Space is measured in nm and energy is in eV.

hbar = 1.054e-34;
m = [0.067 0.022]*9.109e-31*1.602e-37;  % effective mass of InAs and GaAs

numwell = 4;   % number of wells
width = 6.5;  depth = 0.85;  spacing = 3;   % well parameters


From the well parameters we will find the boundaries of each well, then add endpoints far out to each side.

x = cumsum( [0 repmat([width spacing],1,numwell)] );
x = [ -10*spacing, x(1:end-1), x(end)+9*spacing ];


Now we create the potential U as a piecewise constant function. We use a syntax that is compact but opaque: create a cell array of the constant values, one per subinterval.

LW = 'linewidth'; lw = 1.6;
vals = [repmat([depth 0],1,numwell) depth];  % [ 0 -depth 0 ... -depth 0 ]
vals = mat2cell(vals,1,ones(1,2*numwell+1));  % convert to cell
U = chebfun(vals,x);
plot(U,LW,lw), ylabel('potential')
xlim(x([1 end]))


The effective mass is also piecewise constant, so we go through the same motions.

vals = [repmat(m,1,numwell) m(1)];
vals = mat2cell(vals,1,ones(1,2*numwell+1));  % convert to cell
emass = chebfun(vals,x);


Now we are ready to define the Schrodinger operator.

N = chebop(@(psi) -hbar^2./(2*emass).*diff(psi,2) + U.*psi, x);
N.lbc = 0; N.rbc = 0;


We find the numwell lowest allowed energies and their wavefunctions. Rather than getting one isolated wavefunction "bump" per well, we see that the wells influence one another strongly.

[Psi,E] = eigs(N,numwell,0);
energies = diag(E)
subplot(2,1,1)
plot(U,LW,lw), ylabel('potential'), xlim(x([1 end]))
subplot(2,1,2)
plot(Psi,LW,lw), ylabel('wavefunction'), xlim(x([1 end]))

energies =
0.2376
0.2417
0.2469
0.2513


If we look at probability, we find that the first four modes extend significantly over all wells. This is "delocalization" or quantum tunnelling, and it means the device can conduct electricity.

clf
semilogy(Psi.^2,LW,lw), ylabel('probability'),
axis( [x([2 end-1]) 1e-3 1e-1]  )

Warning: Negative data ignored
Warning: Negative data ignored
Warning: Negative data ignored


In practice, though, the potential wells cannot be fabricated identically. A brief experiment shows that the delocalization deteriorates even with small fluctuations in the well depths. Here we perturb by 2% variance.

vals = [repmat([depth 0],1,numwell) depth];  % [ 0 -depth 0 ... -depth 0 ]

randn('state',1138)
vals(2:2:end) = vals(2:2:end) + 0.017*randn(1,numwell);

vals = mat2cell(vals,1,ones(1,2*numwell+1));  % convert to cell
U = chebfun(vals,x);


Now we find that the wavefunctions extend significantly over just one or two wells.

N.op = @(psi) -hbar^2./(2*emass).*diff(psi,2) + U.*psi;
[Psi,E] = eigs(N,numwell,0);
energies = diag(E)
semilogy(Psi.^2,LW,lw), ylabel('probability'),
axis( [x([2 end-1]) 1e-3 1e-0]  )

energies =
0.2038
0.2321
0.2413
0.2530
Warning: Negative data ignored
Warning: Negative data ignored
Warning: Negative data ignored


References: