1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
% some physical constants
hbar=1;
im_one=0+1i;
%three_levels;
two_levels;
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
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)
if ((i==j))
ret=1;
else
ret=0;
endif
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
% no need for above kron_delta since the same conditon checked in the outer if statement
Lt=-im_one/hbar*Lt;
endif
endfor
if ((p == s))
Lt+=-im_one*w_jk;
endif
L(p,s)=Lt;
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;
else
L(1,i)=0;
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_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq)
rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq)
%rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq)
%L*rhoLiouville
|
