RESPES in quanty

I would like to calcuate resonant X-ray photoemission spectra in Quanty. The resonance in question is between a normal XPS process and and auger-like process, e.g. for the d-f resonance of Yb: \[4d^{10}4f^13 \xrightarrow{h\nu} 4d^{10}4f^12\] \[4d^{10}4f^13 \xrightarrow{h\nu} 4d^{9}4f^14\xrightarrow{e^2/r} 4d^{10}4f^12\]

The first process is trivial in quanty: XPSSpectra , XPSResponse = CreateSpectra(finalstate.Hamiltonian, T_4f, groundstate.psiList[1], emin_-2_emax_10_ne_1200_gamma_0.1) XPSSpectra.Print(file_sm4d4f_xps.dat) Where T_4f is a list of annihilation operators on the 4f shell. For the second process I have:

RIXSSpectra, RIXSResponse = CreateResonantSpectra(

  augerstate.Hamiltonian, finalstate.Hamiltonian, T_4d4f_xyz[1], T_4f, groundstate.psiList[1], {
      -- {"restrictions1", augerstate.StartRestrictions}, {"restrictions2", finalstate.StartRestrictions},
      -- {"E0_1", augerstate.energy-groundstate.energy}, {"E0_2", finalstate.energy-groundstate.energy},
      {"Emin1",-15}, {"Emax1",35}, {"NE1",120}, {"Gamma1",0.1},
      {"Emin2",-2}, {"Emax2",10}, {"NE2",1200}, {"Gamma2",0.1},
  })

RIXSSpectra.Print(file_sm4d4f_rixs.dat)

