Simulating a Silicon Quantum Dot in a Uniform Magnetic Field

As of version 5.6 of the COMSOL Multiphysics® software, in the Semiconductor Module, we have introduced the capability of handling multicomponent wave functions to the Schrödinger Equation physics interface. Following the discussion in a previous blog post, we continue to explore this functionality using a model of a silicon quantum dot in a uniform magnetic field as an example.

Introduction to Quantum Dots

Quantum dots are an essential ingredient in nanotechnology with potential applications in solar cells, LEDs, displays, photodetectors, and quantum computing. One of the papers on the last application area is published by Jock et al. (Ref. 1) on the topic of spin-orbit qubits. In the Supplementary Note 1 of the paper, the authors provided the formulation describing a silicon quantum dot in a uniform magnetic field and showed the numerical solution in Supplementary Figure 1, which will be reproduced by this model.

The Formulation

Equation 1 in the Supplementary Note 1 of Ref. 1 gives the single-electron Hamiltonian for a silicon quantum dot in a uniform magnetic field \mathbf{B}, without including spin-orbit coupling:

(1)

where m_\perp and m_\parallel are the effective masses in the lateral and vertical directions, respectively; V is the confinement potential energy of the quantum dot; \mu_B is the Bohr magneton; \mathbf{\sigma} is the vector of Pauli matrices; the gyromagnetic ratio tensor is assumed to be a scalar of value 2 following the paper; and the kinetic momentum \mathbf{P} is given by

(2)

where e is the elementary charge, \mathbf{A} is the vector potential given by \mathbf{A}(\mathbf{r})=\frac{1}{2}\mathbf{B}\times\mathbf{r}, and there is no minus sign in front of the imaginary unit i because all physics interfaces in COMSOL Multiphysics take the engineering sign convention of exp(-i k x + i \omega t) instead of the physics sign convention of exp(i k x – i \omega t).

The confinement potential energy term V(\mathbf{r}) is given by Eq. 9 in the paper:

(3)

where the first two terms represent an anisotropic harmonic trapping potential in the lateral directions with trap angular frequencies \omega_x and \omega_y in the x– and y directions, respectively; the third term describes an electric field F_z in the z direction; q is the particle charge (q=-e for electrons); and the last term is a potential barrier of height U_0 at the silicon-oxide interface at z = 0.

Parameters and Variables

The modeling parameters are shown in the table below:

All parameters are given in the paper, except the magnetic field, for which a guess value of 1[T] is used.

The vector potential \mathbf{A}(\mathbf{r})=\frac{1}{2}\mathbf{B}\times\mathbf{r} is specified using variables, assuming the uniform magnetic field points in the y direction:

Physics Settings

The Schrödinger Equation physics interface is set up with two wave function components, psiu and psid, for the spin-up and spin-down components, respectively:

The default Effective Mass domain condition is used to specify the lateral and vertical effective masses:

Even though, in this model, the two wave function components share the same effective masses, the user interface allows different values for different wave function components for general cases.

Three Electron Potential Energy domain conditions are used to set up the confinement potential energy term V(\mathbf{r}) (Eq. 3) for the lateral harmonic trapping:

For the electric field:

For the oxide barrier:

The contribution to the kinetic momentum from the vector potential (the second term in Eq. 2) is implemented using the Lorentz Force domain condition:

Finally, the magnetic coupling of the spin-up and spin-down components (the last term in Eq. 1) is implemented with the Zeroth Order Hamiltonian domain condition:

Note that the factor of \hbar^2/2m_e assumed by the domain condition should be compensated in the inputs, so that the net result is the last term in Eq. 1:

(4)

In the screenshot above for the Settings window of the Zeroth Order Hamiltonian domain condition, the first row of the Hamiltonian input table implements the (1,2) element of the matrix in Eq. 4, and the second row implements the (2,1) element.

A simple swept mesh is created and an eigenvalue study is used to find the first few eigenstates.

Understanding the Results

By plotting the real and imaginary parts of the spin-up wave function component of the first two eigenstates, we find that they have similar magnitudes and shapes, except for some sign flips:

Spin-Up Wave Function Component
Ground state, real part	Ground state, imaginary part
First excited state, real part	First excited state, imaginary part

The spin-down wave function components also have this trend among themselves and as compared to the spin-up components.

To understand this observation, we evaluate the wave function components near the peak density point:

We find that for the ground state (the first row of the result table shown above), up to an overall scaling factor of about 3.6e11, the amplitude of the spin-up component is 1+i=(1)(1+i), and the amplitude of the spin-down component is 1-i=(-i)(1+i). Therefore, the vector formed by the two components is proportional to \left(\begin{array}{c}1\\-i\end{array}\right), which is recognized as the spin-down eigenstate of the y-spin operator S_y. This is consistent with the intuitive picture that the lower-energy state of an electron in a magnetic field has its spin magnetic moment parallel to the magnetic field, and thus the spin is antiparallel to the magnetic field.

Similarly, for the first excited state, up to an overall scaling factor of about 3.6e11, the amplitude of the spin-up component is seen to be 1-i=(1)(1-i), and the amplitude of the spin-down component is 1+i=(+i)(1-i). Therefore, the vector formed by the two components is proportional to \left(\begin{array}{c}1\\+i\end{array}\right), which is recognized as the spin-up eigenstate of the y-spin operator S_y. This is consistent with the intuitive picture that the higher-energy state of an electron in a magnetic field has its spin magnetic moment antiparallel to the magnetic field, and thus the spin is parallel to the magnetic field.

This observation can be further confirmed by comparing the energy difference between the two computed eigenstates, with the expected energy difference between a spin-up and a spin-down electron in a uniform magnetic field (2\, \mu_B B), using global evaluations:

We see a very good agreement of the computed and expected energy difference between the first two eigenstates — both evaluate to the same value of 0.116 meV.

Probability Density and Kinetic Momentum Density

The plot below shows the probability density (blue-gray gradient) and kinetic momentum density (red arrows) of the ground state on the xz-plane. It compares well with Supplementary Figure 1 in the paper.

The plot below shows the probability density (isosurfaces) and kinetic momentum density (arrows) of the eighth excited state for the model thumbnail. The isosurfaces of the probability density are made partially transparent using the Transparency subnode, available as of COMSOL Multiphysics version 5.6.

Now It’s Your Turn…

In this post and the previous blog post on the k • p method for a strained wurtzite GaN band structure, we have discussed the functionalities of the Schrödinger Equation physics interface for multicomponent wave functions. We would love to hear how you use them for your simulation work.

Try modeling a silicon quantum dot in a uniform magnetic field by clicking the button below, which will take you to the Application Gallery, where you can download the MPH-file.

Reference

1. R.M. Jock et al., “A silicon metal-oxide-semiconductor electron spin-orbit qubit,” Nature Communications, 9:1768, 2018.