-- Eight fermionic modes with four electrons. H0 is a hopping over all eight modes and
-- H1 = U n_0 n_1 is a Hubbard interaction between mode 0 and mode 1, so the self-energy
-- lives on the two modes 0 and 1 and the other six are the bath.
NF = 8
h = Matrix.New({
 { 0.35,-0.62, 0.18, 0.41,-0.27, 0.09, 0.33,-0.15},
 {-0.62,-0.21, 0.54,-0.13, 0.46, 0.28,-0.37, 0.22},
 { 0.18, 0.54, 0.77, 0.31,-0.19,-0.44, 0.12, 0.38},
 { 0.41,-0.13, 0.31,-0.58, 0.24, 0.16, 0.49,-0.26},
 {-0.27, 0.46,-0.19, 0.24, 1.12,-0.35, 0.21, 0.14},
 { 0.09, 0.28,-0.44, 0.16,-0.35,-0.93, 0.29, 0.43},
 { 0.33,-0.37, 0.12, 0.49, 0.21, 0.29, 0.66,-0.18},
 {-0.15, 0.22, 0.38,-0.26, 0.14, 0.43,-0.18, 0.05}})
U  = 2.5
H0 = Matrix.ToOperator(h)
H1 = U * NewOperator(NF,0,{{0,-0,1}}) * NewOperator(NF,0,{{1,-1,1}})
psi = Eigensystem(H0 + H1, {NF,0,{"11111111",4,4}}, 1)

-- The augmented propagator of the operators a_m and q_m = [a_m,H1], m = 0,1. Its first
-- block is the Green's function itself and its second block is what turns it into the
-- self-energy.
Cr = {NewOperator(NF,0,{{ 0,1}}), NewOperator(NF,0,{{ 1,1}})}
An = {NewOperator(NF,0,{{-0,1}}), NewOperator(NF,0,{{-1,1}})}
Cr[3] = H1*Cr[1] - Cr[1]*H1  -- q^dag = [H1,a^dag]
Cr[4] = H1*Cr[2] - Cr[2]*H1
An[3] = An[1]*H1 - H1*An[1]  -- q     = [a,H1]
An[4] = An[2]*H1 - H1*An[2]
opt = {{"Tensor",true},{"NTri",24},{"DenseBorder",0}}
Sp, Gp = CreateSpectra(H0+H1, Cr, psi, opt)
Sm, Gm = CreateSpectra(H0+H1, An, psi, opt)
Gtilde = Gp + ResponseFunction.InvertEnergy(Gm)   -- G(w) = Gp(w) - Gm(-w)^T

-- The Hartree-Fock part is not part of what Downfold returns, it is an expectation value
SigmaHF = Matrix.Zero(2)
for i = 1, 2 do
  for j = 1, 2 do
    SigmaHF[i][j] = psi * (An[i+2]*Cr[j] + Cr[j]*An[i+2]) * psi
  end
end

-- and the rest of the self-energy is the lower right block of the inverse of the augmented
-- propagator, inverted again
R = Matrix.New({{0,0,1,0},{0,0,0,1}})
Sigma, S = ResponseFunction.Downfold(Gtilde, R)

print("Sigma^HF = <{q,a^dag}>")
print(SigmaHF)
print("S, the spectral weight of Sigma - Sigma^HF")
print(S)
print("")
print("  omega      Sigma_00 from Downfold        from G0^-1 - G^-1")
for w = -4, 4, 2 do
  local G0    = Matrix.Sub(Matrix.Inverse(Complex.New(w,0.05)*Matrix.Identity(NF) - h), 2)
  local G     = Matrix.Sub(Gtilde(w,0.1), 2)
  local DSR   = SigmaHF + Sigma(w,0.1)
  local Dyson = Matrix.Inverse(G0) - Matrix.Inverse(G)
  print(string.format("  %5.1f   %9.5f %+9.5f I   %9.5f %+9.5f I", w,
        Complex.Re(DSR[1][1]),   Complex.Im(DSR[1][1]),
        Complex.Re(Dyson[1][1]), Complex.Im(Dyson[1][1])))
end

-- The same self-energy as poles on the grid -4, -3.98, ... 4, taken from Gtilde off the
-- real axis. Every pole of the exact Sigma is split over its two neighbouring grid points.
SigmaM, SM = ResponseFunction.Downfold(Gtilde, R, {{"Method","Mesh"},{"Emin",-4},{"Emax",4},{"NE",400}})
print("")
print("  omega      Sigma_00 - Sigma^HF, default     Method \"Mesh\"")
for w = -4, 4, 2 do
  local a, b = Sigma(w,0.1)[1][1], SigmaM(w,0.1)[1][1]
  print(string.format("  %5.1f   %9.5f %+9.5f I   %9.5f %+9.5f I", w,
        Complex.Re(a), Complex.Im(a), Complex.Re(b), Complex.Im(b)))
end