However this is problematic as the photoexciation operator produces a totally imaginary spectrum in the first axis, and so there can be no Fano resonance (which gives RESPES it's name) when the two spectra are combined.

Are these the wrong operators, or am I making a more fundamental mistake? Many thanks, Ted Trewick

I have pasted a MWE of the calculation below:

NF = 10 + 14 NB = 0

IndexDn_4d = {0,2,4,6,8} IndexUp_4d = {1,3,5,7,9} IndexDn_4f = {10,12,14,16,18,20,22} IndexUp_4f = {11,13,15,17,19,21,23}

-- L.S (spin orbit coupling) Oppldots_4f = NewOperator(“ldots”, NF, IndexUp_4f, IndexDn_4f) Oppldots_4d = NewOperator(“ldots”, NF, IndexUp_4d, IndexDn_4d)

function make_CF_AkmOh(A)

  return {
  {0, 0,                      A[0]}, 
  {4,-4,       (sqrt(5/14))*(A[4])}, 
  {4, 0,                      A[4]}, 
  {4, 4,       (sqrt(5/14))*(A[4])}, 
  {6,-4, (-1)*((sqrt(7/2))*(A[6]))}, 
  {6, 0,                      A[6]}, 
  {6, 4, (-1)*((sqrt(7/2))*(A[6]))}}

end

groundstate = {

  name ='Er 4f11',
  energy = -355606.3853,
  Nfelectrons = 11,
  Ndelectrons = 10,
  Npsi = 364,
  F4f4f_F0F2F4F6 = {0, 16.0911*0.8, 10.0941*0.8, 7.2615*0.8},
  zeta_4f = 0.3022,
  A_Oh = {[0]=0, [4]=0.3, [6]=0.02},

} finalstate = {

  name = 'Er 4f10',
  energy = -355564.1659,
  Nfelectrons = 12,
  Ndelectrons = 10,
  Npsi = 1001,
  F4f4f_F0F2F4F6 = {0, 16.9639*0.8, 10.6870*0.8, 7.7009*0.8},
  zeta_4f = 0.3210,
  A_Oh = {[0]=0, [4]=0.3, [6]=0.02},

}

augerstate = {

  name = 'Er 4d9,4f12',
  energy = -355431.8348,
  Nfelectrons = 12,
  Ndelectrons = 9,
  Npsi = 91,
  F4f4f_F0F2F4F6 = {0, 16.2229*0.8, 10.1820*0.8, 10.1820*0.8},
  zeta_4f = 0.3065,
  F4d4f_F0F2F4 = {0, 18.0665*0.8, 11.5472*0.8},
  G4d4f_G1G3G5 = {21.2249*0.7, 13.3654*0.7, 9.4608*0.7},
  zeta_4d = 3.0449,
  A_Oh = {[0]=0, [4]=0.3, [6]=0.02},

}

for _, Ln in ipairs({groundstate, finalstate}) do

  Ln.OppCFOh   = NewOperator("CF", NF, IndexUp_4f, IndexDn_4f, make_CF_AkmOh(Ln.A_Oh))
  Ln.OppJin    = NewOperator("U",  NF, IndexUp_4f, IndexDn_4f, Ln.F4f4f_F0F2F4F6)
  Ln.Hamiltonian = Ln.OppJin + Ln.zeta_4f*Oppldots_4f + Ln.OppCFOh

  Ln.StartRestrictions = {NF, NB,
      {"1111111111"  .. "00000000000000", Ln.Ndelectrons, Ln.Ndelectrons},  -- 4d
      {"0000000000"  .. "11111111111111", Ln.Nfelectrons, Ln.Nfelectrons},  -- 4f
  }
  Ln.psiList = Eigensystem(Ln.Hamiltonian, Ln.StartRestrictions, Ln.Npsi)

end

do

  local Ln = augerstate

  Ln.OppCFOh = NewOperator("CF", NF, IndexUp_4f, IndexDn_4f, make_CF_AkmOh(Ln.A_Oh))
  Ln.OppJin  = NewOperator("U",  NF, IndexUp_4f, IndexDn_4f, Ln.F4f4f_F0F2F4F6)
  Ln.OppUpd  = NewOperator("U", NF, IndexUp_4d, IndexDn_4d, IndexUp_4f, IndexDn_4f, Ln.F4d4f_F0F2F4, Ln.G4d4f_G1G3G5)

  Ln.Hamiltonian = Ln.OppJin + Ln.zeta_4f * Oppldots_4f + Ln.zeta_4d * Oppldots_4d + Ln.OppCFOh + Ln.OppUpd
  Ln.Hamiltonian = Chop(Ln.Hamiltonian)
  Ln.StartRestrictions = {NF, NB,
      {"1111111111"  .. "00000000000000", Ln.Ndelectrons, Ln.Ndelectrons},
      {"0000000000"  .. "11111111111111", Ln.Nfelectrons, Ln.Nfelectrons},
  }
  -- Ln.psiList = Eigensystem(Ln.Hamiltonian, Ln.StartRestrictions, Ln.Npsi)

end

T_4d4f_xyz = {} T_4d4f_average = 0 T_4f4d_xyz = {} T_4f4d_average = 0 for i, e in ipairs(1.0_0.0_0.0_0.0_1.0_0.0_0.0_0.0_1.0) do

1. - https://www.quanty.org/documentation/standard_operators/transition_operators_for_photon_spectroscopy polarisation_Akm = { {1, -1,( e[1] + I * e[2])/math.sqrt(2)}, {1, 0, e[3]}, {1, 1,(-e[1] + I * e[2])/math.sqrt(2)}} T_4f4d_xyz[i] = NewOperator(“CF”, NF, IndexUp_4d, IndexDn_4d, IndexUp_4f, IndexDn_4f, polarisation_Akm) T_4f4d_average = T_4f4d_average + T_4f4d_xyz[i] T_4d4f_xyz[i] = NewOperator(“CF”, NF, IndexUp_4f, IndexDn_4f, IndexUp_4d, IndexDn_4d, polarisation_Akm) T_4d4f_average = T_4d4f_average + T_4d4f_xyz[i]

end

T_4f = {} for i = 1, 14 / 2 do

  T_4f[2*i - 1] = NewOperator('An', NF, IndexDn_4f[i])
  T_4f[2*i]     = NewOperator('An', NF, IndexUp_4f[i])

end

XPSSpectra , XPSResponse = CreateSpectra(finalstate.Hamiltonian, T_4f, groundstate.psiList[1], emin_-2_emax_10_ne_1200_gamma_0.1) XPSSpectra.Print(file_er4d4f_xps.dat)

RIXSSpectra, RIXSResponse = CreateResonantSpectra(

  augerstate.Hamiltonian, finalstate.Hamiltonian, T_4d4f_xyz[1], T_4f, groundstate.psiList[1], {
      -- {"restrictions1", augerstate.StartRestrictions}, {"restrictions2", finalstate.StartRestrictions},
      -- {"E0_1", augerstate.energy-groundstate.energy}, {"E0_2", finalstate.energy-groundstate.energy},
      {"Emin1",-15}, {"Emax1",35}, {"NE1",120}, {"Gamma1",0.1},
      {"Emin2",-2}, {"Emax2",10}, {"NE2",1200}, {"Gamma2",0.1},
  })

RIXSSpectra.Print(file_er4d4f_rixs.dat)