Numerical Diagonalization for Exciton Problem in Magnetic Fields

Hello! In our last two lessons, we explored the fascinating and complex intermediate magnetic field regime where $γ \approx 1$ . We saw that the non-integrable nature of the Hamiltonian gives rise to quantum chaos, which manifests as observable signatures like avoided crossings and specific level spacing statistics.

A critical question remains: how do we actually calculate the energy level diagrams that reveal this chaotic behavior? Perturbation theory, which served us well in the low-field limit, breaks down here. The answer lies in tackling the Schrödinger equation head-on with numerical methods.

This lesson is dedicated to exactly that. We will learn how to apply numerical diagonalization to solve the exciton problem in an arbitrary magnetic field strength. Given your background in computational physics and programming, this lesson will connect the abstract quantum mechanics to a concrete computational framework, a process central to modern physics research.

From Operator to Matrix: The Essence of Numerical Diagonalization

The time-independent Schrödinger equation, $H^∣ ψ ⟩ = E ∣ ψ ⟩$ , is an operator eigenvalue equation. The core idea of numerical diagonalization is to transform this into a matrix eigenvalue problem, which can be solved using standard linear algebra techniques.

The general procedure is:

Choose a basis: Select a set of known functions (a basis set) ${∣ ϕ_{i} ⟩}$ that can represent the unknown eigenstates $∣ ψ ⟩$ as a linear combination: $∣ ψ ⟩ = \sum_{i} c_{i} ∣ ϕ_{i} ⟩$ .
Construct the Hamiltonian Matrix: Represent the operator $H^$ as a matrix in this basis. The elements of this matrix are given by $H_{ij} = ⟨ ϕ_{i} ∣ H^∣ ϕ_{j} ⟩$ .
Diagonalize the Matrix: The Schrödinger equation becomes the matrix equation $Hc = E c$ , where $H$ is the matrix and $c$ is a column vector of the coefficients $c_{i}$ . Solving this is equivalent to finding a transformation that makes $H$ diagonal. The diagonal entries will be the energy eigenvalues $E$ , and the corresponding eigenvectors $c$ give us the composition of each eigenstate in our chosen basis.

This following video provides an excellent introduction to this concept, using the example of spin systems. The principles are universal and apply directly to our exciton problem.

Exact Diagonalization for spin systems

To understand the fundamentals of transforming the Schrödinger equation into a matrix diagonalization problem, watch this video from Dr Mitchell's physics channel.

Watch the specified sections to grasp these key points: 00:37 - 02:39: The equivalence between solving the Schrödinger equation and diagonalizing the Hamiltonian matrix. 14:32 - 18:12: How to construct the Hamiltonian matrix by calculating the matrix elements (H_{ij} = \langle \phi_i | \hat{H} | \phi_j angle). 20:47 - 29:23: The rigorous derivation showing that this diagonalization is a basis transformation to the eigenbasis, where the Hamiltonian matrix is diagonal.

Now, let's apply this powerful "recipe" to our specific, and much more complex, exciton problem.

Step 1: Choosing a Suitable Basis

The choice of basis is arguably the most critical step. An ideal basis should provide an accurate representation of the true wavefunctions with the fewest possible basis states, to keep the problem computationally manageable.

For the exciton, which is fundamentally a Coulomb-bound electron-hole pair, a natural choice for the radial part of the wavefunction is a basis built from functions that are solutions to a Coulomb-like problem. While hydrogenic wavefunctions might seem like an option, a more flexible and mathematically complete choice (without needing continuum states) is the set of Coulomb-Sturmian functions.

For the angular part, we use spherical harmonics coupled with the various spin and quasi-spin angular momenta involved. A complete basis state for an exciton in a semiconductor with a complex valence band might look like:

∣ Ψ_{basis} ⟩ = ∣ N, L; (I, S_{h}), J; F, S_{e}; F_{t}, M_{F_{t}} ⟩

where $N$ and $L$ are the radial and orbital quantum numbers from the Coulomb-Sturmian part, and the other quantum numbers represent the coupling of the hole's quasi-spin ( $I$ ), hole spin ( $S_{h}$ ), and electron spin ( $S_{e}$ ).

The following readings from two research papers show how this is done in practice for the exciton in Cu₂O. This gives you a direct look at the methods used in cutting-edge research.

Magnetoexcitons break antiunitary symmetries

This paper by Schweiner et al. provides a concise description of the basis set used to study magnetoexcitons and the subsequent setup of the eigenvalue problem.

Please read the section that begins just after Eq. (9) and ends before Fig. 1. Focus on the description of the basis functions in Eq. (10) and the final formulation of the problem as a generalized eigenvalue problem, Ac = EMc.

Theory of excitons in Cu₂O

For a more detailed and rigorous treatment, this second paper by Schweiner et al. dedicates an entire section to the choice of basis.

Read Section II.C, 'Choice of the basis set'. This section justifies the use of Coulomb-Sturmian functions (Eq. 22) and meticulously builds the complete basis state (Eq. 23) by considering all relevant angular momenta. It culminates in the same generalized eigenvalue problem (Eq. 25) we saw in the previous paper.

Step 2: Constructing the Hamiltonian Matrix

With a basis in hand, the next step is to calculate the matrix elements $H_{ij} = ⟨ ϕ_{i} ∣ H^∣ ϕ_{j} ⟩$ . The Hamiltonian $H^$ is the full operator for an exciton in a magnetic field, including:

Kinetic energy of the electron and hole.
The screened Coulomb interaction.
Terms arising from the magnetic field (diamagnetic and Zeeman terms).
Terms describing the complex valence band structure (Luttinger Hamiltonian).

