COMSOL 6.4 - Solving the Hydrogen Atom

In this tutorial model, we take a look at perhaps the most significant exact results of quantum mechanics — the ground state and the first few excited states of the hydrogen atom — and go over how to reproduce these results using COMSOL Multiphysics.

The hydrogen atom with its single electron is described by the time-independent Schrödinger equation (TISE):

(1)

Here,

, m, e, and ε0 are the reduced Planck constant, electron mass, elementary charge, and the permittivity of vacuum, respectively. The dependent variable to be solved for, ψ(r), is the wave function: a complex scalar field whose norm, | ψ(r) |2, gives the probability of finding the electron at r. The operator

, is the Hamiltonian, which corresponds to the total energy of the system. The Laplacian term gives the kinetic energy of the electron, whereas the potential energy is simply given by the Coulomb potential due to the positively charged nucleus at the origin. The TISE is in the form of an eigenvalue problem for the unknown total energy, E, and can be solved exactly by the usual tools of the trade: separation of variables and series expansions. The resulting eigenstates are labeled by three integers: the quantum numbers, n, l, and m, where n, l ∈ N, l ≤ n and m = −l, …, 0, …, l. Using spherical polar coordinates with θ and ϕ denoting the polar and azimuthal angles, respectively, they can be expressed as

(2)

where Rnl(r) is the radial wave functions and

, expressed in terms of the Legendre polynomials

, are known as the spherical harmonics. The eigenenergy only depends on the principal quantum number, n:

(3)

where RH ≈ 13.6057 eV is known as the Rydberg constant.

Model Definition

In the case of the hydrogen atom, the TISE is simply a PDE for a complex-valued scalar field in 3D, thus allowing for a FEM-based solution of the eigenvalue problem. The Semiconductor Module, an add-on to COMSOL Multiphysics, offers a Schrödinger Equation interface (abbreviated as schr), which we can use for building our hydrogen atom model.

Set up a sphere of radius A0 = 15a0, where a0 =

2/(mee2) ≈ 52.9 pm is the Bohr radius. In the Schrödinger Equation interface, modify the default nodes (Effective Mass and Electron Potential Energy) so that the effective mass is equal to me (this model does not use a periodic lattice) and the electron potential energy is equal to the Coulomb potential of the nucleus. The only boundary condition needed here is ψ → 0 as r → ∞, which corresponds to bound states. Therefore, surround the sphere with a thin layer and define it as an infinite element domain. To improve the computation speed, use a mesh that becomes progressively coarser in the radial direction. Finally, specify the energy scale as eV and estimate −15 eV as the starting point for the eigenvalue search.

Results and Discussion

Table 1 shows the eigenenergies obtained for the first two principal quantum numbers together with their theoretical values.

Table 1: eigenenergies obtained for the first two principal quantum numbers — analytical values vs. model simulation results.
n	Analytical method (RH/n2 (eV))	COMSOL (En (eV) from schr)
1	-13.6057	-13.6108
2	-3.4014	-3.4013

Note that the eigenvalues depend on the mesh used, so we always recommend running a mesh refinement study to get reliable results. Figure 1 shows the shape of the eigenstates, or orbitals, plotted using 3D isosurface plots.

Figure 1: Orbital shapes of unperturbed hydrogen, visualized using a series of displaced 3D isosurface plots.

Adopt the labeling used extensively in atomic physics: 1s, 2s, 2px, …. It is clear from the figure that the s-states are radially symmetric, whereas each of the p-states has cylindrical symmetry along one of three mutually perpendicular axes. For the spherically symmetric s-states, the radial probability density has the simple analytical form

(4)

Figure 2 and Figure 3 show a comparison between the analytical and simulated results.

Figure 2: Radial probability of the 1s state (ground state). Analytical results are shown using solid lines, while numerical result obtained in COMSOL Multiphysics are shown using point markers.

Figure 3: Radial probability of the 2s state (first excited state). Analytical results are shown using solid lines, while numerical result obtained in COMSOL Multiphysics are shown using point markers.

Stark Effect

In the hydrogen atom, the Stark effect can be observed as a possible shift of the energy levels, when an external electric field is applied.The perturbation can be represented by the Hamiltonian VStark = eεextz. To implement the Stark effect in COMSOL Multiphysics, simply add another Electron Potential Energy node to the model. Then, the split energy levels and the corresponding orbital shapes can readily be obtained. To make the effect clearly visible, use an extremely high external field of ∼ 2×108 V/m.

See the Modeling Instructions section for more detailed explanations on the model setup.

Semiconductor_Module/Quantum_Systems/solving_hydrogen_atom

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

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 solving_hydrogen_atom_parameters.txt.

Geometry 1

Sphere 1 (sph1)

In the toolbar, click

In the window for , locate the section.

In the text field, type A0+0.1*A0.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.1*A0

Point 1 (pt1)

In the toolbar, click

and choose .

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Psi100	exp(-sys2.r/a0_const)/(sqrt(pi*a0_const^3))		Ground state
Psi200	0.25(2-sys2.r/a0_const)/(sqrt(2pia0_const^3))exp(-sys2.r/a0_const/2)		n=2, l,m=0 eigenstate