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); rho(:)=rhoLiouville((freqIndex-1)*Nlevels^2+1:(freqIndex)*Nlevels^2); rho=rho.'; 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; rho_size = Nlevels*Nlevels; rhoLiouville_w=zeros(N,1); rhoLiouville_r=zeros(N,1); rhoLiouville_c=zeros(N,1); w=1:Nfreq; w_tmplate=(repmat(w,rho_size,1))(:); rhoLiouville_w=w_tmplate; r=1:Nlevels; r_tmplate=(repmat(r,Nlevels,1))(:); rhoLiouville_r=(repmat(r_tmplate,Nfreq,1))(:)'; c=(1:Nlevels)';% hold column value of rho_rc rhoLiouville_c=repmat(c,Nfreq*Nlevels,1); 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 % Matrices will be very sparse so lets create them in sparse representation L0m=sparse(zeros(rho_size)); % (Nlevels^2)x(Nlevels^2) matrix decay_part_m=sparse(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 = sparse(zeros(rho_size)); % (NxN)x(NxN) matrix polarizability_m.left = sparse(zeros(rho_size)); % (NxN)x(NxN) matrix polarizability_m.right = sparse(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); % speed up avoiding guaranteed to be zero calculations if ( kron_delta_m(k,n) && 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) ... ); endif if ( kron_delta_m(m,n) && kron_delta_m(j,k) ) decay_part_m(p,s) += - g_decay(m,j) ; endif % speed up avoiding guaranteed to be zero calculations if ( kron_delta_m(k,n) ) % calculate unperturbed part (Hamiltonian without EM field) L0m(p,s)=H0(j,m); polarizability_m.linear(p,s) = dipole_elements.linear(j,m); polarizability_m.left(p,s) = dipole_elements.left(j,m) ; polarizability_m.right(p,s) = dipole_elements.right(j,m) ; endif % speed up avoiding guaranteed to be zero calculations if ( kron_delta_m(j,m) ) % calculate unperturbed part (Hamiltonian without EM field) L0m(p,s)+=-H0(n,k)*kron_delta_m(j,m); polarizability_m.linear(p,s) +=-dipole_elements.linear(n,k); polarizability_m.left(p,s) +=-dipole_elements.left(n,k); polarizability_m.right(p,s) +=-dipole_elements.right(n,k); endif 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; % creating building blocks of L by rho_size * rho_size for w3i=1:Nfreq w_iner=modulation_freq(w3i); if ((w_iner == 0)) % calculate unperturbed part (Hamiltonian without EM field) L_sub{w3i}=L0m; else % calculate perturbed part (Hamiltonian with EM field) % in other word interactive part of Hamiltonian L_sub{w3i} = ... -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 endfor % Liouville matrix operator has Nlevels*Nlevels blocks % which governed by the same modulation frequency for p_freq_cntr=1:Nfreq p0=1+(p_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); for s_freq_cntr=1:Nfreq s0=1+(s_freq_cntr-1)*rho_size; 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 w3i=(w_l == modulation_freq); if (any(w3i)) % yey, requested modulation frequency exist % lets do L sub matrix filling % at most we should have only one matching frequency w_iner=modulation_freq(w3i); L(p0:p0+rho_size-1,s0:s0+rho_size-1) = L_sub{w3i}; endif endfor % diagonal elements are self modulated % due to rotating wave approximation L(p0:p0+rho_size-1,p0:p0+rho_size-1)+= -im_one*w_jk*eye(rho_size); 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, E_field) % 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; kappa.linear = sum( sum( transpose(dipole_elements.linear) .* rho ) ); kappa.left = sum( sum( transpose(dipole_elements.left) .* rho ) ); kappa.right = sum( sum( transpose(dipole_elements.right) .* rho ) ); kappa.linear /= E_field.linear( wi ); kappa.left /= E_field.left( wi ); kappa.right /= E_field.right( wi ); 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, E_field); xi_linear = xi.linear; xi_right = xi.right; xi_left = xi.left; endfunction function [dEdz_coef_linear, dEdz_coef_left, dEdz_coef_right] =dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field) % complex absorption coefficient % at given E_field frequency component Nlevels=( size(dipole_elements.linear)(1) ); rho=rhoOfFreq(rhoLiouville, wi, Nlevels); dEdz_coef_linear = sum( sum( transpose(dipole_elements.linear) .* rho ) ); dEdz_coef_left = sum( sum( transpose(dipole_elements.left) .* rho ) ); dEdz_coef_right = sum( sum( transpose(dipole_elements.right) .* rho ) ); endfunction function [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem) % complex absorption coefficient for each 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 ; Nfreq=length(modulation_freq); rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq); for wi=1:Nfreq [dEdz_linear(wi), dEdz_left(wi), dEdz_right(wi)] = dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field); endfor endfunction function relative_transmission=total_relative_transmission(atom_field_problem) % summed across all frequencies field absorption coefficient global atom_properties; [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem); %dEdz_total= dEdz_linear .* conj(dEdz_linear) ... %+dEdz_left .* conj(dEdz_left) ... %+dEdz_right .* conj(dEdz_right); %total_absorption = sum(dEdz_total); %total_absorption = |E|^2 - | E - dEdz | ^2 E_field.linear = atom_field_problem.E_field.linear; E_field.right = atom_field_problem.E_field.right; E_field.left = atom_field_problem.E_field.left; modulation_freq = atom_field_problem.modulation_freq ; pos_freq_indexes=(modulation_freq>0); % WARNING INTRODUCED HERE BUT MUST BE DEFINE IN ATOM PROPERTIES FILE %n_atoms - proportional to a number of interacting atoms n_atoms=500; transmited_intensities_vector = ... abs(E_field.linear + n_atoms*(1i)*dEdz_linear).^2 ... +abs(E_field.right + n_atoms*(1i)*dEdz_right).^2 ... +abs(E_field.left + n_atoms*(1i)*dEdz_left).^2; transmited_intensities_vector = transmited_intensities_vector(pos_freq_indexes); input_intensities_vector = ... abs(E_field.linear).^2 ... +abs(E_field.right).^2 ... +abs(E_field.left).^2; input_intensities_vector = input_intensities_vector(pos_freq_indexes); relative_transmission = sum(transmited_intensities_vector) / sum(input_intensities_vector); endfunction % create full list of atom modulation frequencies from positive light frequency amplitudes function [modulation_freq, E_field_amplitudes] = ... light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(... positive_light_frequencies, positive_light_field_amplitudes ) % we should add 0 frequency as first element of our frequencies list modulation_freq = cat(2, 0, positive_light_frequencies, -positive_light_frequencies); % negative frequencies have complex conjugated light fields amplitudes E_field_amplitudes.left = cat(2, 0, positive_light_field_amplitudes.left, conj(positive_light_field_amplitudes.left) ); E_field_amplitudes.right = cat(2, 0, positive_light_field_amplitudes.right, conj(positive_light_field_amplitudes.right) ); E_field_amplitudes.linear = cat(2, 0, positive_light_field_amplitudes.linear, conj(positive_light_field_amplitudes.linear) ); endfunction % vim: ts=2:sw=2:fdm=indent