Ylm- and spin-resolved 4f photoemission amplitudes in Quanty

Dear Lukasz,

Yes, maybe, a bit, or better said, probably Quanty can deliver what you need. I would not suggest to calculate final states. There are just to many of them. We can however calculate the one particle electron removal Green's function as a matrix, which you then can use to do photoelectron diffraction.

The object to calculate would be

$$G_{\tau,\tau'}(\omega) = \left\langle \psi_0 | a^{\dagger}_{\tau} \frac{1}{\omega - (H -E_0) + \mathrm{i} \Gamma/2} a_{\tau'} | \psi_0 \right\rangle$$

This you then can use as an input for the directional and polarisation dependent photo-electron emission.

The code to calculate this object is

ExamplePESTensor.Quanty

-- Determine working directory (same as the directory of the input script)
directory = ""
 
-- Parameters
-- These values are not material realistic and just taken as a dumy set
-- For a real calculation you need to fit these to the experiment, possibly starting
-- from ab-initio values.
 
-- You are free to choose your units, as long as you are consistent. I would recommend eV for all of them.
 
U  =  6                                            -- spherical part of Coulomb interaction
F2 = 10                                            -- Quadrupole part of Coulomb interaction
F4 =  6                                            -- k=4 part of Coulomb interaction
F6 =  4                                            -- k=6 part of Coulomb interaction
F0 = U + (4/195)*F2 + (2/143)*F4 + (100/5577)*F6   -- Due to Pauli principle the monopole part and the spherical part differ
zeta = 0.1                                         -- 4f spin-orbit coupling strength
 
ne = 12                                            -- Number of 4f electrons
 
NPsi = 13                                          -- Number of low energy eigenstates calculated
                                                   -- You need all states that are thermally populated
                                                   -- The spectra starting from one of these states is not influenced by this number
 
B = {0,0,1 * EnergyUnits.Tesla.value }             -- Possible magnetic field
 
Akm={{0,0,0}}                                      -- Crystal field expanded on spherical Harmonics
                                                   -- The literature has no single most used convention
                                                   -- I like to pick the energies of the irreducible representations
                                                   -- in eV as the parameterisation
                                                   -- See https://www.quanty.org/physics_chemistry/point_groups
                                                   -- for examples.
 
--------------------------------------------------
-- -- -- here starts the real calculation -- -- --
--------------------------------------------------
 
-- define the dimension and meaning of the one particle basis
 
NF=14
NB=0
IndexDn_4f={0,2,4,6,8,10,12}
IndexUp_4f={1,3,5,7,9,11,13}
 
-- We define several operators
 
-- angular momentum operators
OppSx   =NewOperator("Sx"   ,NF, IndexUp_4f, IndexDn_4f)
OppSy   =NewOperator("Sy"   ,NF, IndexUp_4f, IndexDn_4f)
OppSz   =NewOperator("Sz"   ,NF, IndexUp_4f, IndexDn_4f)
OppSsqr =NewOperator("Ssqr" ,NF, IndexUp_4f, IndexDn_4f)
OppSplus=NewOperator("Splus",NF, IndexUp_4f, IndexDn_4f)
OppSmin =NewOperator("Smin" ,NF, IndexUp_4f, IndexDn_4f)
 
OppLx   =NewOperator("Lx"   ,NF, IndexUp_4f, IndexDn_4f)
OppLy   =NewOperator("Ly"   ,NF, IndexUp_4f, IndexDn_4f)
OppLz   =NewOperator("Lz"   ,NF, IndexUp_4f, IndexDn_4f)
OppLsqr =NewOperator("Lsqr" ,NF, IndexUp_4f, IndexDn_4f)
OppLplus=NewOperator("Lplus",NF, IndexUp_4f, IndexDn_4f)
OppLmin =NewOperator("Lmin" ,NF, IndexUp_4f, IndexDn_4f)
 
OppJx   =NewOperator("Jx"   ,NF, IndexUp_4f, IndexDn_4f)
OppJy   =NewOperator("Jy"   ,NF, IndexUp_4f, IndexDn_4f)
OppJz   =NewOperator("Jz"   ,NF, IndexUp_4f, IndexDn_4f)
OppJsqr =NewOperator("Jsqr" ,NF, IndexUp_4f, IndexDn_4f)
OppJplus=NewOperator("Jplus",NF, IndexUp_4f, IndexDn_4f)
OppJmin =NewOperator("Jmin" ,NF, IndexUp_4f, IndexDn_4f)
 
-- spin-orbit coupling
Oppldots=NewOperator("ldots",NF, IndexUp_4f, IndexDn_4f)
 
-- Number operator
OppNTable = {}
for i=1,7 do
  OppNTable[2*i-1]={IndexUp_4f[i],-IndexUp_4f[i],1}
  OppNTable[2*i  ]={IndexDn_4f[i],-IndexDn_4f[i],1}
end
OppN = NewOperator(NF, NB, OppNTable)
 
-- The Coulomb interaction
 
OppF0 =NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {1,0,0,0})
OppF2 =NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0,1,0,0})
OppF4 =NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0,0,1,0})
OppF6 =NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0,0,0,1})
 
-- The crystal-field operator
OppCF = NewOperator("CF", NF, IndexUp_4f, IndexDn_4f, Akm)
 
-- Given the operators and the parameters we create the Hamiltonian
 
H = NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {F0,F2,F4,F6})
  + zeta * Oppldots
  + OppCF
  + B[1] * (OppLx + 2*OppSx) + B[2] * (OppLy + 2*OppSy) + B[3] * (OppLz + 2*OppSz) 
H.Chop()
 
