Optimal control of atom transport for quantum gates in optical lattices
Abstract
By means of optimal control techniques we model and optimize the manipulation of the external quantum state (centerofmass motion) of atoms trapped in adjustable optical potentials. We consider in detail the cases of both non interacting and interacting atoms moving between neighboring sites in a lattice of a doublewell optical potentials. Such a lattice can perform interactionmediated entanglement of atom pairs and can realize twoqubit quantum gates. The optimized control sequences for the optical potential allow transport faster and with significantly larger fidelity than is possible with processes based on adiabatic transport.
pacs:
03.67.a, 34.50.s,I Introduction
Quantum degenerate gases, such as BoseEinstein condensates (BECs) blochreview () or cold Fermi gases ingusciofermi (), trapped in optical lattices, provide a flexible platform for investigating condensed matter models and quantum phase transitions blochqpt (). It has been proposed to use these systems as quantum simulators of solid state systems ciracchains () and for implementing quantum information processing (QIP) brennen (); jaksch (); dorner (). Experiments on neutral atoms have shown some of the ingredients needed for QIP: the preparation of a Mott insulator state with just one particle per well, which is used as the initial state of a quantum register blochqpt (), singlequbit rotation lee07 (), and controlled motion of atoms so as to effect entangling interactions mandelNature03 (); lee07 ().
A general requirement of QIP is accurate control of a quantum system. Often this includes control of degrees of freedom other than the qubit or computational basis, for example the center of mass motion of an ion or atom for which the spin (internal state) represents the qubit. One approach to achieving such accurate control is adiabatic manipulation of the relevant Hamiltonian. Unfortunately adiabaticity limits the speed of operations. One way to overcome this difficulty is to use optimal control methods oc1 (); dorner (). Here we show that such techniques could improve the speed and fidelity of transport of atoms in an optical lattice.
Recent experiments anderliniSwap (); mandelNature03 (); trotzky08a () have shown that quantum gates could be implemented in controllable optical potentials by adjusting the overlap between atoms trapped in neighboring sites of an optical lattice. High fidelity of this dynamic process could be achieved by adiabatically changing the trapping potential. This, however, limits the overall gate speed to be much lower than the trapping frequency dorner (); garcia2003 (). Here we present a detailed numerical analysis of the transport process used to effect a twoqubit quantum gate in anderliniSwap (), which is performed with the controllable doublewell optical potential described in anderlinipra () and find that it gives an accurate description of the evolution measured in the experiment. Then we apply optimal control theory to the transport process of the atoms, both with and without interactions, to show how to increase the speed of the gate. The success of this method in this specific problem demonstrates the promise of optimal control for coherent manipulation of a diverse class of quantum systems.
Ii The experiment
A twoqubit quantum gate with neutral atoms can be realized in optical lattices through a controlled interactioninduced evolution of the wavefunction that depends on the states of the two atoms brennen (); jaksch (). Because atoms in their electronic ground states generally have shortrange interactions, in order to use these contact interactions to produce entanglement, the atomic wavefunctions must be made to overlap. Once the interaction has taken place for a fixed time, the two atoms can be separated again thus finishing the gate. In this paper we consider the control of such motion in a specific setup; however our theory can be applied to more general systems.
ii.1 The DoubleWell Lattice
Neutral Rb atoms are loaded into the sites of a 3D optical lattice obtained by superimposing a 2D optical lattice of doublewells anderlinipra () in the horizontal plane and an independent 1D optical lattice in the vertical direction. The horizontal lattice has a unit cell that can be dynamically transformed between singlewell and doublewell configurations. The horizontal potential experienced by the atoms is anderliniJPhysB ():
(1)  
where and are the spatial coordinates, is the laser wavevector and is the laser wavelength. The potential (1) depends on three parameters:

the strength of the potential wells;

the ratio of vertical to horizontal electric field components;

