Quantum matrix diagonalization visualized

TL;DR

This paper presents a visualization method for the diagonalization process of Hamiltonian matrices in quantum systems, using graphical tools and interactive software to enhance understanding.

Contribution

It introduces a visual approach to quantum matrix diagonalization, including Mathematica code and a JavaScript web app for interactive learning.

Findings

Graphical visualization of Hamiltonian diagonalization process

Interactive web app for step-by-step diagonalization

Enhanced understanding of quantum eigenvalue problems

Abstract

We show how to visualize the process of diagonalizing the Hamiltonian matrix to find the energy eigenvalues and eigenvectors of a generic one-dimensional quantum system. Starting in the familiar sine-wave basis of an embedding infinite square well, we display the Hamiltonian matrix graphically with the basis functions alongside. Each step in the diagonalization process consists of selecting a nonzero off-diagonal matrix element, then rotating the two corresponding basis vectors in their own subspace until this element is zero. We provide Mathematica code to display the effects of these rotations on both the matrix and the basis functions. As an electronic supplement we also provide a JavaScript web app to interactively carry out this process.

Figures

Equations

$$\hat{H}{\psi }_n=E_n{\psi }_n,$$

$$\hat{H}=-\frac{{\hslash }^2}{2M}\frac{d^2}{d{x}^2}+V(x).$$

$$\varphi_{m}(x)=\sqrt{\frac{2}{a}}\sin\Bigl{(}\frac{m\pi x}{a}\Bigr{)},\qquad m=1,2,\ldots,$$

$${\psi }_n=\sum _m{c}_{mn}{\phi }_m,$$

$$\sum _m⟨{\phi }_l|\hat{H}|{\phi }_m⟩c_{mn}=E_n\sum _m{c}_{mn}⟨{\phi }_l|{\phi }_m⟩=E_n{c}_{ln}.$$

$$H_{lm}=⟨{\phi }_l|\hat{H}|{\phi }_m⟩={\int }_0^a{\phi }_l(x)\hat{H}{\phi }_m(x)dx,$$

$$\sum _{m=1}^{\infty }{H}_{lm}{c}_{mn}=E_n{c}_{ln},$$

$$\left[\begin{array}{ccc}{H}_{11}& H_{12}& \cdots \\ H_{21}& H_{22}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{c}{c}_{1n}\\ c_{2n}\\ ⋮\end{array}\right]=E_n\left[\begin{array}{c}{c}_{1n}\\ c_{2n}\\ ⋮\end{array}\right].$$

$$C=\left[\begin{array}{ccc}{c}_{11}& c_{12}& \cdots \\ c_{21}& c_{22}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right].$$

$$C^{T}HC=\left[\begin{array}{ccc}{E}_1& 0& \cdots \\ 0& E_2& \cdots \\ ⋮& ⋮& \ddots \end{array}\right].$$

$$R=\left[\begin{array}{ccccccc}1& & & & & & \\ & \ddots & & & & & \\ & & \cos \theta & \cdots & -\sin \theta & & \\ & & ⋮& 1& ⋮& & \\ & & \sin \theta & \cdots & \cos \theta & & \\ & & & & & \ddots & \\ & & & & & & 1\end{array}\right],$$

$$H_{\text{new}}=R^T{H}_{\text{old}}R,$$

$$C=R_1{R}_2{R}_3\cdots ,$$

$$\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{cc}{H}_{mm}& H_{mn}\\ H_{nm}& H_{nn}\end{array}\right]\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right],$$

$$\theta=\frac{1}{2}\tan^{-1}\biggl{(}\frac{2H_{mn}}{H_{mm}-H_{nn}}\biggr{)}.$$

$$H_{mn}=\frac{1}{2}{n}^2{\pi }^2{\delta }_{nm}+⟨{\phi }_m|\hat{V}|{\phi }_n⟩.$$

Full text

I Introduction

A number of recent papers in educational physics journals have highlighted the usefulness of matrix methods for solving the time-independent Schrödinger equation. Applications have included one-dimensional potential wells,Marsiglio double wells,Jelic ; Dauphinee spherical potentials,Jugdutt periodic potentials,Pavelich ; LeVot ; Pavelich2 the helium atom,Hutchinson ; MasseWalker and other systems of two or three interacting particles.Na

Of course matrix methods date back to the birth of quantum mechanics. The new development is the ever-increasing availability of personal computers equipped with easy-to-use softwareMathematica ; Maple ; Matlab ; SciPy that can diagonalize large matrices almost instantly. Anyone who understands basic linear algebra can now use this software to solve quantum systems that could previously be tackled only by dedicated researchers.

