NRG for Spectral Functions
Please download the dependencies for the following example here.
%% Dynamical impurity spin susceptibility
% Perform the iterative diagonalization calculation, with the same parameters
% chosen in the demonstration of the previous examples.
clear
% Hamiltonian parameters
U = 4e-3; % Coulomb interaction at the impurity
epsd = -U/2; % impurity on-site energy
Gamma = 8e-5*pi; % hybridization strength
% NRG parameters
Lambda = 2.5; % discretization parameter
N = 55; % length of the Wilson chain
Nkeep = 300;
% Wilson chain
[ff,gg] = doCLD([-1 1],[1 1]*Gamma/pi,Lambda,N);
% symmetries
symstr = 'Acharge,SU2spin'; % U(1) charge and SU(2) spin
% symstr = 'Acharge,Aspin'; % U(1) charge and U(1) spin
% Construct local operators
[F,Z,S,I] = getLocalSpace('FermionS',symstr,'NC',1);
[F,Z,S,EF] = setItag('s00','op',F,Z,S,I.E);
% particle number operator
NF = QSpace;
for itF = 1:numel(F)
NF(itF) = contract(F(itF),'!2*',F(itF));
end
% Impurity Hamiltonian
H0 = U/2*sum(NF)*(sum(NF)-1) + epsd*sum(NF) + 1e-33*EF;
% ket tensor for the impurity
A0 = getIdentity(setItag('L00',getvac(EF,1)),1,EF,1,'K00',[1,3,2]);
H0 = contract(A0,'!2*',{A0,H0});
% same hopping amplitude and on-site energies for all flavors
ff = repmat(ff,[1,numel(F)]);
gg = repmat(gg,[1,numel(F)]);
% iterative diagonalization
Inrg = NRG_IterDiagQS(H0,A0,Lambda,ff,F,gg,NF,Z,Nkeep);
T = 1e-8;
Inrg = getRhoFDMQS(Inrg,T); % construct full density matrix
%%
% Compute the discrete spectral data of the dynamical impurity spin susceptibility
% $\chi_s (\omega)$. Here, we use the spin-z operator |S(3)| in case of U(1) spin
% symmetry or |S| in case of SU(2) spin symmetry. Note that in case of SU(2) spin
% symmetry, we compute the sum of spin correlations in all three spatial directions,
% which are all identical. To obtain only the z-component, we therefore have to
% divide by 3. By setting the third input as empty (|[]|), the |getAdiscQS| routine
% understands that the input operator is *bosonic*. For fermionic operators, the
% so-called Z string should be considered, and the sign factor between two terms
% in the anti-commutator is $+$. On the other hand, for bosonic operators, Z operators
% do not involve, and the sign factor is $-$ due to commutator.
% dynamical impurity spin susceptibility
if numel(S)>1 % U(1) spin
[odisc,Adisc] = getAdiscQS(Inrg,S(3),[]);
else % SU(2) spin
[odisc,Adisc] = getAdiscQS(Inrg,S,[]);
Adisc = Adisc/3; % divide by 3 to consider only z-direction
end
%%
% Broaden the discrete data.
% broaden the discrete data to have a continuous curve
[ocont,Acont] = getAcont(odisc,Adisc,log(Lambda),T/5);
%%
% Plot the result.
figure;
hold on;
% positive frequency
plot(ocont(ocont>0),Acont(ocont>0), ...
'LineWidth',1,'LineStyle','-');
% negative frequency
plot(-ocont(ocont<0),-Acont(ocont<0), ...
'LineWidth',1,'LineStyle','--');
hold off;
set(gca,'FontSize',13,'LineWidth',1,'XScale','log','YScale','log');
legend({'$\chi''''_s(\omega > 0)$','$-\chi''''_s(\omega < 0)$'}, ...
'Interpreter','latex');
xlabel('$| \omega |$','Interpreter','latex');
ylabel('$\chi''''_s$','Interpreter','latex');
xlim([1e-11 1]);
grid on;
%%
% Note that the minus sign prefactor in $- \chi''_s (\omega < 0)$ to visualize
% the data for negative frequency. We see that the curve is an odd function, i.e.,
% $\chi''_s (\omega) = - \chi''_s (-\omega)$. It is the generic property of the
% imaginary part of the correlation functions of bosonic operators.
%%
% Identify the peak position of the curve.
[~,maxid] = max(Acont);
disp(ocont(maxid));
%%
% This value is similar to the Kondo temperature $T_\mathrm{K}$ that we obtained
% from the Bethe ansatz solution. (Note that there are various ways of computing
% the Kondo temperature, and these ways give similar but different values.)
TK = sqrt(U*Gamma/2)*exp(-pi*U/8/Gamma + pi*Gamma/2/U);
disp(TK);
disp(ocont(maxid)/TK); % ratio
%%
% Indeed, the peak position of $\chi''_s (\omega)$ is the method of choice to
% determine the Kondo temperature for general quantum impurity systems! This method
% is especially useful when analytical approaches (such as poor man's scaling
% and the Bethe ansatz) are not applicable due to the complexity of the system.