%Fig_6_7.m % %Drude dielectric % clear all; clf; t=cputime; eye=complex(0.,1.); %square root of -1 %plotting parameters + fontsizes FS = 18; %label fontsize 18 FSN = 16; %number fontsize 16 LW = 2; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman') set(0,'DefaultAxesFontSize', FSN) % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman') set(0,'DefaultTextFontSize', FSN) hbar=['\fontname{MT Extra}h\fontname{Times New Roman}']; clight =2.99792458e8; %speed of light in vacuum [m s^-1] echarge =1.60219e-19; %electron charge [C] epsilon0 =8.8541878e-12;%permittivity of free-space [F m^-1] mu0 =4*pi*1e-7; %periability of free space [H m^-1] m0 =9.10956e-31; %bare electron mass [kg] meff =1.01*m0; %effective electron mass in Cu [kg] nelectrons =8.46e28; %electron density in Cu [m^-3] sigma0 =1.0*5.9e7; %DC conductity of Cu at room temperature [S m^-1] tau =sigma0*meff/(nelectrons*(echarge^2));%scattering time [s] delta =eye*1.0/tau; %broadening omegap2 =nelectrons*(echarge^2)/(meff*epsilon0);%plasma frequency squared omegap =sqrt(omegap2)%plasma frequency [rad s^-1] omegapeV =sqrt(omegap2)*6.58211889e-16%plasma energy [eV] omegap2tau =omegap2*tau; omegap2tau2 =omegap2*tau^2; omegamin=2*pi*1000e6; %minimum frequency in plot 2*pi*1000e6 2*pi*1000e12 [rad s^-1] omegamax=2*pi*10000e12; %maximum frequency in plot [rad s^-1] npoints=100000; omega =linspace(omegamin,omegamax,npoints); %frequency [rad s^-1] omega2 =omega.^2; omega2tau2 =omega2*tau^2; klight =omega/clight; epsilon = 1-(sigma0*tau./(epsilon0*(1+omega2tau2))); epsilon=epsilon+(eye*sigma0./(epsilon0*omega.*(1+omega2tau2))); kx =sqrt(omega2.*epsilon)/clight; kkx =sqrt(omega2-omegap2)/clight; nr =clight*real(kx)./(omega); aclight =omega./klight; vgroup=zeros(npoints,1); for i=1:1:npoints-1 vgroup(i) =(omega(i+1)-omega(i))/(real(kx(i+1))-real(kx(i))); end vgroup(npoints)=vgroup(npoints-1); ttl=['\sigma_{0} = ',... num2str(sigma0,'%3.2e'),' S m^{-1}, \it{n}\rm_0 = ',... num2str(nelectrons/1e28,'%3.2f'),'\times10^{28} m^{-3}, \tau = ',... num2str(tau*1e15,'%2.0f'),' fs, \it{E}\rm_p = ',... num2str(omegapeV,'%2.1f'),' eV']; figure(1); %plot imag k plot(log10(imag(kx)),log10(omega./(2e9*pi)),'r','LineWidth',LW); hold on %plot real k plot(log10(real(kx)),log10(omega./(2e9*pi)),'b','LineWidth',LW); %plot k in vacuum plot(log10(klight),log10(omega./(2e9*pi)),':r','LineWidth',LW); %grid on ylabel('Frequency, log_{10}(\it{f}\rm) (GHz)'); xlabel('log_{10}(\it{k}\rm) (m^{-1})'); axis([0, 10,log10(omegamin/(2e9*pi)),log10(omegamax/(2e9*pi))]) title(ttl); hold off