summaryrefslogtreecommitdiff
path: root/xmds2/Shahriar_system/Shahriar_system.xmds
blob: 98da563b7ed2105b3cac138062a4bdb88bb7034b (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
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
<?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>
		*

		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		IMPORTANT simplification: E2 = E3 by magnitude then we can use a field consisting
		of beat note between E1 and E2 oscillating at frequency wa=(w1+w2)/2
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		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=2e10*(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, Efc; // 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="10000"   domain="(-2.0e-6, 4.0e-6)" />
		</transverse_dimensions>
	</geometry>

	<!-- Rabi frequency --> 
	<vector name="E_field" type="complex" initial_space="t">
		<components>E1 Ef</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) );
			// very dirty hack, I assume that E2=E3 then I can use beatnote formula since
			// w1 approcimately equals w2
			// I assign it to new combined field with amplitude Ef
			// E*cos( (wa+d/2)*t ) * E*cos( (wa-d/2)*t ) = 2*E*cos(d*t/2)*cos(wa*t) = Ef*cos(wa*t)
			Ef = 2*E2o*cos(delta*t/2);  // E3o assumed to be equal E2o

			]]>
		</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>

	<sequence>
		<!--For this set of conditions ARK45 is faster than ARK89-->
		<integrate algorithm="ARK45" tolerance="1e-5" interval="10e-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</dependencies>
          <boundary_condition kind="left">
            <![CDATA[
							r11 = 1; r22 = 0; r33 = 0; 
							r12 = 0; r13 = 0; 
							r23 = 0; 
            ]]>
          </boundary_condition>
					<![CDATA[
						E1c = conj(E1);
						Efc = conj(Ef);

						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 - Ef*i*r13 + E1c*i*r32;
						dr13_dt = -(E1c*i*r11) - Efc*i*r12 + (-G - gp - 2*gt - d1*i)*r13 + E1c*i*r33;
						dr22_dt = gt - 2*gt*r22 - Ef*i*r23 + Efc*i*r32 + G*r33;
						dr23_dt = -(E1c*i*r21) - Efc*i*r22 + (-G - 2*gt - da*i)*r23 + Efc*i*r33;
						dr33_dt = 2*gp*r11 + E1*i*r13 + Ef*i*r23 - E1c*i*r31 - Efc*i*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] ;
					dEf_dz = i*eta*conj(r23) -Lt[Ef] ;
          ]]>
			</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 If_out</moments>
				<![CDATA[
				I1_out = mod2(E1);
				If_out = mod2(Ef);
				]]>
			</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: 
-->