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?