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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
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
|