-- In the following we determine the eigenstates and set the chemical potential
-- to the center of the PES/IPES gap and set the ground-state energy to zero.
-- If the gap is negative you need to increase U.
 
StartRestrictions     = {NF, NB, {"11111111111111",ne,ne}}
StartRestrictionsmin  = {NF, NB, {"11111111111111",ne-1,ne-1}}
StartRestrictionsplus = {NF, NB, {"11111111111111",ne+1,ne+1}}
 
psiList = Eigensystem(H, StartRestrictions, NPsi, {{"Denseborder",3500}})
if Npsi==1 then psiList = {psiList} end
if ne ~= 0 then
  psimin = Eigensystem(H, StartRestrictionsmin, 1, {{"Denseborder",3500}})
end
if ne ~= 14 then
  psiplus = Eigensystem(H, StartRestrictionsplus, 1, {{"Denseborder",3500}})
end
 
if ne==0 then
  print("Shifts not implemented for ne=0")
elseif ne==14 then
  print("Shifts not implemented for ne=14")
else
  E0     = Chop(psiList[1] * H * psiList[1])
  E0min  = Chop(psimin     * H * psimin)
  E0plus = Chop(psiplus    * H * psiplus)
  De = (E0min - E0plus) / 2
  H = Chop(H + De * OppN - E0 - De * ne)
  E0     = Chop(psiList[1] * H * psiList[1])
  E0min  = Chop(psimin     * H * psimin)
  E0plus = Chop(psiplus    * H * psiplus)
  print("PES / IPES gap is ",E0min + E0plus - 2*E0)
  if E0min + E0plus - 2*E0 < 0 then
    print("With these parameters your gap is negative.")
    print("This is unphysical and an error in parameter space.")
    print("Increase U or choose a different starting valence.")
  end
end
 
-- print some expectation values
oppList={OppSsqr, OppLsqr, OppJsqr, OppSx, OppLx, OppJx, OppSy, OppLy, OppJy, OppSz, OppLz, OppJz, Oppldots, OppF2, OppF4, OppF6};
 
print("Initial state (before PES/IPES) state properties")
PrintExpectationValues(psiList, oppList, H)
if ne~=0 then
  print("Lowest electron removal state")
  PrintExpectationValues(psimin, oppList, H)
end
if ne~=14 then
  print("Lowest electron addition state")
  PrintExpectationValues(psiplus, oppList, H)
end
 
-- Calculate the electron removal Green's function
 
-- define the 14 operators removing one of the 4f electrons
TPES = {}
for i=0,13 do
  TPES[i+1] = NewOperator("An", NF, i)
end
 
-- next we generate the 14*14 tensor
-- <psi_0 | CR(tau) 1/(w-H+E0 +i 0^+) AN(tau') | psi_0>
-- The output S contains a spectrum object with 14*14 = 196 spectra.
--  https://www.quanty.org/documentation/language_reference/objects/spectra/start
-- The output G contains a response function object as a 14 by 14 matrix.
--  https://www.quanty.org/documentation/language_reference/objects/responsefunction/start
-- Both are stored as a table for the different initial states
-- You would need a Bolzman sum over the components of S or G
-- and then transform the Green's functions for the removal of an f-electron into the directional dependent intensity
-- for the f to d and or the f to g channels.
-- see also the example on the M45 PES in NiO
--  https://git.quanty.org/Quanty/tutorials/-/blob/main/08_Materials/NiO/01%20NiO%20Crystal%20field%20theory/18%20PES%20M45/18_PES_M45.Quanty?ref_type=heads
--  where we summed the transition operator to get the directional and polarisation dependent PES
--  instead of calcualing the spectra for one direction you can start from the Green's function and then calculate the directional dependent
--  PES from there
S={}
G={}
for i=1,NPsi do
  S[i], G[i] = CreateSpectra(H, TPES, psiList[i], {{"Emin",-2}, {"Emax",28}, {"NE",3000}, {"Gamma",0.1}, {"Tensor",true}})
end
 
-- In order to visualize the Green's function as a 14 by 14 matrix we can plot a matrix
-- Note you need a relative new version of Quanty to make the plots
-- You can download the develop source code here
--  https://www.quanty.org/download/git/develop
 
pl={}
Ymin = -6
Ymax =  6
dYtick = 2
Emin   = -2
Emax   =  10
dETick = 2
 
for i=1,NPsi do
  pl[i] = Graphics.Plot(S[i],{{"Nrow",14},{"Ncolumn",14},
                              {"Frame",{{"Ymin",Ymin},{"Ymax",Ymax},{"dYTick",dYtick},{"Xmin",Emin},{"Xmax",Emax},{"dXTick",dETick},{"XLabel","Energy"},{"YLabel","Intensity"},{"FontSize",0.01}}}})
  SVGPlG = Graphics.ToSVG(pl[i],{{"RelativeSize",true}})
  file = io.open(directory.."GPES_"..i..".svg",'w')
  file:write(SVGPlG)
  file:close()
end

Note that we store spectra as spectra objects in S and as response function objects in G. The response functions can be evaluated at a particular energy, or seen as a sum over poles. The poles correspond to eigenstates of a Krylov basis. The Krylov basis only corresponds to the full basis if you make the Krylov basis large enough. That is generally very expensive to calculate.

You can find an example how to get the polarization and directional intensity in this example for NiO. Note that we there generated a polarisation and directional dependent operator. As you want all directions it would be much faster to calculate the one particle electron removal Green's function as a matrix and do the polarisation and directional dependent calculation from this. (Note you can sum and rotate response function objects so this is relatively cheap and straight forward to do in Quanty).

Hope this helps, Best wishes, Maurits