COMSOL 6.4 - Self-Consistent Schrödinger–Poisson Results for a GaAs Nanowire

Self-Consistent Schrödinger–Poisson Results for a GaAs Nanowire

This benchmark model simulates a GaAs nanowire using the self-consistent Schrödinger–Poisson theory to compute the electron density and the confining potential profiles. The predefined Schrödinger–Poisson multiphysics coupling feature is combined with the dedicated Schrödinger–Poisson study type to provide streamlined procedure for model setup and automated creation of self-consistent iterations with tunable parameters for optimizing the convergence under different conditions. The computed electron density and confining potential profiles compare well with the figures in the reference paper, with both profiles reproducing the distinct Friedel-type spatial oscillations at low temperatures.

Introduction

As the semiconductor fabrication technology progresses to miniaturize devices, the effect of quantum confinement becomes more and more important. This tutorial follows the Schrödinger–Poisson theory described in Ref. 1 to compute the self-consistent electron density distribution and confining potential profile for a GaAs nanowire.

Model Definition

The nanowire is assumed to be infinitely long and axially symmetric by the authors of Ref. 1. Thus a 1D axisymmetric geometry is used. Not all simulation parameters are detailed in the paper, especially the material properties. When a parameter is not available, a typical value found in the literature is used in the model.

The radius of the nanowire is assumed to be 50 nm. The electron effective mass is assumed to be 0.067 times free electron mass (as suggested by the Fermi temperature computation result shown in the paper), and the dielectric constant is assumed to be 12.9. The Fermi energy level in the model is set to 0 V and the electric potential at the wall to −0.7 V, to be consistent with the Fermi level pinning boundary condition described in the paper.

We model the case of 2·1018 cm−3 uniform ionized dopants at the temperature of 10 K, to compare with the Fig. 2 and 3 in the paper.

The numbers above are entered as global parameters in the model.

Following the approach of the paper, we first solve for the Thomas–Fermi approximate solution. Then use it as the initial condition for the fully coupled Schrödinger–Poisson equation.

See the comments in the section Modeling Instructions for more discussions on the model construction, solution process, and result visualization.

Results and Discussion

Figure 1 and Figure 2 show the electron density and confining potential, which compare well with Fig. 2 (lower panel) and 3 in Ref. 1, respectively.

Figure 1: Zoomed-in plot to show the Friedel-type spatial oscillations.

Figure 2: Zoomed-out plot with partial contributions from azimuthal quantum numbers.

Reference

1. J.H. Luscombe, A.M. Bouchard, and M. Luban, “Electron confinement in quantum nanostructures: Self-consistent Poisson-Schrödinger theory,” Phys. Rev. B, vol. 46, no. 16, p. 10262, 1992.

Semiconductor_Module/Quantum_Systems/schrodinger_poisson_nanowire

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

Since we will first solve for the Thomas–Fermi approximate solution, select the study type. We will add the dedicated study later.

In the tree, select > .

Click

Geometry 1

The Model Wizard starts the COMSOL Desktop at the node. 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 .

Load global parameters and local (spatially dependent) variables from files.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file schrodinger_poisson_nanowire_parameters.txt.

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file schrodinger_poisson_nanowire_variables.txt.

The cross section of the nanowire is a simple radial line segment in the 1D axisymmetric geometry.

Geometry 1

