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
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
|
1;
function ret=decay_total(g_decay,i)
% calculate total decay for particular level taking in account all branches
ret=sum(g_decay(i,:));
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)
% this function create from Liouville density vector
% the density matrix with given modulation frequency
rho=zeros(Nlevels);
rho(:)=rhoLiouville((freqIndex-1)*Nlevels^2+1:(freqIndex)*Nlevels^2);
rho=rho.';
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;
rho_size = Nlevels*Nlevels;
rhoLiouville_w=zeros(N,1);
rhoLiouville_r=zeros(N,1);
rhoLiouville_c=zeros(N,1);
w=1:Nfreq;
w_tmplate(:)=repmat(w,rho_size,1);
rhoLiouville_w=w_tmplate';
r=1:Nlevels;
r_tmplate(:)=repmat(r,Nlevels,1);
rhoLiouville_r(:)=repmat(r_tmplate',Nfreq,1);
c=(1:Nlevels)';% hold column value of rho_rc
rhoLiouville_c=repmat(c,Nfreq*Nlevels,1);
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);
kron_delta_m=eye(Nlevels);
% 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.linear = zeros(rho_size); % (NxN)x(NxN) matrix
polarizability_m.left = zeros(rho_size); % (NxN)x(NxN) matrix
polarizability_m.right = 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_m(k,n)-H0(n,k)*kron_delta_m(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_m(j,m)*kron_delta_m(k,n) ...
- kron_delta_m(m,n)*kron_delta_m(j,k)*g_decay(m,j) ;
polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta_m(k,n)-dipole_elements.linear(n,k)*kron_delta_m(j,m) );
polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta_m(k,n)-dipole_elements.left(n,k)*kron_delta_m(j,m) );
polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta_m(k,n)-dipole_elements.right(n,k)*kron_delta_m(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;
% creating building blocks of L by rho_size * rho_size
for w3i=1:Nfreq
w_iner=modulation_freq(w3i);
if ((w_iner == 0))
% calculate unperturbed part (Hamiltonian without EM field)
L_sub{w3i}=L0m;
else
% calculate perturbed part (Hamiltonian with EM field)
% in other word interactive part of Hamiltonian
L_sub{w3i} = ...
-im_one/hbar*polarizability_m.linear * E_field.linear(w3i) ...
-im_one/hbar*polarizability_m.left * E_field.left(w3i) ...
-im_one/hbar*polarizability_m.right * E_field.right(w3i) ...
;
endif
endfor
% Liouville matrix operator has Nlevels*Nlevels blocks
% which governed by the same modulation frequency
for p_freq_cntr=1:Nfreq
p0=1+(p_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);
for s_freq_cntr=1:Nfreq
s0=1+(s_freq_cntr-1)*rho_size;
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
w3i=(w_l == modulation_freq);
if (any(w3i))
% yey, requested modulation frequency exist
% lets do L sub matrix filling
% at most we should have only one matching frequency
w_iner=modulation_freq(w3i);
L(p0:p0+rho_size-1,s0:s0+rho_size-1) = L_sub{w3i};
endif
endfor
% diagonal elements are self modulated
% due to rotating wave approximation
L(p0:p0+rho_size-1,p0:p0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
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)
% calculate susceptibility for the field at given frequency index
Nlevels=( size(dipole_elements.linear)(1) );
rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
kappa.linear=0;
kappa.left=0;
kappa.right=0;
for i=1:Nlevels
for j=1:Nlevels
kappa.linear += dipole_elements.linear(j,i) * rho(i,j);
kappa.left += dipole_elements.left(j,i) * rho(i,j);
kappa.right += dipole_elements.right(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
function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq)
% calculates rhoLiouville vector assuming steady state situation and normalization of rho_ii to 1
Nlevels=sqrt( size(L0m)(1) );
Nfreq=length(modulation_freq);
% now we create Liouville indexes list
[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
% Liouville operator matrix construction
L=Liouville_operator_matrix(
N,
L0m, polarizability_m,
E_field,
modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
);
%use the fact that sum(rho_ii)=1 to constrain solution
[rhoLiouville_dot, L]=constrain_rho_and_match_L(
N, L,
modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
%solving for density matrix vector
rhoLiouville=L\rhoLiouville_dot;
endfunction
function [xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( atom_field_problem)
% find steady state susceptibility at particular modulation frequency element
% at given E_field
global atom_properties;
L0m = atom_properties.L0m ;
polarizability_m = atom_properties.polarizability_m ;
dipole_elements = atom_properties.dipole_elements ;
E_field = atom_field_problem.E_field ;
modulation_freq = atom_field_problem.modulation_freq ;
freq_index = atom_field_problem.freq_index ;
rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
xi=susceptibility(freq_index, rhoLiouville, dipole_elements);
xi_linear = xi.linear;
xi_right = xi.right;
xi_left = xi.left;
endfunction
% vim: ts=2:sw=2:fdm=indent
|