The full expression for $H^$ is given by Equation (3) in the first paper you read (LINK). Calculating the matrix elements involves extensive angular momentum algebra (using tools like the Wigner-Eckart theorem, 3-j, 6-j, and 9-j symbols) and evaluating radial integrals. The paper LINK provides the explicit formulas for these matrix elements in its appendices (C and D). You don't need to study them in detail, but it's valuable to see the level of complexity involved in translating the theoretical operator into computable numbers.

Step 3: Dealing with Computational Cost - The "Curse of Dimensionality"

In theory, our basis set is infinite. In practice, we must truncate it by setting a maximum for the quantum numbers (e.g., $N + L \leq N_{max}$ ). Even with a modest $N_{max}$ , the number of basis states can be in the thousands or tens of thousands.

The computational cost of diagonalizing an $M \times M$ matrix scales as $O (M^{3})$ . This rapid growth, known as the "curse of dimensionality," means that direct diagonalization quickly becomes intractable.

Modern computational methods in physics part 1: Diagonalization

This video gives an intuitive feel for why the size of the quantum state space—and thus the matrix we need to handle—grows so rapidly.

Watch the section from 06:20 to 10:04. Although it discusses a lattice of spins, the principle is the same: the Hilbert space dimension grows exponentially, making brute-force computation very difficult for large systems or high-resolution basis sets.

Making the Problem Tractable: Block Diagonalization

How can we fight this exponential growth? The most powerful weapon is symmetry.

If our Hamiltonian $H^$ commutes with another operator $Q^$ (i.e., $[H^, Q^] = 0$ ), then the quantum number $q$ associated with $Q^$ is conserved. This has a profound implication: the Hamiltonian cannot connect states with different values of $q$ .

If q_{i} \neq = q_{j}, then ⟨ ϕ_{i} ∣ H^∣ ϕ_{j} ⟩ = 0

If we organize our basis states by their value of $q$ , the huge Hamiltonian matrix breaks down into a block diagonal structure. All non-zero elements are contained within smaller square blocks along the diagonal. Each block corresponds to a specific value of the conserved quantum number.

A schematic of a block diagonal matrix. Each colored block can be diagonalized independently, significantly reducing the computational effort compared to diagonalizing the entire matrix at once.

For the magneto-exciton, conserved quantities could include:

Parity ( $Π$ ): Conserved if no external electric field is present.
z-component of total angular momentum ( $M_{J}$ ): Conserved if the magnetic field is aligned with the z-axis.

By separating the basis states into sectors of definite $(M_{J}, Π)$ , we can diagonalize each small block independently, which is vastly more efficient than diagonalizing the one giant matrix.

Exact Diagonalization for spin systems

The concept of block diagonalization is a cornerstone of computational physics. This video provides a detailed and clear explanation.

This is a key part of the lesson. Please watch carefully: 32:29 - 45:32: This long segment thoroughly explains how conserved quantum numbers lead to a block diagonal structure. It shows how this simplifies the problem into diagonalizing smaller, independent sub-matrices. 45:44 - 55:19: Here, the concepts are applied to a concrete 2-spin example. This demonstrates the entire process: identifying the conserved quantity (total S_z), constructing the matrix, seeing the block structure, and then solving a small 2x2 sub-block.

Test your understanding!

Suppose you are calculating the energy spectrum of an exciton in a magnetic field $B = B z^$ with no electric field. The basis size for $N_{ma x} = 10$ is 10,000 states. However, you know that both parity ( $Π = \pm 1$ ) and the z-component of angular momentum $M_{J}$ are conserved.

How does this knowledge help you, and roughly what does the computational task look like now?

Show answer

Instead of diagonalizing one enormous 10,000 x 10,000 matrix, you can split the problem into many smaller ones. You would first group your 10,000 basis states by their parity ( $Π$ ) and their value of $M_{J}$ .

For example, you would construct a matrix just for the states with $(Π = +, M_{J} = 0)$ , another for $(Π = +, M_{J} = 1)$ , another for $(Π = -, M_{J} = 0)$ , and so on. Let's say there are 200 states in the $(Π = +, M_{J} = 0)$ sector. You would then only need to build and diagonalize a 200x200 matrix for this sector.

The total computational effort would be the sum of the costs for each small block (e.g., $C \times (20 0^{3} + 19 5^{3} + \dots)$ ), which is drastically less than the cost of the full problem ( $C \times 1000 0^{3}$ ). This makes the calculation feasible.

Conclusion

Today, we've demystified how the complex energy level diagrams in the intermediate magnetic field regime are actually computed. We've laid out the complete numerical recipe for what is often called "exact diagonalization"—a term that means finding the exact eigenvalues for a given truncated basis.

Key Takeaways:

Solving the Schrödinger equation for a non-integrable system like a magneto-exciton requires transforming it into a matrix eigenvalue problem.
The process involves three main steps: choosing a suitable basis (e.g., Coulomb-Sturmians), constructing the Hamiltonian matrix, and numerically diagonalizing it.
The computational cost is high due to the "curse of dimensionality," but can be massively reduced by exploiting symmetries to block diagonalize the Hamiltonian.
This method allows us to calculate the entire energy spectrum, providing the raw data—the energy levels—that we need to analyze for signatures of quantum chaos.

Preview of the Next Lesson

We now have the tool to generate the energy levels. In our last lesson, we learned what to look for in those levels (avoided crossings). The next logical step is to perform a quantitative analysis. In the next lesson, we will learn how to analyze level spacing statistics (e.g., Wigner-Dyson distribution) to quantitatively identify the transition to chaos. This will complete our study of the intermediate field regime by connecting the numerical output to the rigorous statistical tests for quantum chaos.

Recommended resources