We typically treat these miraculous diagonalization routines as black boxes, but there are good pedagogical reasons to peek inside and see how they work—or at least how they could work. One reason is that students are naturally curious about how things work. Most students of quantum mechanics have taken only a first course in linear algebra and have therefore learned only one matrix diagonalization algorithm: solve the characteristic polynomial for the eigenvalues, then plug these back into the eigenvalue equation and solve the resulting linear systems for the eigenvectors, one by one. But this process is practical only for very small matrices. How, then, do computers diagonalize matrices with hundreds, or even thousands, of rows and columns?

Besides satisfying students’ curiosity, looking inside a practical matrix diagonalization algorithm can help build students’ geometrical intuition for the high-dimensional vector spaces in which quantum states live. We express these states, and the operators that act on them, with respect to a particular basis. When we diagonalize a Hamiltonian matrix we are rotating our basis vectors. In a two-dimensional vector space this rotation is confined to a single plane, but in larger vector spaces the number of independent rotation planes grows roughly in proportion to the dimension squared.

In this paper we describe a diagonalization algorithm (originally due to Jacobijacobi ) that consists of successive rotations of basis vectors. We then show how to implement the algorithm in a visual way, allowing students to see the effects of each individual rotational step on both the Hamiltonian matrix and the basis functions. We provide Mathematica code for this implementation in Fig. 1, and we provide a JavaScript web app version as an electronic supplement.supplement Our hope is that this software will make the connection between wavefunctions and matrix representations more vivid for students, and help them understand the diagonalization process not as a mere symbol-pushing calculation but as a geometrical transformation.

II Putting the problem into matrix form

Consider a quantum particle in one dimension, trapped in a potential $V(x)$ of arbitrary shape such that its low-lying energy eigenstates are localized between $x=0$ and $x=a$. To find these energy eigenstates $\psi_{n}(x)$ and their corresponding energies $E_{n}$ we must solve the time-independent Schrödinger equation,

[TABLE]

where $\hat{H}$ is the Hamiltonian operator,

[TABLE]

We proceed by introducing a set of normalized basis functions,

[TABLE]

conveniently chosen to be eigenstates of the kinetic energy term of the Hamiltonian—that is, eigenstates of a hypothetical infinite square well extending from 0 to $a$. We can then expand each low-lying eigenstate in terms of these basis functions:

[TABLE]

where the first subscript on the coefficient $c_{mn}$ indicates which component, while the second subscript indicates which eigenstate. Inserting this expansion on both sides of Eq. (1) and then taking the inner product with an arbitrary basis function $\varphi_{l}$, we obtain

[TABLE]

The inner product on the left defines the $lm$ matrix element of $\hat{H}$,

[TABLE]

so Eq. (5) is simply

[TABLE]

which in matrix notation reads

[TABLE]

Solving this matrix eigenvalue equation is equivalent to solving the original operator eigenvalue equation, Eq. (1). Because the basis functions $\varphi_{m}(x)$ are real, the Hamiltonian matrix is real and symmetric. In numerical solutions one must truncate the matrix to a finite number of rows and columns, but this truncation is rarely a practical hindrance to solving one-dimensional problems. The most computationally intensive part of the calculation is not solving the eigensystem but rather computing the matrix elements in the first place, using Eq. (6).remark1 The method outlined above can also be readily generalized to other sets of basis functions and to multidimensional problems.Pavelich ; Pavelich2 ; Hutchinson ; MasseWalker ; Na

Without loss of generality we can choose the eigenvectors of Eq. (8) to be purely real and normalized. Then these column vectors form an orthogonal matrix:

[TABLE]

Performing a similarity transformation with $C$ diagonalizes the Hamiltonian:

[TABLE]

III Jacobi rotations

We now turn to the problem of solving the eigensystem of Eq. (8). Our method,NR originally due to Jacobi,jacobi consists of performing a succession of rotations in two-dimensional subspaces of the (truncated) high-dimensional vector space, with each rotation chosen to zero out a single off-diagonal element of $H$.

We accomplish each rotation with a matrix of the form

[TABLE]

consisting of an identity matrix except for the four entries that accomplish a rotation by $\theta$ in a particular plane; let us say these entries are in rows (and columns) $m$ and $n$. Under this rotation, the Hamiltonian transforms according to

[TABLE]

and we will choose the angle $\theta$ so that the $mn$ element of $H_{\textrm{new}}$ is zero. We then repeat this process with different choices of $m$ and $n$, zeroing out the off-diagonal elements of $H$ in succession. Subsequent rotations will again make $H_{mn}$ nonzero, because each transformation changes all the elements in rows (and columns) $m$ and $n$. Still, the off-diagonal elements of $H$ get smaller and smaller as we repeatedly cycle through all of the possible $mn$ pairs. (See Ref. NR, for a proof of convergence.) We continue until $H$ is diagonal to our desired precision.

