2
Quantum Project Optics Group 3
Hayden McGuinness, Anthony Clark, Anika Panzer 16 May 2007
1
1.1
Problem 1: Simulation of Quantum Diusion
Theory
The Schrodinger equation for a two level atom driven by a Weiner process, which simulates homodyne detector of the X1 quadrature, is given by 1 i + + + d| = H | dt 2 4
2
| dt +
1 + 2
| dW.
(1)
The equivalent master equation can be found by noting that, to second order, the equation for the density matrix can be written as d = d(| |) = d(| ) | + | d( |) + d(| )d( |). (2)
The compontents of d, namely d|( ) and d( |) are either given by equation (1) or the conjugate of (1). The d(| )d( |) term of (2) simplies greatly due to the fact that dt2 is taken to be zero while from Itos rule dW 2 = dt. Therefore
d(| )d( |) =
1 2
+
1 2
| |
1 2
+
2
dt dt
(3)
1 = 2 +
1 + + 4 +
Adding this term to to the d(| ) | and | d( |) terms (straight forward algebra that wont be shown here) yields the result i d = [H, ]dt + D[ ] + HH [ ]dW where 1 D[ ] = ( + ) 2 (5) (4)
and HH [ ] = + + . (6)
If ei then the equation for the density matrix would change accordingly, which leads to the expression i d = [H, ]dt + D[ei ] + HH [ei ]dW (7)
1.2
Data and results
The simulations were made for dierent time steps and a varying number of trajectories. The code can be found on the link on our webpage or under /projects/project2/problem1. The les that have actually been changed are sdesample.f90 and sderk support.f90. In sderk support.f90, we inserted our dierential equations in component form, while we changed the input, the output , the stepnumber and the stepsize. The parameters , and are also dened in this le. The outputles are ntraj1.txt, ntraj20.txt, ntraj100.txt, ntraj1000.txt where each output is 5 collumns wide contatining time, numerical excited state population, total population, excited state population of the corresponding Bloch equation and dierence between the exact and the analytic solution in this order. To judge how accurate numerical solutions to equation (1) are compared to available analytic solutions. There is no analytic solution to this dierential equation, but for an average over a large number of trajectories the solution should converge to the one for the optical Bloch equation. This behaviour will be investigated in the following. The system is driven exactly on resonance ( =0), so we expect the excited state population to converge to ee = for a large number of trajetcories. For the numerical solution initially all population was in the ground state, was set to 10, a time step size of .01 (where both quantities are normalized by , set to 1.0 here) was used, and the total time of evolution was 3s. First, the step size was varied to check whether we were in the right regime: a small variation of the stepsize shouldnt make a big dierence. This was done for our chosen step size .0.01, the results for tstep=0.005 match those of .01 within the desired precision, so our stepsize is of the right magnitude. For a large number of trajectories, the numerical result should converge to the analytic solution of the optical bloch equations. To study the dependence of the convergence on the number of trajectories, the program was varied such that it averaged over a given number of trajectories, dened by the user, after each time step. The number of trajectories was chosen as 1, 20, 100 and 10000. The convergence of the numerical result is well visible in the plots: with an increasing number of trajectories, the numerical result gets closer to the analytic solution of the OBEs. The graphs correspond to our expectations: we can see damped Rabi oscillations with noise that is due to the stochastic Wiener process. The comparison to the exact analytic solution to the optical bloch equations shows the convergence of the error for increasing number of trajectories. 2 22
3 2 3 (1 e 4 t (cos( t) + sin( t). + 2 4
(8)
1 ntraj1 using 1:2
0.8
0.6
0.4
0.2
0 0 0.5 1 1.5 2 2.5 3
Figure 1: Excited state population for one trajectory, = 10, = 0
1 ntraj20 using 1:2 ntraj20 using 1:5
0.8
0.6
0.4
0.2
0
-0.2 0 0.5 1 1.5 2 2.5 3
Figure 2: Excited state population for 20 trajectories, = 10, = 0, the upper line is the excited state solution, while the lower line is the dierence between the OBE solution and the numerical result
1 ntraj100 using 1:2 ntraj100 using 1:5
0.8
0.6
0.4
0.2
0
-0.2 0 0.5 1 1.5 2 2.5 3
Figure 3: Excited state population for 100 trajectories, = 10, = 0, the upper line is the excited state solution, while the lower line is the dierence between the OBE solution and the numerical result
1.3
Numerical Error
In the previous part of the problem, we studied the convergence of the numerical solution to the OBEs for a varying number of trajectories. Now, the number of trajectories will be kept one and we want to investigate the error dependence on time and stepsize. To nd out about the error convergence, we consider the function 3
1 ntraj10000 using 1:2 ntraj10000 using 1:5
0.8
0.6
0.4
0.2
0
-0.2 0 0.5 1 1.5 2 2.5 3
Figure 4: Excited state population for 10000 trajectories, = 10, = 0, the upper line is the excited state solution, while the lower line is the dierence between the OBE solution and the numerical result
f (t) f (t) = e(t, t) = e (t)(t) + O(t=0.5 ),
(9)
where f (t) is the numerical result to an SDE solver and f(t) is the exact solution to the corresponding dierential equation. In order to see the convergence of the error in dependence on the stepsize, we need a large numer of data at dierent timesteps. The command runscript is used to obtain this. It changes the number of substeps for each integration from 1 up to 512 and puts the data in outn. The collumns in the output les are again t, ee analytic, total population, ee exact and dierence between analytic and exact result. The input was again = 0, =1, =10, while the stepsize was set to .01 and the number of steps to 300. then the number of substeps was raised by the command runscript from 1 to 512. These les are read with an octave command called testscript, that takes the dierence between the SDE integrator solution and the exact solution (collumns 2 and 4) and divides it by the analytic solution. The graph shows error over stepsize, and it is visible that the error rises quickly after going through a minimum at a stepsize of 0.01. Our solution wil therefore be the closest to the exact solution for a stepsize around .01, which was also our initial choice. This result is reasonable as for a larger stepsize we expect a less exact result in general. The line is the best t which does notconnect our data points very well as their behaviour is not linear. This is due to the lack of an analytic result to this dierential equation; we cannot see the expected to the order of 1.5 behaviour of the error.
1 line 2
0.1 1e-05
0.0001
0.001
0.01
0.1
Figure 5: Error over timestep on a logarithmic scale, it is visiblt that the error has a dip slightly before .01 and inceases fast after that, the crosses are the measurement points while the line is the best t
4
To see that the integrator works ne, we check the testscript command on the original code. In this plot the behaviour proportional to t1.5 becomes visible.
00.001 .0001 line 2
0.0001 1e-05
1e-05 1e-06 1e-06 1e-07 1e-07
1e-08 1e-08
Figure 6: Error over timestep on a logarithmic scale, it is visiblt that the error scales with t1.5
1e-09 1e-05
0.0001
0.001
0.01
0.1
2
2.1
Problem 2: Numerical Evolution of a Quantum PDE
Introduction
In this problem we investigate the numerical solutions for the time evolution of a Gaussian wavepacket placed in a quantum harmonic oscillator, and compare these to the analytically obtainable results.
2.2
Split Operator Method
The split operator method for solving the Schrodinger equation is a means for time evolving an initial wavefunction for a known system Hamiltonian. In the case of the Harmonic oscilator with unit mass, the Hamiltonian, 1 p2 + 2 x2 = T (p) + V (x), 2 2
H=
(10)
contains no explicit time dependence, and so the time evolution operator is given by U (t1 , t0 ) = e
iH t
(11)
where t = t1 t0 . Now we make use of the fact that the kinetic term of the Hamiltonian has only momentum dependence and the potential term has only position dependence to rewrite this as a product of exponential operators. Mathematically we want Qn that solve the equation: eT +V =
n
eQn
(12)
5
The Qn that solve this are not unique, in fact it turns out that there are innitely many possible decompositions. The rst decomposition we consider is the tbree-term decomposition e(A+B ) = eA eB eZ , (13)
Expressing the exponential operators on both sides by their Taylor series expansions and expanding the operator Z in powers of , it is straightforward to equate terms of like order in from both sides to make the decomposition valid for any . n (T + V )n = n! n=1 Z=
n
m,l,q =0
m+l Am B l Z q m!l!q !
(14)
n Zn + AB Z2 Z0 = 0 Z1 = 0 Z2 = [A, B ]
(15)
0 : 1 : 2 :
1 2 2 (A
1 = 1 Z0 A + B = A + B Z1 + AB + BA + B 2 ) =
A2 2
(16)
+
B2 2
The lowest order nonzero Zn we get is Z2 = [A, B ], so we see that the lowest nontrivial order of in the series expansion of eZ will be 2 . Then letting A T, B V, and t gives the decomposition for the time evolution operator: e(T +V ) = eT eV + O((t)2 ), Neglecting the eZ operator introduces error of order O(t)3 . To reduce errors to O(t)3 , we can use a decomposition of form: e(A+B ) = e 2 eA e 2 eZ ,
B B
(17)
(18)
Writing the exponentials on both sides in terms of their Taylor series expansions, expanding Z as above and equating powers of gives, 0 : 1 : 2 : 3 :
1 2 2 (A
1 = 1 Z0 A + B = A + B Z1 + AB + BA + B 2 ) = 1 (A2 + AB + BA + B 2 ) Z2 2 + A2 + ABA + AB 2 + BA2 + BAB + B 2 A + B 3 ) +
A2 B 8
Z0 = 0 Z1 = 0 Z2 = 0 (19)
1 3 6 (A
=
A3 6
+
AB 2 4
+
B2 A 4
+
BA2 8
+
ABA 8
+
B2 6
Z3
BAB 6
Z3 =
A 2 B 24
+
AB 2 12
+
B2 A 12
6
BA2 24
+
ABA 12
+
so we see that the error in ignoring the eZ operator is O((t)2 ). It is apparrent that this decompositions yielding higher order error terms can be found by using further symetrized products of exponential operators. The higher order splittings will require smaller time steps for equivilent accuracy, thus lowering the computational resources required, but will take more resources to carry out the additional exponential operators required to carry out each time step. The optimal splitting/timestep combination is therefore found by balancing out these two factors. The Schrodinger equation for a particle in a potential well was numerically evolved for dierent step sizes and for orders two, four and six in the Richardson expansion. The initial wavepacket in the position basis was a Gaussian with the form =
2 2 1 e((xx0) /(4x )) 2x
(20)
where x0 = 0.3 and x = 0.4. There were 128 grid points on which the value of the wavepacket was calculated. The dimensionless time steps in which the global error was calculated were 0.00125, 0.0025, 0.005 and 0.01 for order methods two, four and six. Since the exact solution is know to be periodic with period T = 2/ (in these cases = 2 so T = 1) error calculated was by comparing the numerically obtained solution after 10 cycles to the initial conditions. Global error was calculated for several dierent parameters, expectation value of position, mean error of the wavepacket and for the single largest position value of the wavepacket (i.e. the error in the value of the top of the wavepacket, at x = 0.3). The only change in the given code was a change in the precision of the output from eight digits to 13 digits. This was necessary due to the order six method being accurate beyond eight signicant gures. Since this change is very minor the code will not be presented here. A simple MatLab script was designed to calculate errors and plot the results.
4
Log Log Plot of Error in <X> vs. Step Size
order 2 order 4 order 6
10 10 10 10
5
6
7
Log of error
10 10 10 10 10 10
8
9
10
11
12
13
10
2
Log of step size
Figure 7: Log Log plot of error in the expectation value of position calculated as a function of step size for orders two (blue), four (green) and six (red).
7
10 10 10 10
5
Log Log Plot of Mean X Error vs. Step Size
order 2 order 4 order 6
6
7
8
Log of error
10 10 10 10 10 10
9
10
11
12
13
14
10
2
Log of step size
Figure 8: Log Log plot of error in the mean wavepacket value calculated as a function of step size for orders two (blue), four (green) and six (red).
10 10 10 10
4
Log Log Plot of Single X Point Error vs. Step Size
order 2 order 4 order 6
5
6
7
Log of error
10 10 10 10 10 10
8
9
10
11
12
13
10
2
Log of step size
Figure 9: Log Log plot of error of a single position point calculated as a function of step size for orders two (blue), four (green) and six (red). The error results are shown in loglog plot style for expectation value, mean wavepacket and single point in gures 7,8 and 9 respectively. As would be expected for relatively large times steps, orders two and four show a very linearly relationship between log(error) and log(step size). This suggests very strongly that error as a function of step size is given by a power law relation. For order six the error does not exhibit this relationship, but can be understood from the premise that there is an optimal step size where the error is 8
least. This is due to for too small of step sizes round o error becomes a signicant factor while for too large of step sizes there is the more familiar error of the numerical approximation becoming less accurate. So for the appropriate range of step sizes this behavior is also expected. What is not expected is the power of the power law behavior for orders two and four in the parameter errors. For order two the powers are 4.000398, 4.000315 and 4.000258 for global errors in expectation value, mean wavepacket and single point respectively. For order four the powers are 5.018858, 5.050588 and 5.3242811 for global errors in expectation value, mean wavepacket and single point respectively. From theory, the expected powers should be two and four for orders two and four respectively, at least for one of the parameter errors.That is why three dierent parameter errors were calculated. It might be the case that one expects, for example, expectation value error to scale dierently from that of mean wavepacket value or of single position value, but this seems not to be the case. Simple problems, such as not having enough data or not running the simulation for long enough (more than 10 cycles) seems inadequate due to the clear trend in the data. From our knowledge these results are not explainable.2
3
Problem3: Quantum feedback control of the motion of a quantummechanical particle
In this part of the project we model the feedback control of a particle in a one-dimensional harmonic oscillator where the position of the particle is continuously measured. The stochastic master equation for such a system is given by i d = [H, ]dt + 2k D[x]dt + 2k H[x]dW with the Hamiltonian from problem 2 in the form p2 1 + 2 x2 . 2 2 (21)
H=
(22)
We use the easier stochastic Schrodinger equation projected into the position basis, so the evolution of the wave function is described by i d| >= m| > dt k (x x )2 | > dt + 2k (x x )| > dW in our program. To implement the stochastic process into the solution of the partial dierential equation, the second order operator splitting method implemented above was modied that the spatial parts of the time evolution were generated by the sderk routine. First, we used the SDE integartor to let the state evolve from t to t+t, using the ordinary dierential equation obtained by discarding the momentum term in the Hamiltonian. Then the t routine was used to take the fourier transform of the function, the split operator solver routine was used to apply the drift operator. Then the inverse fourier transform was taken and the time step was
1 This value is signicantly dierent from the others most likely due to the error of the smallest step size entering the region of optimal step size, where error due to numerical approximation becomes comparable to roundo error instead of dominating the error. 2 Unfortunately, we are running out of time on this problem set and must move on the more dicult third problem.
(23)
9
completed using the SDE integrator as before. This of course is then repeated until the input nal time is reached. The measurement process has the eect of heating the particle up. In order to counteract this heating, we want use the state information obtained from measurements, in particular the centroid position x at each time step, to change the potential the particle is in in such a way that the particle is slowed down. The parameter that we have control over in an such an experiment is the strength of the oscilator potential well. Noting that the process of measuring position tends to localize the particle in space, it is reasonable to rely primarily on the centroid position to build the feedback signal, as if the particle were classical, disregarding the energy due to the spatial spread of the particles wavefunction. In this case to cool the atom we need to adjust such that it is small when x is large and large when x is small resulting in eective cooling of the particle. This method of switching high and small suddenly is called the bang- bang method. Its goal is to apply kicks to the particle such that its environmental interaction is disrupted. It maximizes the rate of energy extraction. In order to cool the particle we have to extract energy from the system at a higher rate than the system feeds it in. To do this in a real homodyne experiment, we must rst determine x . The quantity which is actually measured is an electrical signal generated by a photodetector. This can be modeled as dW dy (t) = x dt + 8k
(2...