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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
1;
function ret=decay_total(g_decay,i)
% calculate total decay for particular level taking in account all branches
ret=0;
for k=1:size(g_decay)(1)
ret=ret+g_decay(i,k);
endfor
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, Nfreq)
% this function create from Liouville density vector
% the density matrix with given modulation frequency
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
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;
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
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);
% note that L0 and decay parts depend only on combination of indexes
% jk,mn but repeats itself for every frequency
L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) matrix
decay_part_m=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=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);
% calculate unperturbed part (Hamiltonian without EM field)
L0m(p,s)=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m);
decay_part_m(p,s)= ...
( ...
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) ;
polarizability_m(p,s)= ( dipole_elements(j,m)*kron_delta(k,n)-dipole_elements(n,k)*kron_delta(j,m) );
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;
% Liouville matrix operator has Nlevels*Nlevels blocks
% which governed by the same modulation frequency
for p_freq_cntr=1:Nfreq
for s_freq_cntr=1:Nfreq
p0=1+(p_freq_cntr-1)*rho_size;
s0=1+(s_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);
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
for w3i=1:Nfreq
% at most we should have only one matching frequency
w_iner=modulation_freq(w3i);
if ((w_iner == w_l))
% yey, requested modulation frequency exist
% lets do L sub matrix filling
%such frequency exist in the list of modulation frequencies
if ((w_iner == 0))
% calculate unperturbed part (Hamiltonian without EM field)
L(p0:p0+rho_size-1,s0:s0+rho_size-1)=L0m;
else
% calculate perturbed part (Hamiltonian with EM field)
% in other word interactive part of Hamiltonian
L(p0:p0+rho_size-1,s0:s0+rho_size-1)= ...
-im_one/hbar*polarizability_m* E_field(w3i);
endif
% diagonal elements are self modulated
% due to rotating wave approximation
if ((p0 == s0))
L(p0:p0+rho_size-1,s0:s0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
endif
endif
endfor
endfor
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, Nlevels, Nfreq)
% calculate susceptibility for the field at given frequency index
rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq);
kappa=0;
for i=1:Nlevels
for j=1:Nlevels
kappa+=dipole_elements(j,i)*rho(i,j);
endfor
endfor
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
% vim: ts=2:sw=2:fdm=indent
|
