Unexpected mixed singlet and triplet states
asked by Charles Cardot (2024/12/04 23:29)
Hello,
I was exploring how the Coulomb exchange coupling between different orbitals gives rise to states with different multiplicities. In this simplest non-trivial case (two unpaired electrons, one in each orbital) this should give rise to triplet and singlet states, where the triplet states have a <S^2> = 2 (s=1) and the singlet states have <S^2> = 0 (s=0). The code I used is shown below.
-- Create Index
Index, NFermi = CreateAtomicIndicesDict({"Metal_d","Metal_p","Metal_s"})
print(Index)
-------------------
-- Define Operators
-------------------
OppN_xy = NewOperator("Number", NFermi, {0,1}, {0,1}, {1,1})
OppN_yz = NewOperator("Number", NFermi, {2,3}, {2,3}, {1,1})
OppN_z2 = NewOperator("Number", NFermi, {4,5}, {4,5}, {1,1})
OppN_xz = NewOperator("Number", NFermi, {6,7}, {6,7}, {1,1})
OppN_x2y2 = NewOperator("Number", NFermi, {8,9}, {8,9}, {1,1})
OppN_x = NewOperator("Number", NFermi, {14,15}, {14,15}, {1,1})
OppN_y = NewOperator("Number", NFermi, {10,11}, {10,11}, {1,1})
OppN_z = NewOperator("Number", NFermi, {12,13}, {12,13}, {1,1})
OppN_d = NewOperator("Number", NFermi, Index['Metal_d'], Index['Metal_d'], {1,1,1,1,1,1,1,1,1,1})
OppN_p = NewOperator("Number", NFermi, Index['Metal_p'], Index['Metal_p'], {1,1,1,1,1,1})
OppN_s = NewOperator("Number", NFermi, Index['Metal_s'], Index['Metal_s'], {1,1})
YtoZtriple = YtoZMatrix({"Metal_d","Metal_p","Metal_s"})
OppSz_d = Rotate(NewOperator("Sz", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple)
OppSz_p = Rotate(NewOperator("Sz", NFermi, Index["Metal_p_Up"], Index["Metal_p_Dn"]), YtoZtriple)
OppSz_s = Rotate(NewOperator("Sz", NFermi, Index["Metal_s_Up"], Index["Metal_s_Dn"]), YtoZtriple)
OppLz_d = Rotate(NewOperator("Lz", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple)
OppSz_total = Rotate(NewOperator("Sz", NFermi, {1, 3, 5, 7, 9, 11, 13, 15, 17}, {0, 2, 4, 6, 8, 10, 12, 14, 16}), YtoZtriple)
OppSsqr_d = Rotate(NewOperator("Ssqr", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple)
OppSsqr_p = Rotate(NewOperator("Ssqr", NFermi, Index["Metal_p_Up"], Index["Metal_p_Dn"]), YtoZtriple)
OppSsqr_s = Rotate(NewOperator("Ssqr", NFermi, Index["Metal_s_Up"], Index["Metal_s_Dn"]), YtoZtriple)
OppSsqr_total = Rotate(NewOperator("Ssqr", NFermi, {1, 3, 5, 7, 9, 11, 13, 15, 17}, {0, 2, 4, 6, 8, 10, 12, 14, 16}), YtoZtriple)
OppLsqr_total = Rotate(NewOperator("Lsqr", NFermi, {1, 3, 5, 7, 9, 11, 13, 15, 17}, {0, 2, 4, 6, 8, 10, 12, 14, 16}), YtoZtriple)
OppJsqr_total = Rotate(NewOperator("Jsqr", NFermi, {1, 3, 5, 7, 9, 11, 13, 15, 17}, {0, 2, 4, 6, 8, 10, 12, 14, 16}), YtoZtriple)
print("\nCreate the Coulomb operator between the s- and d-shell\n")
OppUsdG2 = Rotate(NewOperator("U", NFermi, Index["Metal_s_Up"],Index["Metal_s_Dn"], Index["Metal_d_Up"],Index["Metal_d_Dn"], {0}, {1}), YtoZtriple)
print("Create the Coulomb operator within the d-shell\n")
OppUddF0 = Rotate(NewOperator("U", NFermi, Index["Metal_d_Up"],Index["Metal_d_Dn"], {1,0,0}), YtoZtriple)
OppUddF2 = Rotate(NewOperator("U", NFermi, Index["Metal_d_Up"],Index["Metal_d_Dn"], {0,1,0}), YtoZtriple)
OppUddF4 = Rotate(NewOperator("U", NFermi, Index["Metal_d_Up"],Index["Metal_d_Dn"], {0,0,1}), YtoZtriple)
print("Create the Coulomb operator between the p- and d-shell\n")
OppUpdF0 = Rotate(NewOperator("U", NFermi, Index["Metal_p_Up"],Index["Metal_p_Dn"], Index["Metal_d_Up"],Index["Metal_d_Dn"], {1,0}, {0,0}), YtoZtriple)
OppUpdF2 = Rotate(NewOperator("U", NFermi, Index["Metal_p_Up"],Index["Metal_p_Dn"], Index["Metal_d_Up"],Index["Metal_d_Dn"], {0,1}, {0,0}), YtoZtriple)
OppUpdG1 = Rotate(NewOperator("U", NFermi, Index["Metal_p_Up"],Index["Metal_p_Dn"], Index["Metal_d_Up"],Index["Metal_d_Dn"], {0,0}, {1,0}), YtoZtriple)
OppUpdG3 = Rotate(NewOperator("U", NFermi, Index["Metal_p_Up"],Index["Metal_p_Dn"], Index["Metal_d_Up"],Index["Metal_d_Dn"], {0,0}, {0,1}), YtoZtriple)
print("Create the Coulomb operator between the p- and s-shell\n")
OppUspF0 = Rotate(NewOperator("U", NFermi, Index["Metal_s_Up"],Index["Metal_s_Dn"], Index["Metal_p_Up"],Index["Metal_p_Dn"], {1}, {0}), YtoZtriple)
OppUspG1 = Rotate(NewOperator("U", NFermi, Index["Metal_s_Up"],Index["Metal_s_Dn"], Index["Metal_p_Up"],Index["Metal_p_Dn"], {0}, {1}), YtoZtriple)
---------------------
-- Build Hamiltonians
---------------------
nd = 1
scale = 10
G1pd = 2 * scale
G3pd = 4 * scale
F2pd = 2 * scale
F0pd = (1/15)*G1pd + (3/70)*G3pd
G2sd = 0.05
Bz = 0.001 * EnergyUnits.Tesla.value
GS_Hamiltonian = Bz * (2 * OppSz_d + OppLz_d) -- Magnetic term
-- s <-> d exchange term
GS_Hamiltonian = GS_Hamiltonian + G2sd * OppUsdG2
-- p <-> d exchange term
FS_Hamiltonian = GS_Hamiltonian + F0pd * OppUpdF0 + F2pd * OppUpdF2
FS_Hamiltonian = FS_Hamiltonian + G1pd * OppUpdG1 + G3pd * OppUpdG3
-----------------
-- Name operators
-----------------
GS_Hamiltonian.Name = "<E>"
FS_Hamiltonian.Name = "<E>"
OppN_d.Name = "N_d"
OppN_p.Name = "N_p"
OppN_s.Name = "N_s"
OppN_xy.Name = "<d_xy>"
OppN_yz.Name = "<d_yz>"
OppN_z2.Name = "<d_z2>"
OppN_xz.Name = "<d_xz>"
OppN_x2y2.Name = "<d_x2y2>"
OppN_x.Name = "<p_x>"
OppN_y.Name = "<p_y>"
OppN_z.Name = "<p_z>"
OppSz_d.Name = "<Sz_d>"
OppSz_p.Name = "<Sz_p>"
OppSz_s.Name = "<Sz_s>"
OppSz_total.Name = "<Sz_tot>"
OppSsqr_s.Name = "<S^2_s>"
OppSsqr_p.Name = "<S^2_p>"
OppSsqr_d.Name = "<S^2_d>"
OppSsqr_total.Name = "<S^2_tot>"
OppLsqr_total.Name = "<L^2_tot>"
OppJsqr_total.Name = "<J^2_tot>"
GS_oppList = {GS_Hamiltonian, OppSz_d, OppSz_s, OppSz_total, OppSsqr_total}
FS_oppList = {FS_Hamiltonian, OppSz_d, OppSz_p, OppSz_total, OppSsqr_total}
------------------------
-- Calculate Eigenvalues
------------------------
print("\nGround State Hamiltonian")
GS_Npsi=20
print("GS_Npsi", GS_Npsi)
StartRestrictions = {NFermi, 0, {"1111111111 000000 00", nd, nd}, {"0000000000 111111 00", 6, 6},{"0000000000 000000 11", 1, 1}}
GS_psiList = Eigensystem(GS_Hamiltonian, StartRestrictions, GS_Npsi, {{'Zero',1e-12}, {'Epsilon',1e-12}})
GS_psiList = Chop(GS_psiList)
print("\nGround State Hamiltonian Expectation Values")
PrintExpectationValues(GS_psiList, GS_oppList, {{"ColWidth", 10}})
print("\nFinal State Hamiltonian")
FS_Npsi = 60
print("FS_Npsi", FS_Npsi)
StartRestrictions = {NFermi, 0, {"1111111111 000000 00", nd, nd}, {"0000000000 111111 00", 5, 5},{"0000000000 000000 11", 2, 2}}
FS_psiList = Eigensystem(FS_Hamiltonian, StartRestrictions, FS_Npsi, {{'Zero',1e-12}, {'Epsilon',1e-12}})
FS_psiList = Chop(FS_psiList)
print("\nFinal State Hamiltonian Expectation Values")
PrintExpectationValues(FS_psiList, FS_oppList, {{"ColWidth", 10}})
print()
print("\nFinished")
Running this code produces the following output:
=============================================================
==== written by Maurits W. Haverkort, ====
==== Yi Lu, Robert Green, Sebastian Macke, ====
==== Marius Retegan, Martin Brass, and Simon Heinze ====
==== (cc) 1995 - 2022 code and binaries released ====
==== under creative commons CC-BY lisence ====
==== www.quanty.org ====
==== Beta version, be critical and report errors!!! ====
=============================================================
==== Version 0.7 beta ====
==== compiled at: May 16 2023 at 15:19:14 ====
=============================================================
==== When used in scientific publications please cite ====
==== one of the following papers as appropriate with ====
==== respect to the methods used in your publication: ====
==== Phys. Rev. B 85, 165113 (2012) ====
==== Phys. Rev. B 90, 085102 (2014) ====
==== Euro Phys. Lett. 108, 57004 (2014) ====
==== J. of Phys.: Conf. Series 712, 012001 (2016) ====
==== New Journal of Physics 22 (9), 093018 (2020) ====
=============================================================
Program executed on: Wed Dec 4 16:15:00 2024
Running on host : macchiato.phys.washington.edu
number of available processors : 48
maximum number of threads in parallel region: 48
Smallest positive float : 2.225074E-308
Smallest deviation from 1: 2.220446E-16
{ All = { 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 } ,
all = { 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 , 16 , 17 } ,
Metal_d = { 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 } ,
Metal_p_Up = { 11 , 13 , 15 } ,
Metal_p = { 10 , 11 , 12 , 13 , 14 , 15 } ,
Metal_s = { 16 , 17 } ,
Metal_p_Dn = { 10 , 12 , 14 } ,
Metal_s_Up = { 17 } ,
Metal_s_Dn = { 16 } ,
Metal_d_Up = { 1 , 3 , 5 , 7 , 9 } ,
Metal_d_Dn = { 0 , 2 , 4 , 6 , 8 } }
Create the Coulomb operator between the s- and d-shell
Create the Coulomb operator within the d-shell
Create the Coulomb operator between the p- and d-shell
Create the Coulomb operator between the p- and s-shell
Ground State Hamiltonian
GS_Npsi 20
Start of BlockGroundState. Converge 20 states to an energy with relative variance smaller than 1.000000000000000E-12
Start of BlockOperatorPsiSerialRestricted
Outer loop 1, Number of Determinants: 20 20 last variance 7.333925802638944E-05
Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 20 states to an energy with relative variance smaller than 1.000000000000000E-12
Start of BlockOperatorPsiSerial
Ground State Hamiltonian Expectation Values
<E> <Sz_d> <Sz_s> <Sz_tot> <S^2_tot>
1 -0.01 -0.5 -0.5 -1 2
2 -0.01 -2.894e-6 2.894e-6 0 2
3 -0.01 -0.5 -0.5 -1 2
4 -0.01 -2.894e-6 2.895e-6 1.35e-9 2
5 -0.01 0.5 0.5 1 2
6 -0.01 -0.5 -0.5 -1 2
7 -0.01 -2.9e-6 2.888e-6 -1.196e-8 2
8 -0.01 0.5 0.5 1 2
9 -0.01 -0.5 -0.5 -1 2
10 -0.01 -2.894e-6 2.894e-6 0 2
11 -0.01 -0.5 -0.5 -1 2
12 -0.01 0.5 0.5 1 2
13 -0.01 -2.894e-6 2.894e-6 0 2
14 -0.01 0.5 0.5 1 2
15 -0.01 0.5 0.5 1 2
16 0.01 2.894e-6 -2.894e-6 0 0
17 0.01 2.894e-6 -2.894e-6 0 0
18 0.01 2.894e-6 -2.894e-6 0 0
19 0.01 2.894e-6 -2.894e-6 0 0
20 0.01 2.894e-6 -2.894e-6 0 0
Final State Hamiltonian
FS_Npsi 60
Start of BlockGroundState. Converge 60 states to an energy with relative variance smaller than 1.000000000000000E-12
Start of BlockOperatorPsiSerialRestricted
Outer loop 1, Number of Determinants: 60 60 last variance 4.747935119629589E+01
Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 60 states to an energy with relative variance smaller than 1.000000000000000E-12
Start of BlockOperatorPsiSerial
Final State Hamiltonian Expectation Values
<E> <Sz_d> <Sz_p> <Sz_tot> <S^2_tot>
1 -7.05762 -0.5 -0.5 -1 2
2 -7.05762 -2.171e-9 2.171e-9 0 2
3 -7.05762 -0.5 -0.5 -1 2
4 -7.05762 0.5 0.5 1 2
5 -7.05762 -2.171e-9 2.171e-9 0 2
6 -7.05762 -0.5 -0.5 -1 2
7 -7.05762 0.5 0.5 1 2
8 -7.05762 -2.171e-9 2.171e-9 0 2
9 -7.05762 0.5 0.5 1 2
10 -4.20048 -0.5 -0.5 -1 2
11 -4.20048 -0.5 -0.5 -1 2
12 -4.20048 -3.94e-9 3.939e-9 0 2
13 -4.20048 -0.5 -0.5 -1 2
14 -4.20048 -3.939e-9 3.939e-9 0 2
15 -4.20048 0.5 0.5 1 2
16 -4.20048 -0.5 -0.5 -1 2
17 -4.20048 -3.94e-9 3.939e-9 0 2
18 -4.20048 0.32472 0.32472 0.649441 2
19 -4.20048 -0.32472 -0.32472 -0.649441 2
20 -4.20048 -3.939e-9 3.939e-9 0 2
21 -4.20048 0.5 0.5 1 2
22 -4.20048 -0.5 -0.5 -1 2
23 -4.20048 -3.94e-9 3.938e-9 0 2
24 -4.20048 -0.495622 -0.495622 -0.991244 2
25 -4.20048 0.495622 0.495622 0.991244 2
26 -4.20048 -3.939e-9 3.939e-9 0 2
27 -4.20048 0.5 0.5 1 2
28 -4.20048 -3.939e-9 3.939e-9 0 2
29 -4.20048 0.5 0.5 1 2
30 -4.20048 0.5 0.5 1 2
31 0.942381 -0.5 -0.5 -1 2
32 0.942381 -0.5 0.5 0 1
33 0.942381 -0.5 0.5 0 1
34 0.942381 -0.5 -0.5 -1 2
35 0.942381 -0.5 -0.5 -1 2
36 0.942381 -0.5 0.5 0 1
37 0.942381 0.5 0.5 1 2
38 0.942381 0.5 -0.5 0 1
39 0.942381 -0.5 0.5 0 1
40 0.942381 -0.5 -0.5 -1 2
41 0.942381 0.5 0.5 1 2
42 0.942381 0.5 -0.5 0 1
43 0.942381 -0.5 0.5 0 1
44 0.942381 -0.5 -0.5 -1 2
45 0.942381 0.5 -0.5 0 1
46 0.942381 0.5 0.5 1 2
47 0.942381 0.5 -0.5 0 1
48 0.942381 0.5 0.5 1 2
49 0.942381 0.5 -0.5 0 1
50 0.942381 0.5 0.5 1 2
51 10.4934 3.939e-9 -3.939e-9 0 0
52 10.4934 3.939e-9 -3.939e-9 0 0
53 10.4934 3.939e-9 -3.939e-9 0 0
54 10.4934 3.939e-9 -3.939e-9 0 0
55 10.4934 3.939e-9 -3.939e-9 0 0
56 10.4934 3.939e-9 -3.939e-9 0 0
57 10.4934 3.939e-9 -3.939e-9 0 0
58 19.609 2.171e-9 -2.171e-9 0 0
59 19.609 2.171e-9 -2.171e-9 0 0
60 19.609 2.171e-9 -2.171e-9 0 0
Finished
The Hamiltonian setup is a system with only a single *d* electron, and one hole in either the *s* or *p* level. When there is a hole in the *s* orbital there are 20 eigenstate (2×10) with 15 of them being triplet states and 5 of them being singlet states at a slightly different energy (G2sd = 0.05 eV). This is the expected behavior.
However, for the system with a hole in the *p* orbital there are 60 eigenstates(6×10) with some states which are clearly a mix of singlet and triplet (<S^2> = 1 or non integer value). My understanding is that the degeneracy between the singlet and triplet states is broken by Coulomb exchange, but the eigen states that are being found seem to no be obeying that.
I don't have any spin-orbit coupling so spin should be a good quantum number in this toy model. Does anyone have an idea what is giving rise to this behavior, specifically why I'm finding <S^2> values other than 0 and 2?
Answers
Dear Charles,
this is because for a $p^5 d^1$ configuration without spin-orbit coupling, the $^1P$ and $^3P$ are degenerate, and thus any linear combination of these states is a valid eigenstate.
There are a few tests you can do.
(1) If you look at the commutator of $S^2$ and $H$ that this should be zero. You can add the following line to the end of the code
Now you will find in your case that the commutator is of the order of $10^{-7}$. This is due to the magnetic field. If you remove this, or replace it by the operator $J_z$ acting on all shells, the commutator becomes zero.
(2) In order to split the degenerate state you could slightly modify the Hamiltonian in order to lift the degeneracy.
(3) You can add the operator $L^2$ to the printed expectation values. Be careful here, as your definition of $L^2$ is incorrect. You define the operator $L^2$ for a state with $l=4$ (multiplicity of 9), but you want the coupled angular momentum of an s,p, and d sub-shell. You can make the $L^2$ operator as follows
OppLx_d = Rotate(NewOperator("Lx", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple) OppLx_p = Rotate(NewOperator("Lx", NFermi, Index["Metal_p_Up"], Index["Metal_p_Dn"]), YtoZtriple) OppLx_s = Rotate(NewOperator("Lx", NFermi, Index["Metal_s_Up"], Index["Metal_s_Dn"]), YtoZtriple) OppLx = OppLx_d + OppLx_p + OppLx_s OppLy_d = Rotate(NewOperator("Ly", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple) OppLy_p = Rotate(NewOperator("Ly", NFermi, Index["Metal_p_Up"], Index["Metal_p_Dn"]), YtoZtriple) OppLy_s = Rotate(NewOperator("Ly", NFermi, Index["Metal_s_Up"], Index["Metal_s_Dn"]), YtoZtriple) OppLy = OppLy_d + OppLy_p + OppLy_s OppLz_d = Rotate(NewOperator("Lz", NFermi, Index["Metal_d_Up"], Index["Metal_d_Dn"]), YtoZtriple) OppLz_p = Rotate(NewOperator("Lz", NFermi, Index["Metal_p_Up"], Index["Metal_p_Dn"]), YtoZtriple) OppLz_s = Rotate(NewOperator("Lz", NFermi, Index["Metal_s_Up"], Index["Metal_s_Dn"]), YtoZtriple) OppLz = OppLz_d + OppLz_p + OppLz_s OppLsqr_total = OppLx * OppLx + OppLy * OppLy + OppLz * OppLzHope this helped, Maurits
Hello Maurits,
This was very helpful, thank you. I do not understand intuitively why those two multiplets are degenerate, but it's clear that you're correct. I tried your suggestion about testing the commutation of <S^2> with my Hamiltonian and indeed they commute. Adding the small term proportional to <S^2> also worked to lift the degeneracy, and now I can simply classify my final states as singlet or triplet.
Sincerely, Charles Cardot
Dear Charles,
An intuitive understanding of atomic multiplets would be nice… The book of Robert Cowan on “The theory of atomic structure and spectra” describes the multiplets of an atom with one unpaired valence electron and core hole in some detail. That said, there are so many instances I would normally take our numerical tools and just calculate the multiplets and look at expectation values. (i.e. L^2, S^2, and J^2 of each sub-shell and the total atom and or L_a . L_b between two shells if L,S, or J is not an (half) integer quantum number. ) You can do this whilst turning on or of the spin-orbit coupling, direct (F's) and exchange (G's) coulomb interaction.
Best wishes, Maurits