1; function ret=decay_total(g_decay,i) % calculate total decay for particular level taking in account all branches ret=sum(g_decay(i,:)); endfunction function ret=kron_delta(i,j) % kroneker delta symbol if ((i==j)) ret=1; else ret=0; endif endfunction function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels) % this function create from Liouville density vector % the density matrix with given modulation frequency 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 function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq) % 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) 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 function [L0m, polarizability_m]=L0_and_polarization_submatrices( ... Nlevels, ... H0, g_decay, g_dephasing, dipole_elements ... ) % 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 %------------------------- useful_constants; rho_size=Nlevels*Nlevels; % now we create Liouville indexes list [Ndummy, rhoLiouville_w_notused, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,1); kron_delta_m=eye(Nlevels); % note that L0 and decay parts depend only on combination of indexes % jk,mn but repeats itself for every frequency L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) 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.linear = zeros(rho_size); % (NxN)x(NxN) matrix polarizability_m.left = zeros(rho_size); % (NxN)x(NxN) matrix polarizability_m.right = 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_m(k,n)-H0(n,k)*kron_delta_m(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_m(j,m)*kron_delta_m(k,n) ... - kron_delta_m(m,n)*kron_delta_m(j,k)*g_decay(m,j) ; polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta_m(k,n)-dipole_elements.linear(n,k)*kron_delta_m(j,m) ); polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta_m(k,n)-dipole_elements.left(n,k)*kron_delta_m(j,m) ); polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta_m(k,n)-dipole_elements.right(n,k)*kron_delta_m(j,m) ); endfor endfor L0m=-im_one/hbar*L0m - decay_part_m; endfunction function L=Liouville_operator_matrix( ... N, ... L0m, polarizability_m, ... E_field, ... modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ... ) % Liouville operator matrix construction % based on recipe from Eugeniy Mikhailov thesis %------------------------- 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.linear * E_field.linear(w3i) ... -im_one/hbar*polarizability_m.left * E_field.left(w3i) ... -im_one/hbar*polarizability_m.right * E_field.right(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 function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ... N, L, ... modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c) % 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 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 function kappa=susceptibility(wi, rhoLiouville, dipole_elements) % calculate susceptibility for the field at given frequency index Nlevels=( size(dipole_elements.linear)(1) ); rho=rhoOfFreq(rhoLiouville, wi, Nlevels); kappa.linear=0; kappa.left=0; kappa.right=0; for i=1:Nlevels for j=1:Nlevels kappa.linear += dipole_elements.linear(j,i) * rho(i,j); kappa.left += dipole_elements.left(j,i) * rho(i,j); kappa.right += dipole_elements.right(j,i) * rho(i,j); endfor endfor endfunction function index=freq2index(freq, modulation_freq) % convert modulation freq to its index in the modulation_freq vector index=[1:length(modulation_freq)](modulation_freq==freq); endfunction function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq) % calculates rhoLiouville vector assuming steady state situation and normalization of rho_ii to 1 Nlevels=sqrt( size(L0m)(1) ); Nfreq=length(modulation_freq); % now we create Liouville indexes list [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq); % Liouville operator matrix construction L=Liouville_operator_matrix( N, L0m, polarizability_m, E_field, modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ); %use the fact that sum(rho_ii)=1 to constrain solution [rhoLiouville_dot, L]=constrain_rho_and_match_L( N, L, modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c); %solving for density matrix vector rhoLiouville=L\rhoLiouville_dot; endfunction function [xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( atom_field_problem) % find steady state susceptibility at particular modulation frequency element % at given E_field global atom_properties; L0m = atom_properties.L0m ; polarizability_m = atom_properties.polarizability_m ; dipole_elements = atom_properties.dipole_elements ; E_field = atom_field_problem.E_field ; modulation_freq = atom_field_problem.modulation_freq ; freq_index = atom_field_problem.freq_index ; rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq); xi=susceptibility(freq_index, rhoLiouville, dipole_elements); xi_linear = xi.linear; xi_right = xi.right; xi_left = xi.left; endfunction % vim: ts=2:sw=2:fdm=indent