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
|
1;
clear all;
t0 = clock (); % we will use this latter to calculate elapsed time
% load useful functions;
useful_functions;
% some physical constants
useful_constants;
% load atom energy levels and decay description
%four_levels;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nlevels=4;
w1=1e9;
w2=0;
w_hpf=6800;
w3=w_hpf;
w4=w_hpf+.1; % separation of levels |3> and |4> somewhat like Zeeman splitting
w12=w1-w2;
w13=w1-w3;
%
% ----------- |1>
% / \
% E_d / \
% / \ E_p
% / \
% -------- |3> \
% -------- |4> \
% \
% ___________ |2>
% unperturbed Hamiltonian energy levels
levels_energy=[ w1, 0, w3, w4];
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_23=.001;
g_decay=zeros(Nlevels);
g_decay(1,2)=gamma; %upper level decay
g_decay(1,3)=gamma; %upper level decay
g_decay(1,4)=gamma; %upper level decay
g_decay(3,2)=gamma_23; % lower levels mixing
g_decay(2,3)=gamma_23; % lower levels mixing
g_decay(4,2)=gamma_23; % lower levels mixing
g_decay(2,4)=gamma_23; % 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);
g_dephasing(1,4)=g_deph;
g_dephasing(4,1)=g_dephasing(4,1);
% dipole matrix
dipole_elements=zeros(Nlevels);
dipole_elements(1,2)=1;
dipole_elements(2,1)=dipole_elements(1,2);
dipole_elements(1,3)=1;
dipole_elements(3,1)=dipole_elements(1,3);
dipole_elements(1,4)=1;
dipole_elements(4,1)=dipole_elements(1,4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EM field definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ep=0.01; %probe
Epc=conj(Ep);
Ed=.1; %drive
Edc=conj(Ed);
Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
Emc=conj(Em);
%wd=w13;
%wp=w12;
%wm=wd-(wp-wd);
%modulation_freq=[0, wp, wd, wm, -wp, -wd, -wm, wp-wd, wd-wp];
E_field =[0, Ep, Ed, Em, Epc, Edc, Emc, 0, 0 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Nfreq=length(modulation_freq);
%tune probe frequency
detuning_p=0;
N_detun_steps=100;
detuning_p_min=-1;
detuning_p_max=-detuning_p_min;
detuning_freq=zeros(1,N_detun_steps+1);
kappa_p =zeros(1,N_detun_steps+1);
kappa_m =zeros(1,N_detun_steps+1);
w_pf1=1e9;
w_hpf_ground=6800;
wd=w_pf1-w_hpf_ground;
wp0=w_pf1;
wp=wp0+detuning_p_min;
wm=wd-(wp-wd);
%modulation_freq=[0, wp, wd, wm, -wp, -wd, -wm, wp-wd, wd-wp];
%E_field =[0, Ep, Ed, Em, Epc, Edc, Emc, 0, 0 ];
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
[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
rhoLiouville=zeros(N,1);
dipoles=dipole_elements;
Efld=E_field;
clear dipole_elements;
dipole_elements.right=dipoles;
dipole_elements.left=dipoles;
dipole_elements.linear=dipoles;
clear E_field;
E_field.right=Efld;
E_field.left=Efld;
E_field.linear=Efld;
% calculate E_field independent properties of athe atom
% to be used as sub matrix templates for Liouville operator matrix
t1 = clock (); % we will use this latter to calculate elapsed time
[L0m, polarizability_m]=L0_and_polarization_submatrices( ...
Nlevels, ...
H0, g_decay, g_dephasing, dipole_elements ...
);
elapsed_time = etime (clock (), t1);
fprintf (stderr, "elapsed time for polarazability creation is %.3f sec\n",elapsed_time);
fflush (stderr);
t1 = clock (); % we will use this latter to calculate elapsed time
% 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);
elapsed_time = etime (clock (), t1);
fprintf (stderr, "elapsed time for full L creation is %.3f sec\n",elapsed_time);
fflush (stderr);
t1 = clock (); % we will use this latter to calculate elapsed time
%solving for density matrix vector
rhoLiouville=L\rhoLiouville_dot;
elapsed_time = etime (clock (), t1);
fprintf (stderr, "elapsed time for rhoL solving is %.3f sec\n",elapsed_time);
fflush (stderr);
rho_2=rhoOfFreq(rhoLiouville, 2, Nlevels);
rho_1=rhoOfFreq(rhoLiouville, 1, Nlevels);
L_new=L;
rhoLiouville_new=rhoLiouville;
rho_2_new=rho_2;
rho_1_new=rho_1;
% uncomment to update reference file
% save 'L_and_rhoL_referenced.mat' L rhoLiouville rho_1 rho_2
clear L rhoLiouville rho_2 rho_1;
load 'L_and_rhoL_referenced.mat'
diff_with_reference_L=sum(sum(abs(L_new-L)))
diff_with_reference_rhoL=(sum(abs(rhoLiouville_new -rhoLiouville)))
diff_with_rho_1=(sum(sum(abs(rho_1_new -rho_1))))
diff_with_rho_2=(sum(sum(abs(rho_2_new -rho_2))))
elapsed_time = etime (clock (), t0)
|