Interval 1 (i1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Coordinates (nm)
0
R0

Enter the global parameter m for the azimuthal quantum number into the corresponding input field in the Schrödinger Equation physics node, so that later we can use the Auxiliary sweep option in the study settings to solve for multiple values of the quantum number.

Schrödinger Equation (schr)

In the window, under click .

In the window for , locate the section.

In the m text field, type m.

Enter the global parameter meff for the effective mass.

Effective Mass 1

In the window, under > click .

In the window for , locate the section.

In the meff,e,11 text field, type meff.

Since the electron potential energy will be provided by the Schrödinger-Poisson multiphysics coupling, set it to zero in the physics feature.

Electron Potential Energy 1

In the window, click .

In the window for , locate the section.

From the Ve list, choose . In the associated text field, type 0.

Set up the boundary condition for the electrostatics physics at the wall of the nanowire given by the Fermi level pinning at 0.7 V as discussed in the reference paper. Enable Advanced Physics Options to choose the option to obtain accurate displacement field computation for the electrical neutrality check later after solving.

Electrostatics (es)

In the window, under click .

Charge Conservation in Fluids 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Electric Potential 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the V0 text field, type V_R0.

Click the

button in the toolbar.

In the dialog, in the tree, select the checkbox for the node > .

Click .

In the window for , click to expand the section.

From the list, choose .

Add the space charge density contribution from the uniform distribution of ionized dopants.

Space Charge Density 1: Ionized dopants

In the toolbar, click

and choose .

In the window for , type Space Charge Density 1: Ionized dopants in the text field.

Locate the section. From the list, choose .

Locate the section. In the ρv text field, type e_const*Nd.

Add the space charge density contribution from the electrons under the Thomas–Fermi approximation. This contribution will be replaced by a statistically weighted sum of the probability densities of the eigenstates when solving the coupled Schrödinger-Poisson equation.

Space Charge Density 2: Thomas Fermi

In the toolbar, click

and choose .

In the window for , type Space Charge Density 2: Thomas Fermi in the text field.

Locate the section. From the list, choose .

Locate the section. In the ρv text field, type rho_TF.

Create a blank material and enter the dielectric constant of GaAs defined as a global parameter previously.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	epsr	1	Basic

For the section, the default option of Fermi–Dirac statistics with parabolic band approximation is suitable to compute the statistically weighted sum of the probability densities of the eigenstates for the electrons in the nanowire. Enter the previously defined global parameters for the temperature, Fermi energy level, and effective mass. To take into account the pairs of degenerate azimuthal quantum numbers, use the formula 1+(m>0) for the degeneracy factor, which evaluates to 1 for m=0 and 2 for m>0.

Multiphysics

Schrödinger–Poisson Coupling 1 (schrp1)

In the window, under > click .

In the window for , locate the section.

In the T text field, type T.

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

In the md text field, type meff.

In the gi text field, type 1+(m>0).

For the section, leave the charge number at the default of -1, appropriate for electrons. The formula for the is based on the limiting case of high temperature, nondegenerate systems. For systems at lower temperatures, convergence of the iteration can be accelerated with a judicious increase of the tuning parameter α. The default expression for the global error variable (schrp1.max(abs(V-schrp1.V_old)))/1[V] computes the max difference between the electric potential fields from the two most recent iterations, in the unit of V.

Locate the section. In the α text field, type 4.

Create a reasonably fine mesh.

Mesh 1

Edge 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 200.

The first study solves for the Thomas–Fermi approximate solution, which is fully configured in the Electrostatics interface, so we clear the checkboxes for the Schrödinger Equation and the Schrödinger-Poisson Coupling in the study settings.

Study 1: Thomas-Fermi

In the window, click .

In the window for , type Study 1: Thomas-Fermi in the text field.

Step 1: Stationary

In the window, under click .

In the window for , locate the section.

In the column of the table, under , clear the checkbox for .

In the column of the table, under > , clear the checkbox for .

In the toolbar, click

Results

Electric Potential (es)

Now add the dedicated Schrödinger-Poisson study.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Schrödinger–Poisson

For a completely new problem, it is often necessary to use the default search option to find the range of the eigenenergies. Once the range is found, switch to the search option with appropriate settings for the range and number of eigenvalues, in order to ensure that all significant eigenstates are found by the solver. Here the estimated energy range is between -0.15 and 0.05 eV. Disable the contribution to the charge density from the Thomas–Fermi approximation as discussed earlier.

In the window for , locate the section.

From the list, choose .

Find the subsection. In the text field, type -0.15.

In the text field, type 0.05.

Locate the section. Select the checkbox.

In the tree, select > > .

Click

The default iteration method updates a table displaying the history of a global error variable after each iteration during the solution process. The built-in global error variable schrp1.global_err computes the max difference between the electric potential fields from the two most recent iterations, in the unit of V, as already configured in the Schrödinger-Poisson Coupling multiphysics node earlier. (Note that the prefix for the variable, in this case , should match the input field of the Schrödinger-Poisson Coupling multiphysics node.) Setting the tolerance to 1e-6 thus means the iteration ends after the max difference is less than 1 uV. Use the Thomas–Fermi approximate solution as the initial condition for this study.

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

From the list, choose .

Use the Auxiliary sweep functionality to solve for a list of nonnegative azimuthal quantum numbers. The negative ones are taken into account using the formula 1+(m>0) for the degeneracy factor as discussed earlier. The dedicated solver sequence automatically performs the statistically weighted sum of the probability densities of all the eigenstates.

Click to expand the section. Select the checkbox.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
m (Azimuthal quantum number)	0 1 2 3 4 5 6

Get initial value to generate default plots to be modified and plotted while solving.

In the window, click .

In the window for , type Study 2: Schrödinger-Poisson in the text field.

In the toolbar, click

Results

Probability Density (schr)

For the plot group to show the probability density of each eigenstate, point the dataset to the solution that stores the wave functions.

In the window for , locate the section.

From the list, choose .

Modify the plot group to be compared with the figures in the reference paper: Change the generic labels to be more specific to this model. Normalize the electron density with the ionized dopant concentration Nd.

Compare V and n with Thomas-Fermi Approximation

In the window, right-click and choose four times.

In the window for , type Compare V and n with Thomas-Fermi Approximation in the text field.

Click to expand the section. In the text area, type Compare V and n with Thomas-Fermi Approximation.

Locate the section. In the text field, type Electron density / Nd.

In the text field, type Electron potential energy (eV).

Electron Potential Energy

In the window, expand the node, then click .

In the window for , type Electron Potential Energy in the text field.

Locate the section.

Select the checkbox. In the associated text field, type Electron Potential Energy.

Electron Potential Energy from Previous Iteration

In the window, under > click .

In the window for , type Electron Potential Energy from Previous Iteration in the text field.

Locate the section.

Select the checkbox. In the associated text field, type Electron Potential Energy, Previous Iteration.

Electron Density from Weighted Sum

In the window, under > click .

In the window for , type Electron Density from Weighted Sum in the text field.

Locate the section. In the text field, type schrp1.n_sum/Nd.

Select the checkbox. In the associated text field, type Electron Density, Weighted Sum.

Electron Density

In the window, under > click .

In the window for , type Electron Density in the text field.

Locate the section. In the text field, type schrp1.rhoq/schrp1.q/Nd.

Select the checkbox. In the associated text field, type Electron Density.

Electron Density from Previous Iteration

In the window, under > click .

In the window for , type Electron Density from Previous Iteration in the text field.

Locate the section. In the text field, type schrp1.rhoq_old/schrp1.q/Nd.

Select the checkbox. In the associated text field, type Electron Density, Previous Iteration.

Add the electron density and potential energy profile from the Thomas–Fermi approximation for comparison. Since two y-axes are used, add an arrow to indicate that the potential energy profile uses the vertical scale on the right hand side.

Electron Density from Thomas-Fermi Approximation

Right-click and choose .

In the window for , type Electron Density from Thomas-Fermi Approximation in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type n_TF/Nd.

In the text field, type Electron Density, Thomas-Fermi.

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

From the list, choose .

Electron Potential Energy from Thomas-Fermi Approximation

Right-click and choose .

In the window for , type Electron Potential Energy from Thomas-Fermi Approximation in the text field.

Locate the section. Select the checkbox.

Locate the section. In the text field, type -V.

In the text field, type Electron Potential Energy, Thomas-Fermi.

Annotation 1

In the window, right-click and choose .

In the window for , locate the section.

Select the checkbox.

Locate the section. In the text field, type $\longrightarrow$.

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

In the text field, type 0.47.

Locate the section. Clear the checkbox.

Locate the section. Select the checkbox.

Add the plot to be plotted while solving, then click compute and watch the progress of the iteration in the plot and in the table of absolute errors.

Study 2: Schrödinger-Poisson

Solver Configurations

In the window, expand the > node.

Solution 2 (sol2)

In the window, expand the > > > node, then click .

In the window for , click to expand the section.

Select the checkbox.

In the table, enter the following settings:


Plot group	Plot window
Compare V and n with Thomas-Fermi Approximation	Graphics

In the toolbar, click

Results

Compare V and n with Thomas-Fermi Approximation

Add a zoomed-in version of the plot to compare with Fig. 2 in the reference paper. Move the arrow indicator to a proper position.

Compare V and n with Thomas-Fermi Approximation 1

Right-click and choose .

Right-click and choose four times.

In the window for , locate the section.

Select the checkbox.

In the text field, type 0.

In the text field, type 16.

In the text field, type 0.9.

In the text field, type 1.002.

In the text field, type -86.5e-3.

In the text field, type -80.5e-3.

Annotation 1

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type 14.5.

In the text field, type -0.0827.

In the toolbar, click

Optionally, add the contribution to the electron density from each azimuthal quantum number to the original plot, to compare with Fig. 3 in the reference paper. The expression for the plot is slightly complicated but is only needed for visualizing the subcontributions. Use the solution from the node in the solver sequence to provide data for the statistically weighted partial sums. The sum operator is used to do the summation and the withsol operator is used to select the solutions with the desired quantum number m to be summed over the eigenvalue lambda. Detailed discussions on the withsol operator can be found in the COMSOL Multiphysics Reference Manual.

Electron Density from m=0

In the window, right-click and choose .

In the window for , type Electron Density from m=0 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,0),setind(lambda,jj)),jj,1,4)/Nd.

