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
|
<?xml version="1.0"?>
<simulation xmds-version="2">
<name>Shahriar_system</name>
<author>Eugeniy Mikhailov, Simon Rochester</author>
<description>
License GPL.
Solving 3 level atom in double drive configuration
after Shahriar paper about white cavity
with field propagation along spatial axis Z
no Doppler broadening.
All fields detuned from upper level i.e. Raman configuration
*
* .....
* / ....
* / .... \
* / / \
* / /-------- |3>
* E3 / \
* / E2 \
* / / \ E1
* ------ |2> \
* \
* ------- |1>
*
We are solving
dE/dz+(1/c)*dE/dt=i*eta*rho_ij, where j level is higher then i.
Note that E is actually a Rabi frequency of electromagnetic field not the EM field
in xmds terms it looks like
dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space
VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
normalize it to something else look drho/dt equation.
No need to renormalizes eta as long as its express through i
the upper level decay rate in the same units as Rabi frequency.
</description>
<features>
<globals>
<![CDATA[
const double pi = M_PI;
const double c=3.e8;
const double lambda=794.7e-9; //wavelength in m
const double N=1e10*(1e6); //number of particles per cubic m i.e. density
const double Gamma_super=6*(2*M_PI*1e6); // characteristic decay rate of upper level used for eta calculations expressed in [1/s]
const double eta = 3*lambda*lambda*N*Gamma_super/8.0/M_PI; // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
// repopulation rate (atoms flying in/out the laser beam) in [1/s]
const double gt=0.01/2 *(2*M_PI*1e6);
// Natural linewidth of the upper level [1/s]
const double G=2*6 *(2*M_PI*1e6);
// total decay of i-th level branching ratios. Rij branching of i-th level to j-th
const double R31=0.5, R32=0.5;
complex E1c, E2c, E3c; // Complex conjugated Rabi frequencies
complex r21, r31, r32; // density matrix elements
]]>
</globals>
<benchmark />
<arguments>
<!-- Rabi frequency divided by 2 in [1/s] -->
<!--probe-->
<argument name="E1o" type="real" default_value="0.0025*(2*M_PI*1e6)" />
<!--pump fields-->
<argument name="E2o" type="real" default_value="1.0*(2*M_PI*1e6)" />
<argument name="E3o" type="real" default_value="1.0*(2*M_PI*1e6)" />
<!-- Fields detuning in [1/s] -->
<!-- probe field detuning-->
<argument name="d1" type="real" default_value="12*(2*M_PI*1e6)" />
<!-- averaged detuning of pump fields i.e. mid point -->
<argument name="da" type="real" default_value="12*(2*M_PI*1e6)" />
<!-- detuning of pump fields with respect to each other -->
<argument name="delta" type="real" default_value="6*(2*M_PI*1e6)" />
<!-- incoherent pumping rate from level |1> to |3> in [1/s]-->
<argument name="gp" type="real" default_value="2*2.0*(2*M_PI*1e6)" />
</arguments>
<bing />
<fftw plan="patient" />
<openmp />
<auto_vectorise />
</features>
<!-- 'z' and 't' to have dimensions [m] and [s] -->
<geometry>
<propagation_dimension> z </propagation_dimension>
<transverse_dimensions>
<dimension name="t" lattice="1000" domain="(-2.0e-6, 4.0e-6)" />
</transverse_dimensions>
</geometry>
<!-- Rabi frequency -->
<vector name="E_field" type="complex" initial_space="t">
<components>E1 E2 E3</components>
<initialisation>
<![CDATA[
// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
// as well as time
E1=E1o*exp(-pow( ((t-0.0)/1e-6),2) );
E2=E2o;
E3=E3o;
]]>
</initialisation>
</vector>
<vector name="density_matrix" type="complex" initial_space="t">
<components>r11 r22 r33 r12 r13 r23 </components>
<!--
note one of the level population is redundant since
r11+r22+r33=1
-->
<initialisation>
<![CDATA[
// Note:
// convergence is really slow if all populations concentrated at the bottom level |1>
// this is because if r11=1, everything else is 0 and then every small increment
// seems to be huge and adaptive solver makes smaller and smaller steps.
// As quick and dirty fix I reshuffle initial population
// so some of the population sits at the second ground level |2>
// TODO: Fix above. Make the equation of motion for r11
// and express other level, let's say r44
// through population normalization
r11 = 1; r22 = 0; r33 = 0;
r12 = 0; r13 = 0;
r23 = 0;
]]>
</initialisation>
</vector>
<vector name="pump_detunings" type="complex" initial_space="t">
<components> d dc </components>
<!--dc is probably redundant since it is just complex conjugate of d-->
<initialisation>
<![CDATA[
d = exp( i*t*delta );
dc = conj(d);
]]>
</initialisation>
</vector>
<sequence>
<!--For this set of conditions ARK45 is faster than ARK89-->
<integrate algorithm="ARK45" tolerance="1e-5" interval="4e-2">
<!--SIC algorithm seems to be much slower and needs fine 'z' step tuning and much finer time grid-->
<!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04-->
<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->
<samples>200 200</samples>
<operators>
<operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
<integration_vectors>density_matrix</integration_vectors>
<dependencies>E_field pump_detunings</dependencies>
<boundary_condition kind="left">
<![CDATA[
r11 = 1; r22 = 0; r33 = 0;
r12 = 0; r13 = 0;
r23 = 0;
]]>
</boundary_condition>
<![CDATA[
E1c = conj(E1);
E2c = conj(E2);
E3c = conj(E3);
r21=conj(r12);
r31=conj(r13);
r32=conj(r23);
// Equations of motions according to Simon's mathematica code
dr11_dt = gt - 2*(gp + gt)*r11 - E1*i*r13 + E1c*i*r31 + G*r33;
dr12_dt = (-gp - 2*gt - d1*i + da*i)*r12 - i*(E2*d + E3*dc)*r13 + E1c*i*r32;
dr13_dt = -(E1c*i*r11) - i*(E3c*d + E2c*dc)*r12 + (-G - gp - 2*gt - d1*i)*r13 + E1c*i*r33;
dr22_dt = gt - 2*gt*r22 - i*(E2*d + E3*dc)*r23 + i*(E3c*d + E2c*dc)*r32 + G*r33;
dr23_dt = -(E1c*i*r21) - i*(E3c*d + E2c*dc)*r22 + (-G - 2*gt - da*i)*r23 + i*(E3c*d + E2c*dc)*r33;
dr33_dt = 2*gp*r11 + E1*i*r13 + i*(E2*d + E3*dc)*r23 - E1c*i*r31 - i*(E3c*d + E2c*dc)*r32 - 2*(G + gt)*r33;
]]>
</operator>
<operator kind="ex" constant="yes">
<operator_names>Lt</operator_names>
<![CDATA[
Lt = i*1./c*kt;
]]>
</operator>
<integration_vectors>E_field</integration_vectors>
<dependencies>density_matrix</dependencies>
<![CDATA[
dE1_dz = i*eta*conj(r13) -Lt[E1] ;
dE2_dz = i*eta*conj(r23) -Lt[E2] ;
dE3_dz = i*eta*conj(r23) -Lt[E3] ;
]]>
</operators>
</integrate>
</sequence>
<!-- The output to generate -->
<output format="binary" filename="Shahriar_system.xsil">
<group>
<sampling basis="t(1000)" initial_sample="yes">
<dependencies>E_field</dependencies>
<moments>I1_out I2_out I3_out</moments>
<![CDATA[
I1_out = mod2(E1);
I2_out = mod2(E2);
I3_out = mod2(E3);
]]>
</sampling>
</group>
<group>
<sampling basis="t(100)" initial_sample="yes">
<dependencies>density_matrix</dependencies>
<moments>
r11_out r22_out r33_out
r12_re_out r12_im_out r13_re_out r13_im_out
r23_re_out r23_im_out
</moments>
<![CDATA[
// populations output
r11_out = r11.Re();
r22_out = r22.Re();
r33_out = r33.Re();
// coherences output
r12_re_out = r12.Re();
r12_im_out = r12.Im();
r13_re_out = r13.Re();
r13_im_out = r13.Im();
r23_re_out = r23.Re();
r23_im_out = r23.Im();
]]>
</sampling>
</group>
</output>
</simulation>
<!--
vim: ts=2 sw=2 foldmethod=indent:
-->
|
