ResponseFunction.CalculateSelfEnergy($G_0$,$G$) calculates \begin{equation} \Sigma = G_0^{-1} - G^{-1} \end{equation}
As an example we take a one dimensional tight binding Hamiltonian of a 4 site ring for $H_0$ and $G_0$. To this we add to each site a one dimensional chain of length 2. In most real cases $H_1$ would be a two particle interaction. This example allows one to see the effect of the different algorithms clearly.
-- We define H0 as a 4 site tight binding model with periodic boundary conditions -- -- ... (0) --- (3) --- (6) --- (9) --- (0) ... -- H0 = NewOperator(12,0,{{ 0, -3,1},{ 3, -0,1}, { 3, -6,1},{ 6, -3,1}, { 6, -9,1},{ 9, -6,1}, { 9, -0,1},{ 0, -9,1}}) print("The tight binding Hamiltonian for a 4 cite chain is given as") print(H0) print("In matrix form") print(OperatorToMatrix(H0)) -- We define H1 as coupling to the local sites -- -- (0) (3) (6) (9) -- | | | | -- (1) (4) (7) (10) -- | | | | -- (2) (5) (8) (11) -- H1 = NewOperator(12,0,{{ 0, -1,1},{ 1, -0,1}, { 1, -2,1},{ 2, -1,1}, { 3, -4,1},{ 4, -3,1}, { 4, -5,1},{ 5, -4,1}, { 6, -7,1},{ 7, -6,1}, { 7, -8,1},{ 8, -7,1}, { 9,-10,1},{10, -9,1}, { 10,-11,1},{11,-10,1}}) print("The full Hamiltonian we take as the one dimensional chain with side chains") print("The side chains are given by the Hamiltonian") print(H1) print("In matrix form") print(OperatorToMatrix(H1)) -- We now can define the full Hamiltonian H = H0 + H1 print("The full Hamiltonian is H0 + H1") print(H) print("In matrix form") print(OperatorToMatrix(H)) print("We now calcualte the bare (G0) and full (G) vacuum Green's functions") -- In order to define the Green's function G0 and G we need the vacuum state psi0 = NewWavefunction(12,0,{{"000 000 000 000",1}}) -- And the transition operators -- We want a^{dag}_k = sum_{x=1}^N \sqrt{1/N} e^{i k x} a^{dag}_x adag={} x = {0,1,2,3} k = {0/4,1/4,2/4,3/4} for ik=1,4 do adag[ik] = NewOperator(12,0,{{ 0, sqrt(1/4)*exp(I*x[1]*k[ik]*2*Pi)}, { 3, sqrt(1/4)*exp(I*x[2]*k[ik]*2*Pi)}, { 6, sqrt(1/4)*exp(I*x[3]*k[ik]*2*Pi)}, { 9, sqrt(1/4)*exp(I*x[4]*k[ik]*2*Pi)}}) end -- Create Green's function S0, G0 = CreateSpectra(H0, adag, psi0,{{"Tensor",true},{"DenseBorder",0}}) G0.Chop() S, G = CreateSpectra(H , adag, psi0,{{"Tensor",true},{"DenseBorder",0}}) G.Chop() print("G0 is") print(G0) print("in matrix form that becomes") print(ResponseFunction.ToMatrix(G0)) print("G is") print(G) print("in matrix form that becomes") print(ResponseFunction.ToMatrix(G)) print("The self energy can be calcualted as Sigma = G0^{-1} - G^{-1}") Sigma = ResponseFunction.CalculateSelfEnergy(G0,G) Sigma = ResponseFunction.ChangeType(Sigma,"Tri") Sigma.Chop() print(Sigma) print("We now can compare G0(w), G(w) and G0(w-Sigma(w))") omega = 0.1 Gamma = 0.3 print("G(omega)") print(G(omega,Gamma/2)) print("G(omega-Sigma(omega))") print(Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) )) print("G(omega)-Sigma(omega)") print(G(omega,Gamma/2) - Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) )) print("The previous should have been zero but computer math is different. Have a look at") print("0.1+0.2-0.3") print(0.1+0.2-0.3) print("If we \"Chop\" the previous result it becomes zero") print(Chop( G(omega,Gamma/2) - Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) )) )
The tight binding Hamiltonian for a 4 cite chain is given as
Operator: Operator
QComplex = 0 (Real==0 or Complex==1 or Mixed==2)
MaxLength = 2 (largest number of product of lader operators)
NFermionic modes = 12 (Number of fermionic modes (site, spin, orbital, ...) in the one particle basis)
NBosonic modes = 0 (Number of bosonic modes (phonon modes, ...) in the one particle basis)
Operator of Length 2
QComplex = 0 (Real==0 or Complex==1)
N = 8 (number of operators of length 2)
C 0 A 3 | 1.00000000000000E+00
C 3 A 0 | 1.00000000000000E+00
C 3 A 6 | 1.00000000000000E+00
C 6 A 3 | 1.00000000000000E+00
C 6 A 9 | 1.00000000000000E+00
C 9 A 6 | 1.00000000000000E+00
C 9 A 0 | 1.00000000000000E+00
C 0 A 9 | 1.00000000000000E+00
In matrix form
{ { 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } }
The full Hamiltonian we take as the one dimensional chain with side chains
The side chains are given by the Hamiltonian
Operator: Operator
QComplex = 0 (Real==0 or Complex==1 or Mixed==2)
MaxLength = 2 (largest number of product of lader operators)
NFermionic modes = 12 (Number of fermionic modes (site, spin, orbital, ...) in the one particle basis)
NBosonic modes = 0 (Number of bosonic modes (phonon modes, ...) in the one particle basis)
Operator of Length 2
QComplex = 0 (Real==0 or Complex==1)
N = 16 (number of operators of length 2)
C 0 A 1 | 1.00000000000000E+00
C 1 A 0 | 1.00000000000000E+00
C 1 A 2 | 1.00000000000000E+00
C 2 A 1 | 1.00000000000000E+00
C 3 A 4 | 1.00000000000000E+00
C 4 A 3 | 1.00000000000000E+00
C 4 A 5 | 1.00000000000000E+00
C 5 A 4 | 1.00000000000000E+00
C 6 A 7 | 1.00000000000000E+00
C 7 A 6 | 1.00000000000000E+00
C 7 A 8 | 1.00000000000000E+00
C 8 A 7 | 1.00000000000000E+00
C 9 A 10 | 1.00000000000000E+00
C 10 A 9 | 1.00000000000000E+00
C 10 A 11 | 1.00000000000000E+00
C 11 A 10 | 1.00000000000000E+00
In matrix form
{ { 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 } }
The full Hamiltonian is H0 + H1
Operator: Operator
QComplex = 0 (Real==0 or Complex==1 or Mixed==2)
MaxLength = 2 (largest number of product of lader operators)
NFermionic modes = 12 (Number of fermionic modes (site, spin, orbital, ...) in the one particle basis)
NBosonic modes = 0 (Number of bosonic modes (phonon modes, ...) in the one particle basis)
Operator of Length 2
QComplex = 0 (Real==0 or Complex==1)
N = 24 (number of operators of length 2)
C 0 A 3 | 1.00000000000000E+00
C 3 A 0 | 1.00000000000000E+00
C 3 A 6 | 1.00000000000000E+00
C 6 A 3 | 1.00000000000000E+00
C 6 A 9 | 1.00000000000000E+00
C 9 A 6 | 1.00000000000000E+00
C 9 A 0 | 1.00000000000000E+00
C 0 A 9 | 1.00000000000000E+00
C 0 A 1 | 1.00000000000000E+00
C 1 A 0 | 1.00000000000000E+00
C 1 A 2 | 1.00000000000000E+00
C 2 A 1 | 1.00000000000000E+00
C 3 A 4 | 1.00000000000000E+00
C 4 A 3 | 1.00000000000000E+00
C 4 A 5 | 1.00000000000000E+00
C 5 A 4 | 1.00000000000000E+00
C 6 A 7 | 1.00000000000000E+00
C 7 A 6 | 1.00000000000000E+00
C 7 A 8 | 1.00000000000000E+00
C 8 A 7 | 1.00000000000000E+00
C 9 A 10 | 1.00000000000000E+00
C 10 A 9 | 1.00000000000000E+00
C 10 A 11 | 1.00000000000000E+00
C 11 A 10 | 1.00000000000000E+00
In matrix form
{ { 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 } ,
{ 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 } }
We now calcualte the bare (G0) and full (G) vacuum Green's functions
Start of LanczosBlockTriDiagonalize
Start of LanczosBlockTriDiagonalizeRC
Reduce dimension from 4 to 0 in block 1
Start of LanczosBlockTriDiagonalize
Start of LanczosBlockTriDiagonalizeRC
Reduce dimension from 4 to 0 in block 3
G0 is
{ { { { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 2 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , -2 , 0 } ,
{ 0 , 0 , 0 , 0 } } } ,
{ { { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } } ,
BlockSize = { 4 , 4 } ,
type = Tri ,
mu = 0 ,
name = Block Tridiagonal Matrix }
in matrix form that becomes
{ { 2 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , -2 , 0 } ,
{ 0 , 0 , 0 , 0 } }
G is
{ { { { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 2 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , -2 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } } ,
{ { { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } ,
{ { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } ,
{ { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } } ,
BlockSize = { 4 , 4 , 4 , 4 } ,
type = Tri ,
mu = 0 ,
name = Block Tridiagonal Matrix }
in matrix form that becomes
{ { 2 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , -2 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 } ,
{ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 } ,
{ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 } }
The self energy can be calcualted as Sigma = G0^{-1} - G^{-1}
Reduce dimension from 4 to 0 in block 2
{ { { { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } ,
{ { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } } } ,
{ { { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } ,
{ { 1 , 0 , 0 , 0 } ,
{ 0 , 1 , 0 , 0 } ,
{ 0 , 0 , 1 , 0 } ,
{ 0 , 0 , 0 , 1 } } } ,
BlockSize = { 4 , 4 , 4 } ,
type = Tri ,
mu = 0 ,
name = Matrix }
We now can compare G0(w), G(w) and G0(w-Sigma(w))
G(omega)
{ { (-0.55141413226786 - 0.046491809053112 I) , 0 , 0 , 0 } ,
{ 0 , (3.1750878423276 - 2.418978356337 I) , 0 , 0 } ,
{ 0 , 0 , (0.4525371980109 - 0.031240475239916 I) , 0 } ,
{ 0 , 0 , 0 , (3.1750878423276 - 2.418978356337 I) } }
G(omega-Sigma(omega))
{ { (-0.55141413226786 - 0.046491809053112 I) , 0 , 0 , 0 } ,
{ 0 , (3.1750878423276 - 2.418978356337 I) , 0 , 0 } ,
{ 0 , 0 , (0.4525371980109 - 0.031240475239916 I) , 0 } ,
{ -0 , -0 , -0 , (3.1750878423276 - 2.418978356337 I) } }
G(omega)-Sigma(omega)
{ { (-1.1102230246252e-16 - 8.3266726846887e-17 I) , 0 , 0 , 0 } ,
{ 0 , (-8.8817841970013e-16 - 4.4408920985006e-16 I) , 0 , 0 } ,
{ 0 , 0 , (-5.5511151231258e-17 - 4.5102810375397e-17 I) , 0 } ,
{ 0 , 0 , 0 , 4.4408920985006e-16 } }
The previous should have been zero but computer math is different. Have a look at
0.1+0.2-0.3
5.5511151231258e-17
If we "Chop" the previous result it becomes zero
{ { 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } ,
{ 0 , 0 , 0 , 0 } }