1; % calculate total decay for particular level taking in account all branches function ret=decay_total(g_decay,i) ret=0; for k=1:size(g_decay)(1) ret=ret+g_decay(i,k); endfor endfunction % kroneker delta symbol function ret=kron_delta(i,j) if ((i==j)) ret=1; else ret=0; endif endfunction % this function create from Liouville density vector % the density matrix with given modulation frequency function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq) rho=zeros(Nlevels); for r=1:Nlevels for c=1:Nlevels rho(r,c)=rhoLiouville((freqIndex-1)*Nlevels^2+(r-1)*Nlevels+c); endfor endfor endfunction % we unwrap density matrix to Liouville density vector and assign all possible % modulation frequencies as well % resulting vector should be Nlevels x Nlevels x length(modulation_freq) function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq) N=Nfreq*Nlevels*Nlevels; rhoLiouville_w=zeros(N,1); rhoLiouville_r=zeros(N,1); rhoLiouville_c=zeros(N,1); i=0; for w=1:Nfreq; for r=1:Nlevels for c=1:Nlevels i+=1; rhoLiouville(i)=0; rhoLiouville_w(i)=w; % hold frequency modulation index rhoLiouville_r(i)=r; % hold row value of rho_rc rhoLiouville_c(i)=c; % hold column value of rho_rc endfor endfor endfor endfunction % create (Nlevels*Nlevels)x*(Nlevels*Nlevels) % sub matrices of Liouville operator % which repeat themselves for each modulation frequency % based on recipe from Eugeniy Mikhailov thesis function [L0m, polarizability_m]=L0_and_polarization_submatrices( ... rho_size, ... H0, g_decay, g_dephasing, dipole_elements, ... E_field, ... modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... ) %------------------------- useful_constants; % note that L0 and decay parts depend only on combination of indexes % jk,mn but repeats itself for every frequency L0m=zeros(rho_size); % (NxN)x(NxN) matrix decay_part_m=zeros(rho_size); % (NxN)x(NxN) matrix % polarization matrix will be multiplied by field amplitude letter % polarization is part of perturbation part of Hamiltonian polarizability_m=zeros(rho_size); % (NxN)x(NxN) matrix for p=1:rho_size % p= j*Nlevels+k % this might speed up stuff since less matrix passed back and force j=rhoLiouville_r(p); k=rhoLiouville_c(p); for s=1:rho_size % s= m*Nlevels+n m=rhoLiouville_r(s); n=rhoLiouville_c(s); % calculate unperturbed part (Hamiltonian without EM field) L0m(p,s)=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m); decay_part_m(p,s)= ... ( ... decay_total(g_decay,k)/2 ... + decay_total(g_decay,j)/2 ... + g_dephasing(j,k) ... )* kron_delta(j,m)*kron_delta(k,n) ... - kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ; polarizability_m(p,s)= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) ); endfor endfor L0m=-im_one/hbar*L0m - decay_part_m; endfunction % Liouville operator matrix construction % based on recipe from Eugeniy Mikhailov thesis function L=Liouville_operator_matrix( ... N, ... L0m, polarizability_m, ... E_field, ... modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... ) %------------------------- useful_constants; L=zeros(N); % NxN matrix Nfreq=length(modulation_freq); % Lets be supper smart and speed up L matrix construction % since it has a lot of voids. % By creation of rhoLiouville we know that there are % consequent chunks of rho_ij modulated with same frequency % this means that rhoLiouville is split in to Nfreq chunks % with length Nlevels*Nlevels=N/Nfreq rho_size=N/Nfreq; % Liouville matrix operator has Nlevels*Nlevels blocks % which governed by the same modulation frequency for p_freq_cntr=1:Nfreq for s_freq_cntr=1:Nfreq p0=1+(p_freq_cntr-1)*rho_size; s0=1+(s_freq_cntr-1)*rho_size; % we guaranteed to know frequency of final and initial rhoLiouville w1i=rhoLiouville_w(p0); % final w_jk=modulation_freq(w1i); w2i=rhoLiouville_w(s0); % initial w_mn=modulation_freq(w2i); % thus we know L matrix element frequency which we need to match w_l=w_jk-w_mn; % lets search this frequency in the list of available frequencies % but since we not guaranteed to find it lets assign temporary 0 to Liouville matrix element for w3i=1:Nfreq % at most we should have only one matching frequency w_iner=modulation_freq(w3i); if ((w_iner == w_l)) % yey, requested modulation frequency exist % lets do L sub matrix filling %such frequency exist in the list of modulation frequencies if ((w_iner == 0)) % calculate unperturbed part (Hamiltonian without EM field) L(p0:p0+rho_size-1,s0:s0+rho_size-1)=L0m; else % calculate perturbed part (Hamiltonian with EM field) % in other word interactive part of Hamiltonian L(p0:p0+rho_size-1,s0:s0+rho_size-1)= ... -im_one/hbar*polarizability_m* E_field(w3i); endif % diagonal elements are self modulated % due to rotating wave approximation if ((p0 == s0)) L(p0:p0+rho_size-1,s0:s0+rho_size-1)+= -im_one*w_jk*eye(rho_size); endif endif endfor endfor endfor endfunction % now generally rhoL_dot=0=L*rhoL has infinite number of solutions % since we always can resclale rho vector with arbitrary constant % lets constrain our density matrix with some physical meaning % sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1 % we will replace first row of Liouville operator with this condition % thus rhoLiouville_dot(1)=1 function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ... N, L, ... modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c) for i=1:N w2i=rhoLiouville_w(i); m=rhoLiouville_r(i); n=rhoLiouville_c(i); w=modulation_freq(w2i); if ((w==0) & (m==n)) L(1,i)=1; else L(1,i)=0; endif endfor rhoLiouville_dot= zeros(N,1); % sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1 % we will replace first row of Liouville operator with this condition % thus rhoLiouville_dot(1)=1 rhoLiouville_dot(1)=1; endfunction % calculate sucseptibility for the field at given frequency index function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq) rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq); kappa=0; for i=1:Nlevels for j=1:Nlevels kappa+=dipole_elements(j,i)*rho(i,j); endfor endfor endfunction % vim: ts=2:sw=2:fdm=indent