In the text field, type m=0.

Locate the section. Find the subsection. From the list, choose .

From the list, choose .

Electron Density from m=+/-1

Right-click and choose .

In the window for , type Electron Density from m=+/-1 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,1),setind(lambda,jj)),jj,1,4)/Nd.

In the text field, type |m|=1.

Electron Density from m=+/-2

Right-click and choose .

In the window for , type Electron Density from m=+/-2 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,2),setind(lambda,jj)),jj,1,4)/Nd.

In the text field, type |m|=2.

Electron Density from m=+/-3

Right-click and choose .

In the window for , type Electron Density from m=+/-3 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,3),setind(lambda,jj)),jj,1,3)/Nd.

In the text field, type |m|=3.

Electron Density from m=+/-4

Right-click and choose .

In the window for , type Electron Density from m=+/-4 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,4),setind(lambda,jj)),jj,1,3)/Nd.

In the text field, type |m|=4.

Electron Density from m=+/-5

Right-click and choose .

In the window for , type Electron Density from m=+/-5 in the text field.

Locate the section. In the text field, type sum(withsol('sol4',schrp1.ni,setval(m,5),setind(lambda,jj)),jj,1,2)/Nd.

In the text field, type |m|=5.

In the toolbar, click

Finally check the charge neutrality as discussed in the reference paper. The built-in integration operator schrp1.int is used to compute the number of electrons per unit length inside the nanowire. This number is added to the contribution from the surface charge density given by the electric displacement field at the nanowire surface. The sum is then divided by the total number of ionized dopants per unit length. The ratio is subtracted by 1 to obtain the relative error.

Point Evaluation 1: check charge neutrality

In the toolbar, click

In the window for , type Point Evaluation 1: check charge neutrality in the text field.

Locate the section. From the list, choose .

Select Boundary 2 only.

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


Expression	Unit	Description
(schrp1.int(2piRschrp1.n_sum)-2piRes.nD/e_const)/(piR^2Nd)-1	1	Relative error

Click

We see that the relative error for charge neutrality is less than 1e-6.