the phase shift between vertical and horizontal light components.
The angle determines the height of the barrier between adjacent doublewell sites: by changing from 0 to the potential changes from a symmetric doublewell configuration, with a spacing of lambda/2 (lambda/2 lattice), to a singlewell configuration, with a spacing of lambda (lambda lattice). By changing and together one varies the energy offset (tilt) of a well with respect to the neighboring one. The tilt of the double well is zero for = , while it is maximum (with a value depending on ) for = 0 or .
To effect a quantum gate, one varies the three parameters in time so as to move atoms occupying adjacent wells into the same well, allowing them to interact and finally returning them to their original positions.
In Fig. 1 we show four snapshots of the crosssection of the optical potential along the direction of the double wells (), and of the singleparticle wave functions of the two atoms during a particular transport sequence. In the initial configuration each atom is prepared in the ground state of separate wells so that the properly symmetrized initial state is:
(2) 
where and are wavefunctions localized in the left and right well, respectively, and 1 and 2 are the labels of the two (indistinguishable) atoms (see Fig. 1a). For quantum gate operation, we would also include the internal state of the atoms, but here we concentrate only on the external state (centerofmass motion). and are linear combinations of the lowest symmetric and antisymmetric energy eigenfunctions of the singleparticle potential. In Eq. (2) and throughout the text we use the convention that single (two) particle states are labeled with lower (upper) case greek letters.
The potential is changed by lowering the barrier and at the same time lowering the right well with respect to the left one, Fig. 1b)c). The atom initially in the right well remains in the ground state, evolving into the lowest state of the final potential, while the atom initially in the left well evolves into the first excited state , Fig. 1d). When the two atoms are in the same potential well they interact through the usual contact interaction, which can be used to generate the entangling operation needed to realize a twoqubit quantum gate jaksch (); dorner ().
ii.2 Experimental Procedure
The experimental characterization of the transport process is accomplished by performing the potential transformation depicted in Fig. 1 with atomic samples loaded either only in the left sites or only in the right sites of the double wells anderliniJPhysB ().
Briefly, BoseEinstein condensates of atoms with are loaded in the sites of the lattice with an exponential ramp of 200 ms duration. This loading populates only the ground band of the optical potential with mostly one atom per lattice site SebbyPRL07 (). Then the potential is transformed to the lattice in such a way that the atoms eventually occupy either only the right sites or only the left sites of the double wells anderliniJPhysB (); lee07 (). Starting from this initialized state, we perform the transport process illustrated in Fig. 1a)–d). At the end of the process we measure the occupation of the lattice bands. To this purpose, we map the quasimomentum of atoms occupying different vibrational levels of the optical potential onto real momenta lying within different Brillouin zones kastberg95 (); greiner2001 (). This is achieved by switching off the optical potential in 500 s and acquiring an absorption image of the sample after a 13 ms timeofflight. In this way atoms occupying different vibrational levels appear spatially separated, allowing us to measure the amount of population in each vibrational state.
The comparison between these measurements and the theoretical model requires an accurate determination of the evolution of the parameters , and characterizing the optical lattice during the experimental sequences. The parameter , which corresponds to the depth of the potential when it is set in the configuration, is measured by pulsing the lattice and observing the resulting momentum distribution in time of flight Ovchinnikov1998 (). The parameters and , which determine the shape of the doublewell potential and are controlled using two electrooptic modulators (EOMs), are determined from measurements of the polarization of the laser beams after the EOMs as a function of their respective control voltages anderlinipra ().
We perform two series of experimental sequences in order to study the properties of the atomic transport as a function of the duration of the process and of the energy tilt between the two potential wells during the merge. In a first series of measurements the sequence involves converting the lattice from the doublewell to the singlewell configuration by changing , rotating the polarization of the input light using a linear ramp, while leaving constant the light intensity and the setting of the electrooptic modulator EOM dedicated to the control of the phase shift . This sequence is repeated for different durations of the linear ramp from ms to 1.01 ms. In a second series of measurements we consider the dependence of the transport on the tilt of the doublewell potential during the merge. We perform the lattice transformation using always the same duration of ms, the same intensity of the light beam and the same ramp for changing the polarization angle , while the the setting of EOM is kept constant in time during a sequence. We then repeat the sequence for different settings of EOM. The time dependence for the three lattice parameters , and for measurements of the first series and of the second series are shown in Fig. 2a and Fig. 2b, respectively. Fig. 2b shows the evolution of the parameter for two different settings of EOM. The potential parameters are determined using our calibration of the lattice setup, taking into account effects such as different losses on the optical elements for different polarizations of the lattice beams and the dependence of the optical potential on the local polarization of the light lee07 (). These effects are responsible for the change of the potential depth and of the angle during the sequence despite the fact that both the intensity of the light and the settings of EOM are not actively changed.
Iii Theoretical model
Here we describe the theoretical methods that we implement for investigating the dynamics in the system described above, starting with the case of noninteracting particles. Then, we consider the experimental realization of the merging of adjacent lattice sites into a single site shown in Fig. 1 and we compare the results obtained in the experiment with the expectations based on our theoretical model. This stage represents a useful benchmark to evaluate the reliability of the numerical model as well as for gaining insight into the details of the optical potential experienced by the atoms. Finally, we present the technique for optimizing the transport sequence, and we show how we can achieve a significantly higher fidelity at fixed operation time for the atomic motion than by using smooth sequences based on adiabatic evolution.
iii.1 Theoretical Framework
We consider the 1D problem restricted to the x axis by assuming that the optical potential can be separated along the three spatial directions, allowing us to express the atomic wavefunctions as a product of three independent terms. We consider the harmonic approximation of the potential in the and directions, having trap frequencies and respectively, that can be calculated as shown in spielman06 () and we assume that along and the atoms always occupy the lowest vibrational state. This restriction does not put limitations in studying dynamic processes involving low energy states of the doublewell potential since it can be chosen to have nondegenerate vibration frequencies along all three directions, with the lowest frequency always along . We calculate the eigenstates of the system along the direction by solving the eigenvalue equation using the finite difference method thomas (). For the time evolution we consider the integration of the time dependent Schrödinger equation using the CrankNicolson method cn (). This method has the advantage of being unconditionally stable and the error in the results scale quadratically with the number of spacetime grid points in which the Schrödinger equation is solved. The relative error of the data presented is always less than . In Appendix A we present a more detailed description of our numerical methods.
iii.2 Comparison to experimental results
In this section we present the theoretical analysis of the transport processes described above and we discuss the agreement between the model and the experimental measurements. We start by considering the time evolution of the Hamiltonian spectrum during the two sequences a and b shown in Fig. 3.
At time the spectrum is made of nearly degenerate doublets of almost equally spaced pairs of harmonic oscillator states, while at time the levels are similar to those of a single harmonic oscillator^{1}^{1}1We have verified that the results for the 1D spectrum are in good agreement with full calculations in 2D (restricted to states with vibrational excitation along the direction). Additional energy levels are present in the 2D spectrum, associated with states with vibrational excitation along . However, those states can be neglected for studying the dynamical process considered here since their energy is always higher than the three lowest states of the 1D spectrum..
Fig. 3 shows the time evolution of the eigenstates of the singleparticle Hamiltonian. The atoms initially prepared in the two local ground states in the right and left wells ( and ) evolve into the instantaneous eigenstates ending in the ground and first excited state of the final configuration, respectively. This approach requires the sequence to be performed slowly with respect to the timescale associated with the gaps between the relevant energy levels. The optimal “speed” in the parameter space can be calculated using the LandauZener theory for avoided level crossings.
For gaining quantitative insight into the properties of the transport we perform numerical simulations for the sequences used in the experiments, also taking into account possible deviations of the parameters from the experimental calibrations, and we compare the results with the experimental measurements. The relevant quantities for our analysis will be the overlap of the energy eigenstates of the final potential with the evolved state where indicates the initial well occupation:
(3) 
where the operator is the singleparticle timeevolution operator from time to time . In the experiment can be measured as the population of each energy level at the end of the process.
We now consider how the atoms evolve when the parameters change according to sequence a of Fig. 2 as a function of the total time , focusing on atoms starting in ^{2}^{2}2Both in the experiments and in the simulations the evolution of the atom initially in the right well, i.e. in state , shows a weaker dependence on the properties of the sequence and is less instructive. For instance, in the simulations for ms the population in the ground state is of order of for a broad range of parameters..
In Fig. 4 we show the final population of the ground , first and second excited state measured in the experiments and calculated for four values of which differ from the one of Fig. 2 a by a constant offset ^{3}^{3}3We do not consider variations in and due to the small dependence of the transport process on those parameters within the range associated with the accuracy of our calibrations.. The results obtained by the model are in reasonable agreement with the experimental observations; we find best quantitative matching for = 0.03, for which the rms deviation between model and theory is reduced from 0.13(at ) to 0.08.
Now we consider the sequence b of Fig. 2, performed for different ending values of the parameter around . As shown in Fig. 5, both the experiment and the model show a strong dependence on for the transport of the atom starting in the left site of the double well. The transport into the first excited state has a maximum theoretical fidelity of 0.95 for = 0.474. Less negative values of , i.e. increasing tilts, lead to a decrease of fidelity due to the increase in the fraction of population ending in the second excited state. Values of closer to 0.5, i.e. more symmetric configurations of the double well, lead to decrease of fidelity associated with larger fractions of population ending in the ground state. Also in this case the experimental data and the theoretical model are in satisfactory agreement. For these data the deviation between theory and experiment is more sensitive to the value of the phase shift . We find best agreement by shifting the value determined from the calibration by an offset = 0.015, which reduces the rms deviation from 0.4 to 0.15. The axis for the experimental data in Fig. 5 has been corrected by the offset .
Thus, while showing the reliability of the model in describing the dynamic process, the comparison between theoretical and experimental results also allows one to refine the calibration of the parameters characterizing the optical potential.
Finally we find that adding an offset of = 0.016 to the calibration brings the data from both sequences to a good agreement with the theory and reduces the rms deviation from 0.19 to 0.11. This is three times larger than the uncertainty of the offset from our EOM calibration but is still consistent with measurements of the lattice tilt from SebbyPRL07 (). This discrepancy might be explained by the birefringence in the vacuum cell windows, which is not accounted for in our model. Inclusion of this offset should improve both the predictivity of the model and the experimental optimization of the collisional gate based on the numerical technique described below.
Iv Optimized transport
In this section we employ optimal control theory to obtain fast and highfidelity gates. Our aim is to find a temporal dependence of the control parameters , , that improves the fidelity even for a shorter sequence duration, when the adiabatic sequences presented above yield a poor fidelity. Quantum optimal control techniques have been successfully employed in a variety of fields: molecular dynamics oc1 (); oc2 (); oc3 (), dynamics of ultracold atoms in optical lattices tannor2002 (); vaucher (); romero (), implementation of quantum gates dorner (); montaoc ().
We use the Krotov algorithm krotov () as the optimization procedure. The objective is to find the optimal shapes of the control parameter sequences that maximize the overlap (fidelity) between the evolved initial wave function and a target wave function. The initial and target wave functions are fixed a priori. The algorithm works also for more than one particle. The method consists in iteratively modifying the shape of the control parameters according to a “steepest descent method” in the space of functions (for more details see dorner ()). The method requires evolving each particle’s wave function and an auxiliary wave function backward and forward in time according to the Schrödinger equations. In our simulations we use the CrankNicolson scheme to realize this step as described in Appendix A.
iv.1 Noninteracting case
We optimize the gate for ms^{4}^{4}4We chose this time in order to show the benefits of the optimization procedure for a sequence duration which cannot provide a good fidelity with smooth parameter ramps based on adiabatic evolution. While a shorter duration could be chosen in principle, this choice can be easily experimentally implemented with no major changes in the present experimental apparatus, allowing future prosecution of our studies on this subject. Increasing the total time should improve the best fidelity. choosing as a starting point for the optimization a sequence similar to Fig. 2b where is for simplicity taken constant to the final value . Without optimization the fidelities for the atom initially in the left and right well are and , respectively. The infidelities are shown in Fig. 6 as a function of the number of optimization steps: the algorithm of optimization is proven to yield a monotonic increase in fidelity oc1 (), however it does not guarantee to reach its 100% value. The results for the two atoms give a fidelity above .
The resulting optimized parameter sequences are shown in Fig. 7 and compared to the original sequence without optimization. We find that the optimized sequence for the potential depth differs negligibly from the initial guess. In principle, the algorithm could achieve still higher singleparticle fidelities from different starting points.
In Fig. 8 we show the square absolute value of the wave functions of the two atoms as a function of time, the 1D potential time dependence and the projections of the initially leftwell state onto the lowest four instantaneous energy eigenstates :
(4) 
Notice that . As can be easily seen, the optimal time evolution is much less smooth than the adiabatic one as it takes advantage of quantum interference between nonadiabatic excitation paths to obtain better results.
iv.2 Interaction effects
Up to now we have considered only the singleparticle evolution of the system, i.e. without including any interaction between the particles. This approximation is valid in our transport sequence as long as the two wave functions in nearby wells do not overlap. When the two particles overlap in the same well we must take into account interactions, and we model them with an effective 1D contact potential:
(5) 
where are the coordinates of the two atoms and is an effective coupling strength olshanii () expressed by , where is the scattering length for atoms and is the Planck constant. The spectrum is modified by the interactions: the state with one atom in each well is lower by kHz than the doublyoccupied states.
As in the case without interactions we start with each atom localized in a separate well. Notice that we are considering wavefunctions that are symmetric under the exchange of the coordinates of the two particles.
We consider the twoparticle fidelity:
(6) 
where is the twoparticle evolution operator for the Hamiltonian of the two atoms, which includes interactions. is an eigenstate of the twoparticle Hamiltonian at time , corresponding in the limit to the symmetrized product of the singleparticle wavefunctions localized in each well (see Eq.(2)); the target state is an eigenstate of the twoparticle Hamiltonian at time which, in the limit of vanishing interactions, corresponds to the state
(7) 
The square modulus of and , in the twoatom coordinate representation, are shown in Fig. 9ab. In order to make a comparison between the interacting and non interacting cases we define a twoparticle fidelity also in the noninteracting case:
(8) 
where and are the singleparticle evolution operators for the two atoms without interactions.
non optimized  0.69  0.57  0.22  0.22 
transport optimized  0.99  0.99  0.98  0.93 
interaction optimized  0.98  0.98  0.96  0.97 
In Table 1 we summarize our results for ms obtained with three different sequences: first, the non optimized sequence Fig. 2b; second, the transport optimized case Fig. 7 where we used the optimal control algorithm to optimize the singleparticle populations not taking into account interactions; third, the interaction optimized case where we apply the optimal control algorithm using as the initial guess the transport optimized sequence Fig. 7 and then optimizing including the interactions in the evolution.
The resulting wavefunctions for the non optimized and transport optimized sequences are compared in Fig. 9cd. Without optimal control the twoparticle fidelity with and without interactions is while with (noninteracting) optimization we obtain and . This shows that interactions spoil slightly the efficiency of the transport process as one might expect. Optimal control can subsequently be applied while including interactions in the optimization, producing a control sequence with a fidelity of .
Another consideration is the experimental bandwidth available for feedback control. The optimized control waveforms Fig. 7 were obtained with no restriction on the frequency response of the control, and typically have frequency components on the order of a few times the lattice vibrational spacings (see Fig. 10), i.e. larger than the bandwidth of our control electronics. Clearly, using a filtered version of these waveforms will lead to lower control fidelity and it will be important to increase the experimental bandwidth of the control electronics (currently about 50 kHz). In addition, it may be useful to develop an optimization sequence that includes the limited control bandwidth, although it is likely that frequencies on the order of the vibrational spacing will always be needed.
V Conclusions
We have presented a detailed, numerical analysis of the transport process of neutral atoms in a time dependent optical lattice. We show how to improve the fidelity of the transport process for ms from , using simple adiabatic switching, to , using optimal control theory. We expect better results for longer control times. We analyze the effect of atomatom interactions on the transport process and we show that the optimal control parameter sequences found in the noninteracting case still work when including interaction. We obtained the same transformation as in the case of the adiabatic transport with a better fidelity and in a time shorter by more than a factor of three, which represents a relevant improvement in terms of scalability of the number of gates that can be performed before the system decoheres due to the coupling to its environment. This technique can be easily adapted to other similar transport processes and also extended to atoms in different magnetic states, which can allow the implementation of a fast, highfidelity quantum gate in a real optical lattice setup with the qubits encoded in the atomic internal states anderliniSwap (). In the future, it would be interesting to study the possibility of including the effect of errors in the optimization procedure and thus investigate in more details the robustness and noiseresilience of the optimal control technique.
Acknowledgements.
This work was supported by the European Commission, through projects SCALA (FET/QIPC), EMALI (MRTNCT2006035369) and QOQIP, by the National Science Foundation through a grant for the Institute for Theoretical Atomic, Molecular and Optical Physics at Harvard University and Smithsonian Astrophysical Observatory, and by DTO. SM acknowledges support from ISTEUROSQIP and the Quantum Information program of “Centro De Giorgi” of Scuola Normale Superiore. The computations have been performed on the HPC facility of the Department of Physics, University of Trento. We thank J. SebbyStrabley for experimental support.Appendix A Numerical method
In our numerical simulations we employ a finite difference method (see for example thomas ()) that consists in discretizing the coordinate representation in a homogeneous points mesh in the interval : , , . The number is the lattice spacing. In this discretized representation the eigenvalue equation becomes:
(9) 
where , (3.5 is the conversion coefficient between kHz and lattice recoils) and is the discretized wave function. The discretized second order derivative operator acts on any function as:
(10) 
Expression (9) is second order in . If one arranges the function in an dimensional array, then Eq. (9) can be rewritten as an eigenvalue problem with a Hamiltonian matrix . Since this matrix is tridiagonal, the principal diagonal being and the two subdiagonals being filled with , one can take advantage of its sparse structure for storage and computing. We calculate low energy eigenstates by approximately diagonalizing the Hamiltonian using the JacobiDavidson method jacobi (). This method is an iterative method which is capable of finding a few (typically less than ) eigenstates close to a chosen target energy. The advantages of using this and similar algorithms (Lanczos, Arnoldi) instead of exact diagonalization are that the method is much faster and one does not need to store the whole Hamiltonian. In our eigenvalue problem we take full advantage of this method, given the sparse structure of the Hamiltonian . The error of this approximate method compared to an exact diagonalization method is negligible for our purposes.
To study the time evolution of the wave functions of the atoms we integrate numerically the time dependent Schrödinger equation. Introducing a time slicing in the interval with time interval , the Schrödinger equation has the form:
(11) 
The discretized expression (11) gives an iterative relation to compute the wavefunction at time , from the expression of the wavefunction at time . This is one example of an explicit method: the coefficients can be directly calculated from . Explicit schemes have the great advantages of being extremely fast and easily implemented. However this expansion is only first order in and is not always stable. Therefore, we used the CrankNicolson scheme cn (), an implicit method, which consists in taking a time average of the righthand side of (11) between time and , namely
This method is of the second order in time and space and it is unconditionally stable. The price for all these advantages is that a tridiagonal set of linear equations must be solved to get as shown in Eq. (A). We used common Fortran routines to solve the linear equations problem lapack ().
We solve a 2D time dependent Schrödinger equation in the two coordinates of the atoms by making use of the extension in two dimensions of the CrankNicolson method called the PeacemanRachford method thomas (). This is an implicit method and the integration proceeds in two steps: first the initial wavefunction is integrated in time considering only one direction in the coordinate space, then from the intermediate wavefunction we obtain the final wavefunction by integrating in the other direction. This method is an example of alternate direction implicit schemes.
In our simulations we used and that assures convergence of the results with a relative error which is less than .
References
 (1) For a review, see for example: I. Bloch, J. Phys. B 38, S629 (2005).
 (2) G. Modugno, F. Ferlaino, R. Heidemann, G. Roati, M. Inguscio, Phys. Rev. A 68, 011601(R) (2003).
 (3) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature London 415, 39 (2002).
 (4) J. J. GarcìaRipoll, M. A. MartinDelgado, and J. I. Cirac Phys. Rev. Lett. 93, 250405 (2004).
 (5) G. K. Brennen, C. M. Caves, P. S. Jessen, and I. H. Deutsch, Phys. Rev. Lett. 82, 1060 (1999).
 (6) D. Jaksch, H.J. Briegel, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 82, 1975 (1999).
 (7) T. Calarco, U. Dorner, P. Julienne, C. Williams, and P. Zoller, Phys. Rev. A 70, 012306 (2004).
 (8) P. J. Lee, M. Anderlini, B. L. Brown, J. SebbyStrabley, W. D. Phillips and J. V. Porto, Phys. Rev. Lett. 99, 020402 (2007).
 (9) O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. Hänsch and I. Bloch, Nature 425, 937 (2003).
 (10) A. P. Peirce, M. A. Dahleh, and H. Rabitz, Phys. Rev. A 37, 4950 (1988). R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni, and D. J. Tannor, Chem. Phys. 139, 201 (1989).
 (11) S. Trotzky, P. Cheinet, S.Folling, M. Feld, U. Schnorrberger, A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and I. Bloch, Science 319, 295 (2008).
 (12) M. Anderlini, P. J. Lee, B. L. Brown, J. SebbyStrabley, W. D. Phillips and J. V. Porto, Nature 448 452–456 (2007).
 (13) J. J. GarcìaRipoll, P. Zoller, and J. I. Cirac, Phys. Rev. Lett. 91, 157901 (2003).
 (14) J. SebbyStrabley, M. Anderlini, P. S. Jessen, and J. V. Porto, Phys. Rev. A 73, 033605 (2006).
 (15) M. Anderlini, J. SebbyStrabley, J. Kruse, J.V. Porto, and W.D. Phillips, J. Phys. B 39, S199 (2006).
 (16) J. SebbyStrabley, B. L. Brown, M. Anderlini, P. J. Lee, W.D. Phillips, J. V. Porto, and P. R. Johnson, Phys. Rev. Lett. 98, 200405 (2007).
 (17) A. Kastberg, W. D. Phillips, S. L. Rolston, R. J. C. Spreeuw, and P. S. Jessen, Phys. Rev. Lett. 74, 1542 (1995).
 (18) M. Greiner, I. Bloch, O. Mandel, T. W. Hänsch, and T. Esslinger, Phys. Rev. Lett. 87, 160405 (2001).
 (19) Yu. B. Ovchinnikov, J. H. Müller, M. R. Doery, E. J. D. Vredenbregt, K. Helmerson, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett 83, 284 (1999).
 (20) I. B. Spielman, P. R. Johnson, J. H. Huckans, C. D. Fertig, S. L. Rolston, W. D. Phillips, and J. V. Porto Phys. Rev. A 73, 020702(R) (2006).
 (21) J. W. Thomas: “Numerical Partial Differential Equations”, vol. 1, Springer, New York (1995).
 (22) J. Crank and P. Nicolson, Proc. Cambridge Philos. Soc. 43, 50 (1947).
 (23) S. A. Rice and M. Zhao: “Optical control of Molecular dynamics”, J. Wiley, New York (2000).
 (24) M. Shapiro and P. Brumer: “Principles of the quantum control of molecular processes”, J. Wiley, New York (2003).
 (25) S. E. Sklarz and D. J. Tannor, Phys. Rev. A 66, 053619 (2002).
 (26) B. Vaucher, S. R. Clark, U. Dorner, D. Jaksch, New J. Phys. 9 221 (2007).
 (27) O. RomeroIsart and J. J. GarcìaRipoll, Phys. Rev. A 76, 052304 (2007).
 (28) S. Montangero, T. Calarco and R. Fazio, Phys. Rev. Lett. 99, 170501 (2007).
 (29) V. F. Krotov: “Global methods in optimal control theory”, M. Dekker Inc., New York (1996).
 (30) M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).
 (31) G. L. G. Sleijpen and H. A. van der Vorst, SIAM Review 42, 267 (2000), we used the routine written by G. L. G. Sleijpen and coworkers: http://www.math.uu.nl/people/sleijpen/
 (32) E. Anderson et al.: “LAPACK User’s guide”, SIAM, Philadelphia (1999).