RESPES in quanty

I would like to calculate 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 Er: \[4d^{10}4f^{11} \xrightarrow{h\nu} 4d^{10}4f^{10}\] \[4d^{10}4f^{11} \xrightarrow{h\nu} 4d^{9}4f^{12}\xrightarrow{e^2/r} 4d^{10}4f^{10}\]

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:

AUGERSpectra, AUGERResponse = 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},
    })
AUGERSpectra.Print({{"file", "Sm4d4f_AUGER.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 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
    -- 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"}})


AUGERSpectra, AUGERResponse = 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},
    })
AUGERSpectra.Print({{"file", "Er4d4f_AUGER.dat"}})

-- combine the two spectra i.e. |AUGERSpectra+XPSSpectra|^2 (I've been doing that in python, not the problem here)