Supplementary Material

This is the supplementary material for Surface and Interfacial Dynamics of Polymeric Bilayer Films. It contains matlab codes and script for numerical and approximated analytical calculations of the over-damped capillary wave relaxation time constants for a supported PS film via bilayer model described in the paper. fzerosscpt.m and sscptker.m are function codes for numerical calcuations, tau_fast.m and tau_slow.m for analytical calcuations.

Contents

Bilayer Film Parameters

Susceptibility Map

Over-damped Relaxation Time Constants

Mode amplitude

Plot Relaxation Time Constants

Bilayer Film Parameters

For bilayer system, define layer thicknesses, viscosities, surface/interfacial tensions and mass densities. For the PS system considered in the paper, interfacial tension is zero. Mass densities do not affect capillary wave relaxation time constants. Thus to prevent numerical erros, two arbitrary values are used in the calculation.

d = 80; % total film thickness (nm)
h1 = 70; % bottom layer thickness (nm)
eta1 = 4.4e5/1e9; % bottom layer viscosity (kg/s/nm)
eta2 = 4.4e4/1e9; % top layer viscosity (kg/s/nm)
r1 = 0; % interface tension (kg/s^2)
r2 = 0.0314; % surface tension (kg/s^2)
rho1 = 1.05e-10; % bottom layer mass density (kg/nm^3)
rho2 = 1.05e-10; % top layer mass density (kg/nm^3)

Susceptibility Map

Logrithmic intensity map of

No = 206; % # of p (laplace frequency: Hz) or w (frequency) points
Nk = 50; % # of k (wave vector: nm^-1) points
kmax = 0.1; kmin = 0.001; k = 10.^linspace(log10(kmin),log10(kmax),Nk);
omax = 5; omin = -19; o = linspace(omin,omax,No);
chiType = [2,2]; % specify the chi type (chi11,chi12,chi21,or chi22)
inv_chi = zeros(No,Nk);
for iNk = 1:Nk
 for iNo = 1:No
 chi = sscptker(o(iNo),k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2);
 inv_chi(iNo,iNk) = 1/chi(chiType(1),chiType(2));
 end
end
inv_chi = log(abs(inv_chi));
figure('position',[300 200 500 200]); imagesc(log10(k),o,inv_chi);
set(gca,'xscale','linear','ydir','normal','xlim',[log10(kmin),log10(kmax)],'ylim',[omin,omax]);
ylabel('-i\omega (Hz)'); xlabel('log(k) (nm^{-1})'); colorbar;

Over-damped Relaxation Time Constants

chiType = [2,2]; % specify the chi type (chi11,chi12,chi21,or chi22)
Nk = 40; kmax = 0.12; kmin = 0.001; k = 10.^linspace(log10(kmin),log10(kmax),Nk);
numdiv = 100; % total number of points in p space
xmax = -1e2; xmin = .1; xp = linspace(xmin,xmax,numdiv); % p list
xpStep = abs(xp(1)-xp(2))/1000; % step size of p to search for poles
% predefine poles
% 1st col: k (wave vector)
% 2nd col: the root or zero point of 1/real(chi22)
% 3rd col: slope at this zero point
% 4th col: total intensity of power spectrum
poles = zeros(4,1000);
numpole = 0; % counter for number of poles
for iNk = 1:Nk % calculate tau at each k
 fp = zeros(1,numdiv-1); % initialize fp=1/chi
 for iNumdiv = 1:numdiv-1 % evaluate at fp at each p
 chi = sscptker(xp(iNumdiv),k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2);
 fp(iNumdiv) = 1./real(chi(chiType(1),chiType(2)));
 end
 for iNumdiv = 1:numdiv-2 % at each k, look for root of fp
 if fp(iNumdiv) * fp(iNumdiv+1) <=0
 try
 xprange = ...
 [min(xp(iNumdiv),xp(iNumdiv+1)),max(xp(iNumdiv),xp(iNumdiv+1))];
 pzero = fzero(@fzerosscpt,...
 xprange,[],k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2,chiType);
 numpole = numpole + 1;
 poles(2,numpole) = pzero;
 poles(1,numpole) = k(iNk);
 chiLeft = sscptker(...
 pzero-xpStep,k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2);
 chiRight = sscptker(...
 pzero+xpStep,k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2);
 slope = (1./real(chiRight(chiType(1),chiType(2))) ...
  -1./real(chiLeft(chiType(1),chiType(2))))/2/xpStep;
 poles(3,numpole) = 1.0/slope;
 chi = sscptker(1e-5,k(iNk),d,h1,eta1,eta2,r1,r2,rho1,rho2);
 poles(4,numpole) = 1.0/real(chi(chiType(1),chiType(2)));
  catch
 end
 end
 end
end
poles(:,find(poles(2,:)>=0)) = []; % remove zero and positive poles

Mode amplitude

figure('position',[300 200 500 200]);
plot(poles(1,:), abs( -poles(4,:).*poles(3,:)./poles(2,:) ),'-o');
set(gca,'Xscale','log','Xlim',[kmin,0.1],'YLim',[0.8 1.4]);
xlabel('k (nm^{-1})'); ylabel('Amplitude');

Plot Relaxation Time Constants

Both approximated analytical result and exact numerical result are shown for comparison. Analytical expressions do not have explicit dependence on mass densities. If interfacial tension is non zero, there is another slow mode. However, for the current PS system with zero interfacial tension, such a slow mode does not exist. It can be seen that the approximated analytical result agrees very well with the numerical calculations.

figure('position',[300,200,400,300]); hold on;
plot(poles(1,:),-1.0./poles(2,:),'ro'); % numerical result
tau0_fast = tau_fast(k,eta1,eta2,r1,r2,h1,d); % analytical result (fast mode)
% tau0_slow = tau_slow(k,eta1,eta2,r1,r2,h1,d); % analytical result (slow mode)
hold on; plot(k,tau0_fast,'b-');
set(gca,'Yscale','log','Xscale','log','Xlim',[0.001,0.1]); hold off; box on;
xlabel('k (nm^{-1})'); ylabel('\tau (sec)');
legend('numerical','analytical');

Last Updated: 02/09/06