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
|
<?xml version="1.0"?>
<simulation xmds-version="2">
<name>Genas_system</name>
<author>Eugeniy Mikhailov and Simon Rochester</author>
<description>
License GPL.
Solving 4 level atom in 0->1 configuration,
with field propagation along spatial axis Z
no Doppler broadening
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=1e6*(1e6); //number of particles per cubic m i.e. density
const double Gamma=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/16.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*M_PI*1e6);
// Natural linewidth of upper state in [1/s]
const double g0=1.*(2*M_PI*1e6);
complex Exc, Eyc; // Complex-conjugated Rabi frequency
complex r21, r31, r41, r32, r42, r43, r44; // density matrix elements
]]>
</globals>
<benchmark />
<arguments>
<!-- Real and imaginary parts of complex Rabi frequency in [1/s] -->
<argument name="ExReo" type="real" default_value="(3.+0.001)*(2*M_PI*1.e6)" />
<argument name="ExImo" type="real" default_value="0." />
<argument name="EyReo" type="real" default_value="0." />
<argument name="EyImo" type="real" default_value="(3.-0.001)*(2*M_PI*1.e6)" />
<!-- light detuning in [1/s] -->
<argument name="delta" type="real" default_value="3.0*(2*M_PI*1e6)" />
<!--shift of upper-state M=0 sublevel-->
<argument name="delta0" type="real" default_value="1.*(2*M_PI*1e6)" />
<!--Static B-field Larmor frequency-->
<argument name="OL" type="real" default_value="4.*(2*M_PI*1e6)" />
<!--rf Rabi frequency-->
<argument name="Orf" type="real" default_value="0.1*(2*M_PI*1e6)" />
<!--rf frequency-->
<argument name="orf" type="real" default_value="4.*(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="(-1.5e-7, 2.5e-7)" />
</transverse_dimensions>
</geometry>
<!-- Rabi frequency -->
<vector name="E_field" type="complex" initial_space="t">
<components>Ex Ey</components>
<initialisation>
<![CDATA[
// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
// as well as time
Ex=(ExReo+i*ExImo)*exp(-pow( ((t-0.0)/1e-7),2) );
Ey=(EyReo+i*EyImo)*exp(-pow( ((t-0.0)/1e-7),2) );
]]>
</initialisation>
</vector>
<vector name="density_matrix" type="complex" initial_space="t">
<components>r11 r22 r33 r12 r13 r14 r23 r24 r34 r44</components>
<initialisation>
<![CDATA[
r11 = 1; r22 = 0; r33 = 0; r44 = 0;
r12 = 0; r13 = 0; r14 = 0;
r23 = 0; r24 = 0;
r34 = 0;
]]>
</initialisation>
</vector>
<vector name="rfField" type="real">
<components> rf </components>
<initialisation>
<![CDATA[
rf = Orf*sin(orf*t);
]]>
</initialisation>
</vector>
<sequence>
<integrate algorithm="ARK45" tolerance="0.05e-7" interval="1e1">
<samples>200 200</samples>
<operators>
<operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
<integration_vectors>density_matrix</integration_vectors>
<dependencies>E_field rfField</dependencies>
<boundary_condition kind="left">
<![CDATA[
r11 = 1; r22 = 0; r33 = 0; r44 = 0;
r12 = 0; r13 = 0; r14 = 0;
r23 = 0; r24 = 0;
r34 = 0;
]]>
</boundary_condition>
<![CDATA[
Exc = conj(Ex);
Eyc = conj(Ey);
r21=conj(r12);
r31=conj(r13);
r41=conj(r14);
r32=conj(r23);
r42=conj(r24);
r43=conj(r34);
// Equations of motions according to Simon's mathematica code
dr11_dt = gt - gt*r11 + (-Ey - Ex*i)*r12 + i*(Ex + Ey*i)*r14 + i*(Exc + Eyc*i)*r21 + (-Eyc - Exc*i)*r41 + g0*(r22 + r33 + r44);
dr12_dt = (2*(Eyc - Exc*i)*r11 - i*((2*delta - (g0 + 2*gt)*i - 2*OL - 2*rf)*r12 - 2*(Exc + Eyc*i)*r22 + 2*(Exc - Eyc*i)*r42))/2.;
dr13_dt = (-g0/2. - gt - delta*i + delta0*i)*r13 + i*(Exc + Eyc*i)*r23 + (-Eyc - Exc*i)*r43;
dr14_dt = -(i*((2*delta - (g0 + 2*gt)*i + 2*OL + 2*rf)*r14 - 2*(Exc + Eyc*i)*r24 - 2*(Exc - Eyc*i)*(r11 - r44)))/2.;
dr22_dt = (Ey + Ex*i)*r12 + (Eyc - Exc*i)*r21 - (g0 + gt)*r22;
dr23_dt = (Ey + Ex*i)*r13 + i*(delta0 + i*(g0 + gt + i*OL) - rf)*r23;
dr24_dt = (Ey + Ex*i)*r14 + (Eyc + Exc*i)*r21 - (g0 + gt + 2*i*OL + 2*i*rf)*r24;
dr33_dt = -((g0 + gt)*r33);
dr34_dt = (Eyc + Exc*i)*r31 - i*(delta0 - (g0 + gt)*i + OL + rf)*r34;
dr44_dt = (Ey - Ex*i)*r14 + (Eyc + Exc*i)*r41 - (g0 + gt)*r44;
]]>
</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[
dEx_dz = i*eta*conj(r12-r14) - Lt[Ex] ;
dEy_dz = -eta*conj(r12+r14) - Lt[Ey] ;
]]>
</operators>
</integrate>
</sequence>
<!-- The output to generate -->
<output format="binary" filename="Genas_system.xsil">
<group>
<sampling basis="t(100)" initial_sample="yes">
<dependencies>E_field</dependencies>
<moments>Ex_re_out Ex_im_out Ey_re_out Ey_im_out</moments>
<![CDATA[
Ex_re_out = Ex.Re();
Ex_im_out = Ex.Im();
Ey_re_out = Ey.Re();
Ey_im_out = Ey.Im();
]]>
</sampling>
</group>
<group>
<sampling basis="t(100)" initial_sample="yes">
<dependencies>density_matrix</dependencies>
<moments>
r11_out r22_out r33_out r44_out
r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out
r23_re_out r23_im_out r24_re_out r24_im_out
r34_re_out r34_im_out
</moments>
<![CDATA[
// populations output
r11_out = r11.Re();
r22_out = r22.Re();
r33_out = r33.Re();
r44_out = r44.Re();
// coherences output
r12_re_out = r12.Re();
r12_im_out = r12.Im();
r13_re_out = r13.Re();
r13_im_out = r13.Im();
r14_re_out = r14.Re();
r14_im_out = r14.Im();
r23_re_out = r23.Re();
r23_im_out = r23.Im();
r24_re_out = r24.Re();
r24_im_out = r24.Im();
r34_re_out = r34.Re();
r34_im_out = r34.Im();
]]>
</sampling>
</group>
</output>
</simulation>
<!--
vim: ts=2 sw=2 foldmethod=indent:
-->
|