Once we are finished, the diagonal elements of $H$ will be the eigenvalues we seek. Meanwhile, the product of all the rotation matrices will be (approximately) equal to the matrix $C$ of Eqs. (9) and (10):

[TABLE]

where the subscripts indicate the sequence of the rotations. We can then use Eq. (4) to construct the eigenfunctions. In fact, we can (and will) use successive approximations to $C$ to obtain approximate eigenfunctions, with ever-increasing accuracy, at each stage during the process.

To find the correct angle $\theta$ to zero out the new value of $H_{mn}$, it suffices to write a simplified version of the right-hand side of Eq. (12) for just the rotation’s two-dimensional subspace:

[TABLE]

where $H_{nm}=H_{mn}$ because $H$ is symmetric. Working out either off-diagonal element of this product and setting the result to zero gives

[TABLE]

The order in which we choose $mn$ pairs for Jacobi rotations is not critical. For “manual” diagonalization as we describe below, it is natural to begin with the larger off-diagonal elements of $H$. Alternatively, one can simply iterate through the off-diagonal elements in order by row and column. Both the Mathematica code that we present below and the JavaScript web app that we provide as a supplement leave the order up to the user.

IV Implementation in Mathematica

Figure 1 shows an implementation of the Jacobi rotation algorithm in Mathematica. The code uses natural units in which $\hbar=M=a=1$.

The first line of code sets the number of basis functions used in the expansion (Eq. (4)), which is then the dimension of every vector and matrix. The value 8 is large enough to get the main ideas across without requiring more than a few dozen rotational steps; for increased accuracy we have used values up to around 20.

The next two lines define the sinusoidal basis functions and the potential energy function. For demonstration purposes we use a harmonic oscillator potential, centered at $x=0.5$, with a classical angular frequency of 100 in natural units (sufficient to confine the low-lying energy eigenfunctions within the interval from 0 to 1).

Line 4 generates the Hamiltonian matrix by breaking up the inner product of Eq. (6) into kinetic and potential energy terms:

[TABLE]

The first term is just the eigenvalues of the embedding infinite square well, whose eigenfunctions are $\varphi_{n}(x)$. The Method and AccuracyGoal options speed up the numerical integration somewhat, and eliminate warning messages about matrix elements that are zero due to symmetry.

Line 5 initializes the matrix $C$ of Eq. (9) to the identity matrix, so we can later multiply it by each rotation matrix according to Eq. (13).

Lines 6–11 define a function that performs a single Jacobi rotation to zero out a selected off-diagonal cell $H_{mn}$. This function first finds the rotation angle using Eq. (15), then sets up the rotation matrix in the form of Eq. (11). Line 9 rotates the Hamiltonian matrix according to Eq. (12), while line 10 updates the $C$ matrix according to Eq. (13). For the matrix products we simply use Mathematica’s . operator, making no attempt to optimize for the sparse rotation matrix. (A much more efficient implementation is given in Ref. NR, .)

Lines 12–19 create our visualization scheme: a matrix plot of the Hamiltonian alongside a plot of the current basis functions overlaid on the potential (see Fig. 2). The Dynamic function causes these plots to update automatically each time rotateCell is executed. Lines 12 and 13 define two variables that should be adjusted from time to time as needed: a contrast for the matrix plot and a maximum energy for the potential/eigenfunction plot. Line 15 scales the wavefunctions by an arbitrary factor and shifts each of them upward by its (average) energy. Line 17 plots the Hamiltonian using a color scheme that maps zero to white, with distinct colors (modify according to taste!) for positive and negative values.

After executing lines 1 through 19, one merely has to repeatedly execute line 20, calling the rotateCell function for each Jacobi rotation. The function’s two arguments are the row and column numbers of the off-diagonal element (actually two symmetrically located elements) of $H$ that will be zeroed out. The matrix display will update to show the newly zeroed elements, along with side effects on the other elements in the same rows and columns. Meanwhile, the plots of the two corresponding basis functions will update to reflect their newly mixed shapes as well as their new average energies. Figure 2 shows the graphical display at three different stages during the diagonalization procedure.

The harmonic oscillator potential is a useful first example because we can compare to the known exact solutions. The algorithm can handle arbitrary potential shapes. Some interesting examples are explored in the exercises in Sec. VI.

After calling rotateCell two or three times for each nonzero off-diagonal element, it may be necessary to increase the contrast variable in order to distinguish the colors of the now-small off-diagonal elements from white. Alternatively, it may be more convenient to simply display the matrix numerically, with an instruction such as

Dynamic[MatrixForm[Chop[hMatrix]]]

(where the Chop function sets very small values to zero to simplify the display).

As a test, one can check that the instruction

