CDW谱函数数值计算

代码

main

dk = 0.01;
%kx = 0:dk:2.0943952;
%kx = 0:dk:4.18879028;
kx = 0;
eps = -1:0.005:2;
%ky = -3.62759876;
%ky = 0;
ky = 0:dk:3.62759876;
delta = 0.03;
gammar = 0.02;
% E0 = 0;
[E1,E2,Ek,Ekq]= get_E1_and_E2(kx,ky,delta);
A=[];
for e=-1:0.005:2
    A0 = spectral_function(e,E1,E2,Ekq,delta,gammar);
    A = [A;A0];
end
[X,Y] = meshgrid(ky,eps);
mesh(X,Y,A)

对角化计算Ek

function E = get_energy(kx,ky)
n = length(ky);
n
E = zeros(1,n);
for i=1:n
    H=Hamiltonian_TaS2_k(kx,ky(i));
    [~,Ek]=eig(H);
    E(1,i) = Ek(1,1);
end

计算对应的Ek，Ekq，E1，E2

function [E1,E2,Ek,Ekq]= get_E1_and_E2(kx,ky,delta)
%Ek = -cos(kx);
Ek = get_energy(kx,ky);
%Ekq = -cos(kx+pi);
Ekq = get_energy(kx,ky-3.62759876*ones(1,length(ky)));
E1 = (Ek+Ekq)/2 + sqrt((Ek-Ekq).^2/4 + delta^2);
E2 = (Ek+Ekq)/2 - sqrt((Ek-Ekq).^2/4 + delta^2);

计算谱函数，delta函数用归一化的洛伦兹函数近似。

function A=spectral_function(E0,E1,E2,Ekq,delta,gammar)
C1 = 1+delta^2./(E1-Ekq).^2;
C2 = 1+delta^2./(E2-Ekq).^2;
n = length(E1);
E = E0*ones(1,n);
A1 = (gammar/pi)*((E-E1).^2+gammar^2*ones(1,n)).^(-1);
A2 = (gammar/pi)*((E-E2).^2+gammar^2*ones(1,n)).^(-1);
A = A1./C1 + A2./C2;