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
|
1;
% calculate total decay for particular level taking in account all branches
function ret=decay_total(g_decay,i)
ret=0;
for k=1:size(g_decay)(1)
ret=ret+g_decay(i,k);
endfor
endfunction
% kroneker delta symbol
function ret=kron_delta(i,j)
if ((i==j))
ret=1;
else
ret=0;
endif
endfunction
% 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
% we 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)
function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
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
% Liouville operator matrix construction
function L=Liouville_operator_matrix(
N,
H0, g_decay, g_dephasing, dipole_elements,
E_field,
modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
)
%-------------------------
useful_constants;
L=zeros(N); % NxN matrix
Nfreq=length(modulation_freq);
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;
decay_part=0;
Lt=0;
for w3i=1:Nfreq
w_iner=modulation_freq(w3i);
decay_part=0;
if ((w_iner == w_l))
%such frequency exist in the list of modulation frequencies
if ((w_iner == 0))
% calculate unperturbed part (Hamiltonian without EM field)
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) ;
Lt=L0;
else
% calculate perturbed part (Hamiltonian with EM field)
% in othe word interactive part of Hamiltonian
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 - decay_part;
endif
endfor
if ((p == s))
Lt+=-im_one*w_jk;
endif
L(p,s)=Lt;
endfor
endfor
endfunction
% 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
function [rhoLiouville_dot, L]=constran_rho_and_match_L(
N, L,
modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
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
% calculate sucseptibility for the field at given frequency index
function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)
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
|
