Note
Go to the end to download the full example code.
Exact diagonalization¶
Here we show how to find the eigenvalues and eigenvectors of a many-body Hamiltonian of fermions with Coulomb interactions. We then determine their spin and orbital angular momentum and how this changes when we switch on spin-orbit coupling.
Import the necessary modules.
import numpy as np
import matplotlib.pyplot as plt
import scipy
import edrixs
Parameters¶
Define the orbital angular momentum number
Coulomb interactions¶
The Coulomb interactions in EDRIXS are described by a tensor. Understanding this in full is complicated and requires careful consideration of the symmetry of the interactions. See example 6 for more discussion if desired. EDRIXS can construct the matrix via
Create basis¶
Now we build the binary form of the Fock basis
[[1 1 0 0 0 0]
[1 0 1 0 0 0]
[1 0 0 1 0 0]
[1 0 0 0 1 0]
[1 0 0 0 0 1]
[0 1 1 0 0 0]
[0 1 0 1 0 0]
[0 1 0 0 1 0]
[0 1 0 0 0 1]
[0 0 1 1 0 0]
[0 0 1 0 1 0]
[0 0 1 0 0 1]
[0 0 0 1 1 0]
[0 0 0 1 0 1]
[0 0 0 0 1 1]]
We expect the number of these states to be given by the mathematical combination of two electrons distributed among six states (three spin-orbitals with two spins per orbital).
We predict C(norb=6, noccu=2)=15 states and we got 15, which is reassuring!
Note that in more complicated problems with both valence and core electrons, the edrixs convention is to list the valence electrons first.
Transform interactions into Fock basis¶
edrixs works by initiailly creating a Hamiltonian matrix
generated as
We needed to specify n_fermion = 4
because the
edrixs.build_opers
function can also make two fermion terms.
Diagonalize the matrix¶
For a small problem such as this it is convenient to use the native
scipy diagonalization routine. This returns eigenvalues
e
and eignvectors v
where eigenvalue e[i]
corresponds
to eigenvector v[:,i]
.
15 eignvalues and 15 eigvenvectors 15 elements long.
Computing expectation values¶
To interpret the results, it is informative to compute the expectations values
related to the spin
We again transform these matrices to our Fock basis to build the operators
Recall that quantum mechanics forbids us from knowing all three Cartesian components of angular momentum at once, so we want to compute the squares of these operators i.e.
Remember that the eigenvalues of
We can determine the degeneracy of the eigenvalues numerically and print out the values as follows
e = np.round(e, decimals=6)
degeneracy = [sum(eval == e) for eval in e]
header = "{:<3s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}"
print(header.format("# ", "E ", "S(S+1)", "L(L+1)", "Degen."))
for i, eigenvalue in enumerate(e):
values_list = [i, eigenvalue, S2_val[i], L2_val[i], degeneracy[i]]
print("{:<3d}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:>3d}".format(*values_list))
# E S(S+1) L(L+1) Degen.
0 3.800 2.000 2.000 9
1 3.800 2.000 2.000 9
2 3.800 2.000 2.000 9
3 3.800 2.000 2.000 9
4 3.800 2.000 2.000 9
5 3.800 2.000 2.000 9
6 3.800 2.000 2.000 9
7 3.800 2.000 2.000 9
8 3.800 2.000 2.000 9
9 4.040 0.000 6.000 5
10 4.040 0.000 6.000 5
11 4.040 0.000 6.000 5
12 4.040 0.000 6.000 5
13 4.040 0.000 6.000 5
14 4.400 0.000 -0.000 1
We see
Energy level diagram¶
Let us show our findings graphically
fig, ax = plt.subplots()
for i, eigenvalue in enumerate(np.unique(e)):
art = ax.plot([0, 1], [eigenvalue, eigenvalue], '-', color='C{}'.format(i))
ind = np.where(eigenvalue == e)[0][0]
L = (-1 + np.sqrt(1 + 4*L2_val[ind]))/2
S = (-1 + np.sqrt(1 + 4*S2_val[ind]))/2
message = "L={:.0f}, S={:.0f} ({:.0f})"
ax.text(1, eigenvalue, message.format(L, S, degeneracy[ind]),
horizontalalignment='right',
verticalalignment='bottom',
color='C{}'.format(i))
ax.set_ylabel('Energy')
for loc in ['right', 'top', 'bottom']:
ax.spines[loc].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.set_xticks([])
plt.show()

We see Hund’s rules in action! Rule 1 says that the highest spin
Spin orbit coupling¶
For fun, we can see how this changes when we add spin orbit coupling (SOC). This is a two-fermion operator that we create, transform into the Fock basis and add to the prior Hamiltonian. To make things easy, let us make the SOC small so that the LS coupling approximation is valid and we can still track the states.
Then, we redo the diagonalization and print the results.
e2, v2 = scipy.linalg.eigh(H2)
e2 = np.round(e2, decimals=6)
degeneracy2 = [sum(eval == e2) for eval in e2]
print()
message = "With SOC\n {:<3s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}"
print(message.format("#", "E", "S(S+1)", "L(L+1)", "J(J+1)", "degen."))
J2_val_soc = edrixs.cb_op(J2, v2).diagonal().real
L2_val_soc = edrixs.cb_op(L2, v2).diagonal().real
S2_val_soc = edrixs.cb_op(S2, v2).diagonal().real
for i, eigenvalue in enumerate(e2):
values_list = [i, eigenvalue, S2_val_soc[i], L2_val_soc[i], J2_val_soc[i],
degeneracy2[i]]
print("{:<3d}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:8.3f}".format(*values_list))
With SOC
# E S(S+1) L(L+1) J(J+1) degen.
0 3.673 1.927 1.927 -0.000 1.000
1 3.750 2.000 2.000 2.000 3.000
2 3.750 2.000 2.000 2.000 3.000
3 3.750 2.000 2.000 2.000 3.000
4 3.827 1.802 2.396 6.000 5.000
5 3.827 1.802 2.396 6.000 5.000
6 3.827 1.802 2.396 6.000 5.000
7 3.827 1.802 2.396 6.000 5.000
8 3.827 1.802 2.396 6.000 5.000
9 4.063 0.198 5.604 6.000 5.000
10 4.063 0.198 5.604 6.000 5.000
11 4.063 0.198 5.604 6.000 5.000
12 4.063 0.198 5.604 6.000 5.000
13 4.063 0.198 5.604 6.000 5.000
14 4.427 0.073 0.073 -0.000 1.000
and we make an equivalent energy level diagram.
fig, ax = plt.subplots()
for i, eigenvalue in enumerate(np.unique(e2)):
art = ax.plot([0, 1], [eigenvalue, eigenvalue], '-', color='C{}'.format(i))
ind = np.where(eigenvalue == e2)[0][0]
J = (-1 + np.sqrt(1+4*J2_val_soc[ind]))/2
message = "J={:.0f} ({:.0f})"
ax.text(1, eigenvalue, message.format(J, degeneracy2[ind]),
horizontalalignment='right',
verticalalignment='bottom',
color='C{}'.format(i))
ax.set_ylabel('Energy')
for loc in ['right', 'top', 'bottom']:
ax.spines[loc].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.set_xticks([])
plt.show()

It is clear that we have split the
Total running time of the script: (0 minutes 0.138 seconds)