In this tutorial, we are going to model optical intersubband transitions between the conduction subbands of a quantum well. In the conduction band, intersubband transitions only couple to the electric field component across the structure (along the growth direction/perpendicular to the plane of the quantum well etc)[ref]. This means that the medium is anisotropic which leads to interesting properties like a shift in the absorption frequency due to a 'depolarisation shift'[ref] as well as standing wave effects[ref]. I could go on about this for a while [my ref] but instead, let's cut straight to the action.
First we will model a 'normal' quantum well structure
import aestimo.aestimo as solver
import aestimo.config as ac
ac.messagesoff = True # turn off logging in order to keep notebook from being flooded with messages.
import aestimo.database as adatabase
import numpy as np
import matplotlib.pyplot as plt
import copy
from pprint import pprint
%matplotlib notebook
WARNING aestimo logs automatically to aestimo.log in the current working directory. INFO using psi_at_inf_cython module WARNING aestimo_eh logs automatically to aestimo.log in the current working directory.
class Structure(object): pass
s1 = Structure() # this will be our datastructure
# TEMPERATURE
s1.T = 15.0 #Kelvin
# COMPUTATIONAL SCHEME
# 0: Schrodinger
# 1: Schrodinger + nonparabolicity
# 2: Schrodinger-Poisson
# 3: Schrodinger-Poisson with nonparabolicity
# 4: Schrodinger-Exchange interaction
# 5: Schrodinger-Poisson + Exchange interaction
# 6: Schrodinger-Poisson + Exchange interaction with nonparabolicity
s1.computation_scheme = 6
# Non-parabolic effective mass function
# 0: no energy dependence
# 1: Nelson's effective 2-band model
# 2: k.p model from Vurgaftman's 2001 paper
s1.meff_method = 2
# Non-parabolic Dispersion Calculations for Fermi-Dirac
s1.fermi_np_scheme = True #needed only for aestimo_numpy2.py
# QUANTUM
# Total subband number to be calculated for electrons
s1.subnumber_e = 3
# APPLIED ELECTRIC FIELD
s1.Fapplied = 0.00 # (V/m)
# GRID
# For 1D, z-axis is choosen
s1.gridfactor = 0.1 #nm
s1.maxgridpoints = 200000 #for controlling the size
# REGIONS
# Region input is a two-dimensional list input.
# | Thickness (nm) | Material | Alloy fraction | Doping(cm^-3) | n or p type |
s1.material =[
[ 20.0, 'AlGaAs', 0.3, 0.0, 'n'],
[ 11.0, 'GaAs', 0, 2e18, 'n'],
[ 20.0, 'AlGaAs', 0.3, 0.0, 'n'],
]
structure1 = s1
# Initialise structure class
model1 = solver.StructureFrom(structure1,adatabase) # structure could also be a dictionary.
#calculate QW states
result1 = solver.Poisson_Schrodinger(model1)
#solver.save_and_plot(result,model) # Write the simulation results in files
solver.QWplot(result1,figno=None) # Plot QW diagram
solver.logger.info("Simulation is finished.")
INFO Total layer number: 3 INFO Total number of materials in database: 23 WARNING The calculation of Vxc depends upon m*, however when non-parabolicity is also considered m* becomes energy dependent which would make Vxc energy dependent. Currently this effect is ignored and Vxc uses the effective masses from the bottom of the conduction bands even when non-parabolicity is considered elsewhere. INFO calculation time 3.39166 s
INFO Simulation is finished.
Now we will calculate the intersubband transitions using the appropriately named 'intersubband_optical_transitions' module. The module contain several different models starting with semiclassical models and finishing with a perturbative quantum mechanical model that accounts for collective interactions between all of the transitions. The theory behind the intersubband transitions is a little more complicated than might be expected, see [?] for details.
import aestimo.intersubband_optical_transitions as isbt
# set dielectric constants
case = 2
if case==1: #scalar dielectric constants
eps_b = 12.90
eps_z = 12.90
elif case==2: #z-dependent dielectric constants
eps_b = 10.364
eps_gaas = 10.364 # @ 16um
eps_algaas = 8.2067
#eps_z = isbt.eps_background_GaAs(model,eps_gaas,eps_algaas)
eps_z = np.zeros(model1.n_max)
position = 0.0 # keeping in nanometres (to minimise errors)
for layer in model1.material:
startindex = solver.round2int(position*1e-9/model1.dx)
position += layer[0] # update position to end of the layer
finishindex = solver.round2int(position*1e-9/model1.dx)
#
matType = layer[1]
if matType == 'GaAs':
eps_z[startindex:finishindex] = eps_gaas
elif matType == 'AlGaAs':
eps_z[startindex:finishindex] = eps_algaas
elif case==3: #w-dependent dielectric constants
#because the zeroth axis is assumed to be the z-axis, our eps_z array must be 2d
pass
#currently the matrix model doesn't cope with frequency dependent dielectric constants
#therefore the classical model would be the best approach (model2) although it seems
#to over-estimate the coupling between the transitions.
#Alternatively, we could resolve the matrix at each frequency (for each value of the
#background dielectric constant which would be accurate but may be quite computationally
#intensive.
elif case==4: #z-dependent and w-dependent dielectric constants
pass
#currently the matrix model doesn't cope with frequency dependent dielectric constants
#therefore the classical model is the best approach (model2) although it seems
#to over-estimate the coupling between the transitions.
#Alternatively, we could resolve the matrix at each frequency (for each value of the
#background dielectric constant which would be accurate but may be quite computationally
#intensive.
#Set thickness of effective medium
#Lperiod = sum([layer[0] for layer in model.material])*1e-9 #m
Lperiod = model1.x_max #m
# Define Linewidth function (or constant)
def linewidth(freq): return 0.1*freq #define linewidth in THz
# Optical Intersubband Transitions
transitions_table,(hdr,units) = isbt.transitions(result1,Lperiod,eps_z,linewidth)
isbt.print_levels(result1)
isbt.print_transitions(transitions_table,hdr,units)
INFO the energy levels\population are (meV) (m**-2): INFO 1019.9 2.1643e+16 INFO 1092.7 3.5687e+14 INFO 1197.5 2.2983e-21 INFO T = 15K INFO the energy levels gaps are INFO (meV) (THz) (um) (wavno) INFO 72.8 17.6 17 587 INFO 178 42.9 6.98 1.43e+03 INFO 105 25.3 11.8 845 INFO Summary of Intersubband Transitions INFO j 0 1 2 INFO ilevel 0 0 1 INFO flevel 1 2 2 INFO dE meV 72.8 178 105 INFO freq THz 17.6 42.9 25.3 INFO lambda um 17 6.98 11.8 INFO wavno cm-1 587 1.43e+03 845 INFO dN 1e11cm-2 @15K 21.3 21.6 0.357 INFO z_if nm -2.65 0.321 -3 INFO f 0.921 0.0328 1.69 INFO Leff nm 9.29 15.4 6.83 INFO S_if nm 0.821 0.203 0.776 INFO S_if_b nm 0.822 0.203 0.777 INFO wp THz 16.2 12.7 2.44 INFO R THz 2.06 0.393 0.362 INFO Lperiod nm 51 51 51 INFO y_if THz 1.76 4.29 2.53 INFO eps_w 10.3 10.3 10.3
# Save default rcParams so we can reset them later
# WARNING: Do not run this cell after changing rcParams, as it will overwrite the
# defaults that we are trying to preserve.
rcdef = plt.rcParams.copy()
newparams = {'savefig.dpi': 120, 'legend.fontsize': 'smaller'}
#newparams = {'axes.labelsize': 14, 'axes.linewidth': 1, 'savefig.dpi': 150,
# 'lines.linewidth': 1.5, #'figure.figsize': (8, 3),
# 'figure.subplot.wspace': 0.4,
# 'ytick.labelsize': 12, 'xtick.labelsize': 12,
# 'ytick.major.pad': 5, 'xtick.major.pad': 5,
# 'legend.fontsize': 12, 'legend.frameon': False,
# 'legend.handlelength': 1.5}
# Update the global rcParams dictionary with the new parameter choices
# Before doing this, we reset rcParams to its default again, just in case
plt.rcParams.update(rcdef)
plt.rcParams.update(newparams)
f1=isbt.plotting_absorption(model1,result1,transitions_table,eps_b,eps_z,linewidth)
INFO calc_wR_multiplasmon: using eigh() solver for Hermitian matrix INFO calc_wR_multiplasmon2: using eigh() solver for Hermitian matrix INFO calc_wR_multiplasmon3: using eigh() solver for Hermitian matrix
In the plot above, the vertical blue lines mark the energy gaps between the quantum well's energy levels. We can also see that apart from the 'naive' model, the three different models are almost in agreement for this case. The result shows that the absorption strength of the lower transition has been shifted to a higher frequency which is an example of the 'depolarisation shift'. It also appears that the first and second transitions have almost merged, in fact we often find that when transitions are close, the higher energy transition will 'steal' the transition strength of the lower one due to the coulombic interactions between the different transitions.
Let's model an even wider quantum well to show what happens when there are many transtions interacting.
structure2 = copy.copy(structure1)
# DOPING
Nd = 1e18
structure2.material =[
[ 10.0, 'AlGaAs', 0.3, 0.0, 'n'],
[ 30.0, 'GaAs', 0, Nd, 'n'],
[ 10.0, 'AlGaAs', 0.3, 0.0, 'n']]
structure2.subnumber_e = 7
# Initialise structure class
model2 = solver.StructureFrom(structure2,adatabase) # structure could also be a dictionary.
INFO Total layer number: 3 INFO Total number of materials in database: 23
This example also allows me to show that sometimes the Poisson-Schrodinger algorithm needs to be optimised. Using the standard parameters, we get weird results. If we run aestimo on the command-line we can see how the energy-state solutions evolve as the algorithm iterates; if we looked at the output for this simulation, we would see that the solutions are drastically changing between each iteration, the algorithm has gone into oscillations rather than tending towards the correct solution. There will probably be a warning that aestimo has reached the maximum number of iterations.
ac.damping = 0.5 # standard value of this parameter
#calculate QW states
result2 = solver.Poisson_Schrodinger(model2)
#solver.save_and_plot(result,model) # Write the simulation results in files
solver.QWplot(result2,figno=None) # Plot QW diagram
solver.logger.info("Simulation is finished.")
WARNING The calculation of Vxc depends upon m*, however when non-parabolicity is also considered m* becomes energy dependent which would make Vxc energy dependent. Currently this effect is ignored and Vxc uses the effective masses from the bottom of the conduction bands even when non-parabolicity is considered elsewhere. WARNING Have reached maximum number of iterations INFO calculation time 27.6719 s
INFO Simulation is finished.
Now, we change one of the parameters in the configuration module so that the algorithm is more 'damped'. This allows the algorithm to find the correct solution.
# Optimise stability of Poisson-Schrodinger Algorithm
ac.damping = 0.3
#calculate QW states
result2 = solver.Poisson_Schrodinger(model2)
#solver.save_and_plot(result,model) # Write the simulation results in files
solver.QWplot(result2,figno=None) # Plot QW diagram
solver.logger.info("Simulation is finished.")
WARNING The calculation of Vxc depends upon m*, however when non-parabolicity is also considered m* becomes energy dependent which would make Vxc energy dependent. Currently this effect is ignored and Vxc uses the effective masses from the bottom of the conduction bands even when non-parabolicity is considered elsewhere. INFO calculation time 12.5607 s
INFO Simulation is finished.
# set dielectric constants
case = 2
if case==1: #scalar dielectric constants
eps_b2 = 12.90
eps_z2 = 12.90
elif case==2: #z-dependent dielectric constants
eps_b2 = 10.364
eps_gaas = 10.364 # @ 16um
eps_algaas = 8.2067
eps_z2 = np.zeros(model2.n_max)
position = 0.0 # keeping in nanometres (to minimise errors)
for layer in model1.material:
startindex = solver.round2int(position*1e-9/model2.dx)
position += layer[0] # update position to end of the layer
finishindex = solver.round2int(position*1e-9/model2.dx)
#
matType = layer[1]
if matType == 'GaAs':
eps_z2[startindex:finishindex] = eps_gaas
elif matType == 'AlGaAs':
eps_z2[startindex:finishindex] = eps_algaas
#Set thickness of effective medium
#Lperiod = sum([layer[0] for layer in model.material])*1e-9 #m
Lperiod2 = model2.x_max #m
# Optical Intersubband Transitions
transitions_table2,(hdr,units) = isbt.transitions(result2,Lperiod2,eps_z2,linewidth)
isbt.print_levels(result2)
isbt.print_transitions(transitions_table2,hdr,units)
#plot output
isbt.plotting_absorption(model2,result2,transitions_table2,eps_b2,eps_z2,linewidth)
plt.gca().set_xlim(xmin=0,xmax=30)
plt.gca().set_ylim(bottom=0,top=0.09)
INFO the energy levels\population are (meV) (m**-2): INFO 1002.7 1.4672e+16 INFO 1015.4 1.1073e+16 INFO 1039 4.2554e+15 INFO 1071.5 3.8099e+08 INFO 1111.1 1.881e-05 INFO 1156.3 1.2895e-20 INFO 1205 3.0188e-37 INFO T = 15K INFO the energy levels gaps are INFO (meV) (THz) (um) (wavno) INFO 12.7 3.07 97.8 102 INFO 36.3 8.78 34.1 293 INFO 68.8 16.6 18 555 INFO 108 26.2 11.4 874 INFO 154 37.1 8.07 1.24e+03 INFO 202 48.9 6.13 1.63e+03 INFO 23.6 5.72 52.5 191 INFO 56.1 13.6 22.1 452 INFO 95.7 23.2 12.9 772 INFO 141 34.1 8.8 1.14e+03 INFO 190 45.9 6.54 1.53e+03 INFO 32.5 7.85 38.2 262 INFO 72.1 17.4 17.2 582 INFO 117 28.4 10.6 946 INFO 166 40.1 7.47 1.34e+03 INFO 39.6 9.59 31.3 320 INFO 84.8 20.5 14.6 684 INFO 134 32.3 9.28 1.08e+03 INFO 45.2 10.9 27.4 364 INFO 93.9 22.7 13.2 758 INFO 48.7 11.8 25.4 393 INFO Summary of Intersubband Transitions INFO j 0 1 2 3 INFO ilevel 0 0 0 0 INFO flevel 1 2 3 4 INFO dE meV 12.7 36.3 68.8 108 INFO freq THz 3.07 8.78 16.6 26.2 INFO lambda um 97.8 34.1 18 11.4 INFO wavno cm-1 102 293 555 874 INFO dN 1e11cm-2 @15K 3.6 10.4 14.7 14.7 INFO z_if nm -6.43 -0.00719 -0.627 0.0432 INFO f 0.927 3.32e-06 0.0478 0.000358 INFO Leff nm 21.6 27.2 30.6 33.5 INFO S_if nm 2.07 0.573 0.269 0.156 INFO S_if_b nm 2.07 0.574 0.269 0.156 INFO wp THz 4.64 7.03 7.87 7.52 INFO R THz 0.964 0.0031 0.442 0.0382 INFO Lperiod nm 50 50 50 50 INFO y_if THz 0.307 0.878 1.66 2.62 INFO eps_w 9.28 9.28 9.28 9.28 INFO j 4 5 6 7 INFO ilevel 0 0 1 1 INFO flevel 5 6 2 3 INFO dE meV 154 202 23.6 56.1 INFO freq THz 37.1 48.9 5.72 13.6 INFO lambda um 8.07 6.13 52.5 22.1 INFO wavno cm-1 1.24e+03 1.63e+03 191 452 INFO dN 1e11cm-2 @15K 14.7 14.7 6.82 11.1 INFO z_if nm -0.214 0.0812 -6.48 0.0331 INFO f 0.0124 0.00236 1.75 0.000109 INFO Leff nm 36.2 40.3 14 21.5 INFO S_if nm 0.102 0.0694 1.72 0.469 INFO S_if_b nm 0.102 0.0695 1.72 0.469 INFO wp THz 7.24 6.85 7.94 8.15 INFO R THz 0.225 0.0981 1.82 0.0183 INFO Lperiod nm 50 50 50 50 INFO y_if THz 3.71 4.89 0.572 1.36 INFO eps_w 9.28 9.28 9.28 9.28 INFO j 8 9 10 11 INFO ilevel 1 1 1 2 INFO flevel 4 5 6 3 INFO dE meV 95.7 141 190 32.5 INFO freq THz 23.2 34.1 45.9 7.85 INFO lambda um 12.9 8.8 6.54 38.2 INFO wavno cm-1 772 1.14e+03 1.53e+03 262 INFO dN 1e11cm-2 @15K 11.1 11.1 11.1 4.26 INFO z_if nm -0.747 0.112 -0.313 -6.57 INFO f 0.0943 0.0031 0.0328 2.48 INFO Leff nm 26.5 31.2 35.4 10.3 INFO S_if nm 0.223 0.129 0.0844 1.69 INFO S_if_b nm 0.223 0.129 0.0846 1.69 INFO wp THz 7.34 6.77 6.36 7.3 INFO R THz 0.539 0.0977 0.318 1.71 INFO Lperiod nm 50 50 50 50 INFO y_if THz 2.32 3.41 4.59 0.785 INFO eps_w 9.28 9.28 9.28 9.28 INFO j 12 13 14 15 INFO ilevel 2 2 2 3 INFO flevel 4 5 6 4 INFO dE meV 72.1 117 166 39.6 INFO freq THz 17.4 28.4 40.1 9.59 INFO lambda um 17.2 10.6 7.47 31.3 INFO wavno cm-1 582 946 1.34e+03 320 INFO dN 1e11cm-2 @15K 4.26 4.26 4.26 3.81e-07 INFO z_if nm 0.0962 -0.811 0.201 -6.7 INFO f 0.00118 0.136 0.0119 3.14 INFO Leff nm 17.8 23.2 29.7 8.38 INFO S_if nm 0.442 0.208 0.115 1.7 INFO S_if_b nm 0.442 0.208 0.115 1.7 INFO wp THz 5.56 4.87 4.3 0.00242 INFO R THz 0.0374 0.402 0.119 0.000577 INFO Lperiod nm 50 50 50 50 INFO y_if THz 1.74 2.84 4.01 0.959 INFO eps_w 9.28 9.28 9.28 9.28 INFO j 16 17 18 19 INFO ilevel 3 3 4 4 INFO flevel 5 6 5 6 INFO dE meV 84.8 134 45.2 93.9 INFO freq THz 20.5 32.3 10.9 22.7 INFO lambda um 14.6 9.28 27.4 13.2 INFO wavno cm-1 684 1.08e+03 364 758 INFO dN 1e11cm-2 @15K 3.81e-07 3.81e-07 1.88e-20 1.88e-20 INFO z_if nm 0.182 -0.891 -6.86 0.307 INFO f 0.00499 0.187 3.75 0.0157 INFO Leff nm 15.5 21.1 7.21 14.3 INFO S_if nm 0.431 0.201 1.74 0.422 INFO S_if_b nm 0.431 0.201 1.74 0.422 INFO wp THz 0.00178 0.00153 5.8e-10 4.12e-10 INFO R THz 2.3e-05 0.000141 1.4e-10 9.06e-12 INFO Lperiod nm 50 50 50 50 INFO y_if THz 2.05 3.23 1.09 2.27 INFO eps_w 9.28 9.28 9.28 9.28 INFO j 20 INFO ilevel 5 INFO flevel 6 INFO dE meV 48.7 INFO freq THz 11.8 INFO lambda um 25.4 INFO wavno cm-1 393 INFO dN 1e11cm-2 @15K 1.29e-35 INFO z_if nm -7.12 INFO f 4.36 INFO Leff nm 6.47 INFO S_if nm 1.8 INFO S_if_b nm 1.8 INFO wp THz 1.6e-17 INFO R THz 3.96e-18 INFO Lperiod nm 50 INFO y_if THz 1.18 INFO eps_w 9.28
INFO calc_wR_multiplasmon: using eigh() solver for Hermitian matrix INFO calc_wR_multiplasmon2: using eigh() solver for Hermitian matrix INFO calc_wR_multiplasmon3: using eigh() solver for Hermitian matrix
(0, 0.09)
So when there are multiple transitions, we can see that the matrix model gives a very different result compared to treating the transitions independently. In general, we find that the frequency of the absorption will be at a higher frequency than expected and most of the transition strength will be concentrated into a single transition.