Eigensystem[hMatrix]

yields the same eigenvalues and eigenvectors as the Jacobi rotation process. Going through so many steps to accomplish something Mathematica can do in a single line might seem foolish, but we feel that the visual step-by-step process provides insight that one cannot obtain from a black-box function call.

V Implementation as a web app

As an electronic supplementsupplement we provide an alternate visual implementation of the Jacobi rotation algorithm in the form of a JavaScript web app. This version has the drawback of making the code less accessible and harder to modify. On the other hand, its graphical user interface is more intuitive and appealing. Instead of typing the row and column numbers of a chosen off-diagonal matrix element, the user simply clicks or taps that element in the visual display. The two corresponding basis functions are immediately highlighted. Then the user presses and turns a graphical dial control by any desired angle, as both the matrix display and the basis functions update in real time to show the effects of the basis vector rotation. To facilitate diagonalizing larger matrices, the interface also provides shortcut buttons to select the largest remaining off-diagonal element and to rotate by the correct angle to zero out the selected element.

The web app provides a menu of several built-in potential functions, with adjustable parameters. There is also an option to draw an arbitrary potential function.

VI Student exercises

We hope that students will be curious enough to use the code we have provided for open-ended exploration, asking and answering their own questions. In case further direction might be helpful, here are some suggested exercises.

Derive Eq. (15) for the Jacobi rotation angle. 2. 2.

Explain the checkerboard appearance of the harmonic oscillator Hamiltonian matrix in Fig. 2. 3. 3.

Explain why the initial value of $H_{13}$ for the harmonic oscillator Hamiltonian matrix is positive. Suppose that you choose $H_{13}$ as the first matrix element to zero out via a Jacobi rotation. What is the rotation angle in degrees? What is the rotation matrix? Explain why the rotation causes the average energy of the $n=1$ basis state to decrease. 4. 4.

When you perform a Jacobi rotation on basis vectors $m$ and $n$, all the matrix elements in the corresponding rows and columns can (potentially) change. Why, then, don’t all the other basis functions change? Explain carefully. 5. 5.

Suppose that you’re in a hurry to find the ground-state energy and wavefunction, but you don’t need to know anything about the excited states. How should you go about choosing $mn$ pairs for Jacobi rotations? What if instead of the ground state you want to find only a particular excited state? 6. 6.

For the harmonic oscillator example described above, compare the energy eigenvalues to the exact values that you would obtain from an analytic treatment of the quantum harmonic oscillator. What is the main reason for the inaccuracies? What change should you make to obtain greater accuracy? 7. 7.

Suppose that the particle in question is an electron and that $a=1$ nm. What, then, is the size of one natural unit of energy, in electron-volts? (Hint: What combination of $\hbar$, $M$, and $a$ has units of energy?) 8. 8.

Modify the code in Fig. 1 to solve the “quantum bouncer” potential: $V(x)=\infty$ for $x<0$ and $V(x)=\alpha x$ for $x>0$, where $\alpha$ is a constant for which a reasonable value, in natural units, is 500. Explain why the initial Hamiltonian matrix for this system is qualitatively very different from that for the harmonic oscillator. Note the final energy eigenvalues when nMax = 8, then check them by increasing nMax and either repeating the Jacobi rotation process or simply using the Eigensystem function. Briefly discuss the resulting energy level structure and the shapes of the energy eigenfunctions. 9. 9.

Modify the code in Fig. 1 to solve a symmetric double-well potential. A simple way to do this is to use Mathematica’s If function to set $V(x)$ to a positive constant within a narrow interval around $x=0.5$, and set $V(x)=0$ everywhere else. Alternatively, you could devise a smooth double-well potential using a quartic polynomial. Either way, make the central barrier high enough that several low-lying energy eigenvalues are below the barrier height. Explain the general features of the initial Hamiltonian matrix, then carry out the solution and discuss the resulting eigenvalues and eigenfunctions. 10. 10.

Sometimes a Jacobi rotation will put the basis states into the “wrong” order, with a higher-energy state having a lower $n$ value than a lower-energy state. (See Fig. 2 for an example.) Write a Mathematica function, similar in many ways to rotateCell, that swaps the order of any two basis states, so that you can use this function at any time to put the states back into the “right” order. 11. 11.

Write a Mathematica function to automatically carry out an arbitrary number of successive Jacobi rotations, either zeroing out the largest remaining off-diagonal Hamiltonian matrix element during each step, or simply cycling through all the off-diagonal elements in order by row and column. (It’s a nice touch to use the Pause function to slow down the process, so you can watch the display update with each step.) After testing your code, use it to explore the rate at which the Jacobi rotation algorithm converges. That is, how many rotations per off-diagonal element are needed before all of these elements are smaller than a given threshold?