diff options
| -rw-r--r-- | liouville.m | 140 | ||||
| -rw-r--r-- | useful_functions.m | 113 |
2 files changed, 154 insertions, 99 deletions
diff --git a/liouville.m b/liouville.m index bc9d524..b232a56 100644 --- a/liouville.m +++ b/liouville.m @@ -3,9 +3,7 @@ useful_functions; % some physical constants -hbar=1; -im_one=0+1i; - +useful_constants; % load atom energy levels and decay description three_levels; @@ -16,113 +14,57 @@ field_description; Nfreq=length(modulation_freq); -% now we create Liouville indexes list -% we unwrap density matrix and assign all posible -% frequencies as well -% resulting vector should be Nlevels x Nlevels x length(modulation_freq) -N=length(modulation_freq)*Nlevels*Nlevels; -rhoLiouville=zeros(N,1); -rhoLiouville_w=rhoLiouville; -rhoLiouville_r=rhoLiouville; -rhoLiouville_c=rhoLiouville; -i=0; -for w=1:length(modulation_freq) - for r=1:Nlevels - for c=1:Nlevels - i+=1; - rhoLiouville(i)=0; - rhoLiouville_w(i)=w; - rhoLiouville_r(i)=r; - rhoLiouville_c(i)=c; - endfor - endfor -endfor -% Liouville operator matrix -L=zeros(N); % NxN matrix -Li=zeros(N); % NxN Liouville interactive -L0=zeros(N); % NxN Liouville from unperturbed hamiltonian - -for p=1:N - for s=1:N - j=rhoLiouville_r(p); - k=rhoLiouville_c(p); - m=rhoLiouville_r(s); - n=rhoLiouville_c(s); - % we garanted to know frequency of final and initial rhoLiouville - w1i=rhoLiouville_w(p); - w2i=rhoLiouville_w(s); - w_jk=modulation_freq(w1i); - 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 wrequency in the list of available frequencyes - % but since we not garanteed to find it lets assign temporary 0 to Liouville matrix element - L(p,s)=0; - decay_part=0; - Lt=0; - for w3i=1:Nfreq - w_iner=modulation_freq(w3i); - decay_part=0; - if ((w_iner == w_l)) - %such frequency exist in the list of modulation frequencies - if ((w_iner == 0)) - L0=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m); - decay_part=\ - ( 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) ; - Lt=L0; - else - Li= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i); - Lt=Li; - endif - %Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero - % no need for above kron_delta since the same conditon checked in the outer if statement - Lt=-im_one/hbar*Lt - decay_part; - endif - endfor - if ((p == s)) - Lt+=-im_one*w_jk; - endif - L(p,s)=Lt; - endfor -endfor +%tune probe frequency +detuning_p=0; +N_detun_steps=15; +detuning_p_min=-15; +detuning_p_max=-detuning_p_min; +detuning_freq=zeros(1,N_detun_steps+1); +kappa_p =zeros(1,N_detun_steps+1); +detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps; +for detuning_p_cntr=1:N_detun_steps+1; +wp0=w12; +detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1); +wp=wp0+detuning_p; +modulation_freq=[0, wp, wd, -wp, -wd, wp-wd, wd-wp]; -% 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 +% now we create Liouville indexes list +[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq); +rhoLiouville=zeros(N,1); + +% Liouville operator matrix construction +L=Liouville_operator_matrix( + N, + H0, g_decay, g_dephasing, dipole_elements, + E_field, + modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c + ); -rhoLiouville_dot=rhoLiouville*0; -% 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; +%use the fact that sum(rho_ii)=1 to constrain solution +[rhoLiouville_dot, L]=constran_rho_and_match_L( + N, L, + modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c); %solving for density matrix vector rhoLiouville=L\rhoLiouville_dot; -rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq) -rho_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq) -rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq) +rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq); +rho_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq); +rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq); %rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq) -%L*rhoLiouville +kappa_p(detuning_p_cntr)=sum(sum(rho_1)); +detuning_freq(detuning_p_cntr)=detuning_p; + +%kappa_p_re=real(kappa_p); +%kappa_p_im=imag(kappa_p); + +endfor +figure(1); plot(detuning_freq, real(kappa_p)); +figure(2); plot(detuning_freq, imag(kappa_p)); diff --git a/useful_functions.m b/useful_functions.m index 748a677..5522afb 100644 --- a/useful_functions.m +++ b/useful_functions.m @@ -28,3 +28,116 @@ function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq) 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 + + +% Liouville operator matrix construction +function L=Liouville_operator_matrix( + N, + H0, g_decay, g_dephasing, dipole_elements, + E_field, + modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c + ) + %------------------------- + useful_constants; + L=zeros(N); % NxN matrix + + Nfreq=length(modulation_freq); + + for p=1:N + for s=1:N + j=rhoLiouville_r(p); + k=rhoLiouville_c(p); + m=rhoLiouville_r(s); + n=rhoLiouville_c(s); + % we garanted to know frequency of final and initial rhoLiouville + w1i=rhoLiouville_w(p); + w2i=rhoLiouville_w(s); + w_jk=modulation_freq(w1i); + 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 wrequency in the list of available frequencyes + % but since we not garanteed to find it lets assign temporary 0 to Liouville matrix element + L(p,s)=0; + decay_part=0; + Lt=0; + for w3i=1:Nfreq + w_iner=modulation_freq(w3i); + decay_part=0; + if ((w_iner == w_l)) + %such frequency exist in the list of modulation frequencies + if ((w_iner == 0)) + % calculate unperturbed part (Hamiltonian without EM field) + L0=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m); + decay_part=\ + ( 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) ; + Lt=L0; + else + % calculate perturbed part (Hamiltonian with EM field) + % in othe word interactive part of Hamiltonian + Li= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i); + Lt=Li; + endif + %Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero + % no need for above kron_delta since the same conditon checked in the outer if statement + Lt=-im_one/hbar*Lt - decay_part; + endif + endfor + if ((p == s)) + Lt+=-im_one*w_jk; + endif + L(p,s)=Lt; + 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]=constran_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 + |