Self-Consistent Schrödinger–Poisson Solution in a Doped Quantum Well
Introduction
Solving the Schrödinger–Poisson (SP) equation provides crucial insight into the behavior of devices such as High Electron Mobility Transistors (HEMTs) [1] or Resonant Tunneling Diodes (RTDs) [2]. The goal of this example is to show how to solve the SP equation in the Wolfram Language for a Quantum Well (QW). The QW is made from the junction 
, 
being the material sandwiched between 
. The 
region is 
-doped, meaning it has excess conduction electrons, all of which are assumed to occupy the ground state of the QW—an assumption valid at low temperatures, e.g. as low as 
in some cases [3]. This means that the charge carriers will be distributed according to the probability density of the ground state. However, this redistribution of the charge in turn changes the potential energy experienced by the electrons confined in the QW. So, arriving at a wavefunction that is consistent with the potential energy field, which has been modified by the wavefunction itself, is what is referred to as a self-consistent solution. Getting such a solution for the SP equation consists of an iterative approach. First, the Schrödinger equation is solved for the finite barrier potential, and the resulting wavefunction is used to calculate the charge density. Then, Poisson's equation is solved to obtain the electric potential—and therefore the potential energy—which is subsequently inserted back into the Schrödinger equation. Repeating this process iteratively until the energy eigenvalue converges yields a self-consistent solution of the SP equation. The calculations shown here are based on an example on the third chapter of the book Quantum Wells, Wires and Dots by Paul Harrison [4].
The QW for our problem has a length of 100 Å, and its surrounded by 200 Å barriers. The potential barrier that defines the QW is proportional to the difference between the energy gaps of each material. The confinement potential is then defined as a Piecewise function that depends on the energy gap, which in turn depends on the aluminum concentration x.
Eg[x_] := 1426 + 1247x
𝒱[z_] := Piecewise[{{0, 200 <= z <= 300}, {0.66*(Eg[0.2] - Eg[0]), True}}]Note that in the definition of 
, the 0.2 in the Eg[0.2] - Eg[0] term comes from the fact that the aluminum concentration in the barriers made of 
is 0.2. A plot of the finite barrier potential is shown below. This would be the only potential energy term in the Schrödinger equation if no self-consistent solution were sought.
Plot[𝒱[z], {z, 0, 500}, ...]Note that 
is chosen as the spatial variable instead of the more common 
. This is because 
often denotes the growth axis—i.e. the crystallographic direction perpendicular to the substrate's surface—during the QW fabrication by epitaxy.
Geometry
Next, the geometry is defined, in this case is a very simple one: a line. It is possible to control the discretization of the geometry by creating a mesh and choosing the "MaxCellMeasure" accordingly.
Needs["NDSolve`FEM`"]Ω = Line[{{0}, {500}}]mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 10];
MeshRegion[mesh]Parameters
After setting up the geometry, some parameters need to be defined. In particular, as mentioned above, the material at the center of the QW (
) is 
-doped, which, in short, means that it has excess conduction electrons. These extra conduction electrons are present because impurities have been introduced into the GaAs material—i.e. atoms that do not naturally belong in the 
crystal structure, such as silicon (
), which have one more valence electron than gallium (
). Each impurity atom donates one free electron when ionized. Therefore, it is possible to determine the number of conduction electrons in the QW by knowing the concentration of donor atoms introduced into the material.
The density of dopants in the 
material is assumed to be constant, with a magnitude of 
. And it is only nonzero in the inside of the well.
d[z_] := Piecewise[{{QuantityMagnitude[Quantity[2*10^18, "Centimeters"^(-3)], "SIBase"], 200 <= z <= 300},
{0, True}}]Plot[d[z], {z, 0, 500}, Rule[...]]Integrating the dopant density along the 
axis yields the total number of electrons. Note that a factor of 
is introduced in the integral such that the units are consistent, since the units of the geometry are angstroms.
𝒩 = Subsuperscript[∫, -∞, ∞] d[z]10^-10 ⅆz//NNote the assumption that all the donors in the 
region are ionized, meaning their electrons are in the conduction band of 
. However, a more sophisticated model could also be considered, one that accounts for partially ionized donors or even nonuniform donor distributions.
Furthermore, the electron charge is defined below, along with the charge density function, which depends on the wavefunction. As mentioned earlier, the electrons will be distributed in the QW according to the probability density function 
such that the charge density generated by the electrons is of the form 
. However, the ionized donors have a remaining positive charge. This means that the total charge in a region is the combination of the charge contributed by the electrons and the positive charge resulting from donors that have lost an electron. This is where the subtraction in the definition for 
below comes from.
q = -QuantityMagnitude[Quantity[1, "ElementaryCharge"], "SIBase"];sfψ = QuantityMagnitude[Quantity[1, (1/Sqrt["Angstroms"])], (1/Sqrt["Meters"])]ρ[z_, ψunscaled_] := Module[{ψ},
ψ = sfψ ψunscaled[z];
q(𝒩 ψψ - d[z])
]Note also that in the definition of 
a scaling factor for the wavefunctions is introduced. Since the geometry is defined in units of angstroms, the units of the wavefunction will not be the traditional 1/
, but in 
. So a factor is introduced accordingly, in a way that leaves the value for the charge density in SI base units.
Setting Up the Model
Next, the variables and parameters for both the Schrödinger and Poisson equations are defined. An important thing to note here is that, for quantum mechanical problems—which often involve small dimensions, on the order of angstroms in this case—defining the geometry or mesh using SI base units is not a good idea. This is because creating a mesh or geometry with dimensions on the order of 
might cause numerical issues. For more details on how large differences in the scale of the geometry and the material parameters can spoil the solution, take a look at this example. A better alternative is to define the geometry in the target unit, just like the geometry defined above in angstroms, and then use the "ScaleUnits" parameter in the PDEComponent function of your choice. In this particular case, the parameter is set up like this: "ScaleUnits" -> {"Meters" -> "Angstroms"}. Note that if "ScaleUnits" is used, all the parameters that have units must include those units. Internally, all the units are converted to SI base units and then another conversion takes place. If "ScaleUnits" is specified, the quantities with meters in them are converted to angstroms. This ensures that the PDE coefficients have units consistent with the geometry.
To generate the PDEs of interest, the functions SchrodingerPDEComponent and ElectrostaticPDEComponent will be used. When the ElectrostaticPDEComponent is used without specifying any polarization vector 
, you obtain a Poisson equation whose solution is the electrostatic potential. Then, the variables and parameters needed are defined below.
schvars = {ψ[z], {z}};
schpars = <|"SchrodingerPotential" -> Quantity[𝒱[z], "Millielectronvolts"], "ReducedPlanckConstant" -> Quantity[1, "ReducedPlanckConstant"], "ScaleUnits" -> {"Meters" -> "Angstroms"}, "Mass" -> 0.066Entity["Particle", "Electron"][EntityProperty["Particle", "Mass"]]|>;poissvars = {𝒱ρ[z], {z}};
poisspars = <|"RelativePermittivity" -> 13.18, "VolumeChargeDensity" -> Quantity[ρ[z, fun], ("Coulombs"/"Meters"^3)], "ScaleUnits" -> {"Meters" -> "Angstroms"}|>;Note that here all parameters, except those without physical dimensions, use the InputForm Quantity[magnitude,units]. Also, the value chosen for the "VolumeChargeDensity" parameter will become clearer later.
Since the doped region is the central part of the QW, it can be expected that very far from the well, where there is no excess charge, the electric potential will decay to zero. A practical way to model this decay is to set the electric potential to zero at the external boundaries of the defined geometry. This can be done with the ElectricPotentialCondition.
poissBc = {ElectricPotentialCondition[z == 0 || z == 500, poissvars, poisspars, <|"ElectricPotential" -> 0|>]}Solving the Problem Step by Step
Now, everything is set up to get a self-consistent SP solution. The process consists of solving the Schrödinger eigensystem with only the finite barrier potential, and with that first solution computing the charge density. Next, that charge density is used to compute the solution of Poisson's equation; the resulting electric potential is used to recompute the electric potential energy by just multiplying by the charge 
, which goes back into the Schrödinger equation to be solved again. This process is repeated until convergence in the energy eigenvalue is reached. A flowchart for this algorithm is shown in the image below.
Flowchart for the SP solution algorithm.
A program can be written to perform this algorithm with just one function call. But, for clarity, the first two steps of that process will be shown first.
The first step is to solve the Schrödinger equation with just the finite barrier confinement potential. As mentioned earlier, the focus is only on the first eigenvalue and first energy eigenvalue.
ℋ = SchrodingerPDEComponent[schvars, schpars]{{en}, {fun}} = NDEigensystem[ℋ, ψ, {z}∈mesh, 1]Plot[fun[z]fun[z], {z, 0, 500}, Rule[...]]Also, the charge density computed using the solution from the Schrödinger equation can be plotted.
Plot[(ρ[z, fun] /Abs[q]), {z, 0, 500}, ...]The charge density plot shows that the most negative charge is in the center of the well region, while positive charge is found near the barriers. This is expected since there are fewer electrons near the barriers, and therefore the background positive charge from the ionized donors is observed. Also, in the barrier regions, which are not doped, there is a negative charge that quickly goes to zero, following the decaying behavior of the wavefunction. This is also expected since there is a nonzero chance of finding electrons outside the QW; consequently, a portion of the electrons is found in the barrier regions.
Next, with the charge density available, you can solve the Poisson equation. First, the Poisson operator is defined with the ElectrostaticPDEComponent, and then NDSolveValue computes the solution using the previously defined boundary conditions.
poissOp = ElectrostaticPDEComponent[poissvars, poisspars]Note that the parameters of the Poisson equation use the function for the charge density defined earlier, which has the already obtained wavefunction as an input.
poissparspoissSol = NDSolveValue[{poissOp == 0, poissBc}, 𝒱ρ, z∈mesh]Before plotting the solution for Poisson's equation, its necessary to think about the units introduced earlier in the PDE coefficients. As you might remember, the "ScaleUnits" parameter was specified to {"Meters" -> "Angstroms"}. So the solution's length units are angstroms instead of meters. For that reason, a scale parameter is introduced to get the solution in SI base units.
sf = QuantityMagnitude[Quantity[1, ("Angstroms"^2"Kilograms"/"Amperes""Seconds"^3)], ("Millielectronvolts"/"Coulombs")]Note that the electric potential units have been transformed from 
in SI base form (
), but with the length unit expressed in angstroms instead of meters (
), into 
.
Plot[q sf poissSol[z], z∈Ω, ...]Now, the potential energy obtained from Poisson's equation can be added to the already existing confinement potential of the QW, and you can solve the Schrödinger equation again, following the iterative algorithm discussed earlier. But instead of doing each step manually, several helper functions are defined to do the calculation in a more convenient way.
Helper Functions
First, a helper function to convert the energy eigenvalues to 
.
toMilliElectronVolts[x_] := UnitConvert[Quantity[x, "Kilograms"("Angstroms"^2/"Seconds"^2)], "Millielectronvolts"]Also, a helper function is defined that takes the potential energy obtained from the Poisson's equation solution. Here, the parameter "SchrodingerPotential" is updated in each function call. Then the new PDE is generated, and the corresponding eigensystem is solved.
solveSchrodinger[Uρ_] := Module[{pars, ℋ, en, fun},
schpars["SchrodingerPotential"] = Quantity[𝒱[z], "Millielectronvolts"] + Uρ;
ℋ = SchrodingerPDEComponent[schvars, schpars];
{{en}, {fun}} = NDEigensystem[ℋ, ψ, {z}∈mesh, 1];
{en, fun}
]Likewise, a helper function that solves the Poisson equation is defined, and it takes as a input the wavefunction obtained in the previous step. Here, the redefined parameter is "VolumeChargeDensity":
solvePoisson[ψn_] := Module[{op},
poisspars["VolumeChargeDensity"] = Quantity[ρ[z, ψn], ("Coulombs"/"Meters"^3)];
op = ElectrostaticPDEComponent[poissvars, poisspars];
NDSolveValue[{op == 0, poissBc}, 𝒱ρ, z∈mesh]
]Finally, a helper function with a script that applies the SP algorithm introduced earlier is defined. A while loop is used to iterate until the difference between two following eigenvalues is less than a certain tolerance, which here is defined to be 
.
schrodingerPoissonStep[previousStep_] := Module[
{en, fun, old𝒱ρ, 𝒱ρ, schVal, schFun},
{en, fun, old𝒱ρ} = previousStep;
𝒱ρ = solvePoisson[fun];
{schVal, schFun} = solveSchrodinger[q sf 𝒱ρ[z]Quantity[1, "Millielectronvolts"]];
{schVal, schFun, 𝒱ρ}
]Note that the constants q and sf are defined outside the solveSchrodinger helper function. They are used as multiplicative factors in the input to the function, even though they are not passed explicitly as arguments.
schrodingerPoissonSolver[] := Module[
{en, fun, steps, energies, wavefunction, 𝒱ρ, tol = 10^-5},
{en, fun} = solveSchrodinger[0Quantity[1, "Millielectronvolts"]];
steps = NestWhileList[
schrodingerPoissonStep, {en, fun, None}, Abs[#1[[1]] - #2[[1]]] > tol&, 2];
{energies, wavefunction, 𝒱ρ} = Transpose[steps];
{toMilliElectronVolts[energies], Last[wavefunction], Last[𝒱ρ]}
]Solving the Whole Problem
Now, with a simple function call it is possible to get the self-consistent solution for the SP equation. In particular, the energies for each iteration are obtained, and the final solutions for both the Schrödinger's and Poisson's equations, respectively.
{energies, fun, 𝒱sol} = schrodingerPoissonSolver[]First, the energies for each iteration are plotted, where the final energy eigenvalue is higher than the initial guess, proving that an SP analysis is necessary to get a more accurate energy eigenvalue for this problem. Furthermore, the iteration process quickly converges after three iterations, showing minimal changes in the energy eigenvalue for later iterations.
ListLinePlot[energies, Rule[...]]Finally, the wavefunction obtained from the SP solution, along with the modified confinement potential of the QW—i.e. the sum of the band-edge potential and the potential energy from the Poisson equation—is plotted. A scaling factor of 10000 is applied for visualization purposes, so that the wavefunction can be displayed alongside the QW potential energy.
scaleWaveFunction = 10000;
Plot[{scaleWaveFunction fun[z] fun[z], 𝒱[z] + q sf 𝒱sol[z]}, {z, 0, 500}, ...]Conclusions
A self-consistent solution for the Schrödinger–Poisson equation describing an 
-doped QW was obtained. The change in the energy eigenvalue is appreciable, indicating that it is necessary to include the effect of charge redistribution to obtain an accurate energy value, compared to solving only the Schrödinger equation. On the other hand, a slight modification of the confinement potential—shown in orange in the plot above—is observed, but its effect on the wavefunction appears negligible.
References
1. Wu, S., Webster, R. T., & Anwar, A. F. M. (1999). "Physics-Based Intrinsic Model for AlGaN/GaN HEMTs." MRS Internet Journal of Nitride Semiconductor Research, 4(S1), 775–780. https://doi.org/10.1557/S1092578300003409.
2. Sun, J. P., & Haddad, G. I. (1998). "Self‐Consistent Scattering Calculation of Resonant Tunneling Diode Characteristics." VLSI Design, 6(1–4), 83–86. https://doi.org/10.1155/1998/78412.
3. Goñi, A. R., Haboeck, U., Thomsen, C., Eberl, K., Reboredo, F. A., Proetto, C. R., & Guinea, F. (2002). "Exchange Instability of the Two-Dimensional Electron Gas in Semiconductor Quantum Wells." Physical Review B, 65(12), 121313. https://doi.org/10.1103/PhysRevB.65.121313.
4. Harrison, P. (2005). Quantum Wells, Wires and Dots. John Wiley & Sons, Ltd. https://doi.org/10.1002/0470010827.