diff options
Diffstat (limited to 'liouville.m')
| -rw-r--r-- | liouville.m | 190 |
1 files changed, 190 insertions, 0 deletions
diff --git a/liouville.m b/liouville.m new file mode 100644 index 0000000..baff45a --- /dev/null +++ b/liouville.m @@ -0,0 +1,190 @@ + +levels=1:3; + +% ----------- |1> +% / \ +% E_d / \ +% / \ E_p +% / \ +% -------- |3> \ +% \ +% ___________ |2> + + +% some physical constants +hbar=1; +im_one=0+1i; + +Nlevels=3; +w12=1e9; +w_hpf=6800; +w13=w12-w_hpf + + +% unperturbed Hamiltonian energy levels +levels_energy=[ w12, 0, w_hpf]; +levels_energy=levels_energy*hbar; +H0=zeros(Nlevels); +H0=diag(levels_energy); +%for i=1:Nlevels + %H0(i,i)=levels_energy(i); +%endfor + +% decay matrix g(i,j) correspnds to decay from i-->j +gamma=6; +gamma_13=.01 +g_decay=zeros(Nlevels); +g_decay(1,2)=gamma; %upper level decay +g_decay(1,3)=gamma; %upper level decay +g_decay(3,2)=gamma_13; % lower levels mixing +g_decay(2,3)=gamma_13; % lower levels mixing + +%defasing matris +g_deph=0; +g_dephasing=zeros(Nlevels); +g_dephasing(1,2)=g_deph; +g_dephasing(2,1)=g_dephasing(1,2); +g_dephasing(1,3)=g_deph; +g_dephasing(3,1)=g_dephasing(1,3); + + + +% dipole matrix +dipole_elements=zeros(Nlevels); +dipole_elements(1,2)=1; +dipole_elements(2,1)=dipole_elements(1,2); +dipole_elements(3,1)=1; +dipole_elements(1,3)=dipole_elements(3,1); + + +%EM field definition +Ep=100.0; +Epc=conj(Ep); +Ed=0; +Edc=conj(Ed); +wd=w13; +wp=w12; +modulation_freq=[0, wp, wd, -wp, -wd, wp-wd, wd-wp]; +E_field =[0, Ep, Ed, Epc, Edc, 0, 0 ]; +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) +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 + +\ + +N=length(rhoLiouville); +% Liouville operator matrix +L=zeros(N); % NxN matrix +Li=zeros(N); % NxN Liouville interactive +L0=zeros(N); % NxN Liouville from unperturbed hamiltonian + +function ret=decay_total(g_decay,i) + ret=0; + for k=1:size(g_decay)(1) + ret=ret+g_decay(i,k); + endfor +endfunction + +function ret=kron_delta(i,j) + ret=((i==j)); +endfunction + +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; + for w3i=1:length(modulation_freq) + w_iner=modulation_freq(w3i); + 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) ; + decay_part=decay_part*hbar/im_one; + L0=decay_part; + 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 + Lt+=-im_one*w_jk*kron_delta(w_iner,w_jk); + L(p,s)=Lt; + endif + endfor + endfor +endfor + + +% 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; + endif +endfor + +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; + + +%solving for density matrix vector +rhoLiouville=L\rhoLiouville_dot; + +% 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 + + +rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq) +%rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq) + +%L*rhoLiouville + |