summaryrefslogtreecommitdiff
path: root/useful_functions.m
blob: a1135f786fa9d03bba3a35f4ea03631d8a4d351a (plain)
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
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	


% 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
function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
		rho_size, ...
		H0, g_decay, g_dephasing, dipole_elements, ...
		E_field, ...
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
		)
	%-------------------------
	useful_constants;
	% note that L0 and decay parts depend only on combination of indexes
	% jk,mn but repeats itself for every frequency
	L0m=zeros(rho_size); % (NxN)x(NxN) 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

% Liouville operator matrix construction
% based on recipe from Eugeniy Mikhailov thesis
function L=Liouville_operator_matrix( 
		N, 
		L0m, polarizability_m,
		E_field, 
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
		)
	%-------------------------
	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


% 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]=constrain_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

% vim: ts=2:sw=2:fdm=indent