Correcting errors in a quantum gate with pushed ions via optimal control
aa r X i v : . [ qu a n t - ph ] A ug Correcting errors in a quantum gate with pushed ions via optimal control
Uffe V. Poulsen
Lundbeck Foundation Theoretical Center for Quantum System ResearchDepartment of Physics and Astronomy, University of Aarhus, DK 8000 Aarhus C, Denmark
Shlomo Sklarz and David Tannor
Department of Chemical Physics, Weizmann Institute of Science, 76100 Rehovot , Israel
Tommaso Calarco
Institut f¨ur Quanteninformationsverarbeitung, Albert-Einstein-Allee 11, D-89069 Ulm, DeutschlandEuropean Centre for Theoretical Studies in Nuclear Physics and Related Areas, I-38050 Villazzano (TN), Italy
We analyze in detail the so-called “pushing gate” for trapped ions, introducing a time dependentharmonic approximation for the external motion. We show how to extract the average fidelity for thegate from the resulting semi-classical simulations. We characterize and quantify precisely all typesof errors coming from the quantum dynamics and reveal for the first time that slight nonlinearitiesin the ion-pushing force can have a dramatic effect on the adiabaticity of gate operation. By meansof quantum optimal control techniques, we show how to suppress each of the resulting gate errorsin order to reach a high fidelity compatible with scalable fault-tolerant quantum computing.
I. INTRODUCTION
Trapped ultracold ions have represented a major can-didate for the implementation of scalable quantum in-formation processing since the beginning of this researchfield. The first proposal of an ion-based quantum com-puter by Cirac and Zoller in 1995 [1] has been followedby a great variety of other schemes based on ions [2], onother quantum optical systems like neutral atoms [3, 4]and on solid-state systems [5] as well. With the progressof experimental techniques and the demonstration of en-tangling quantum gates based on several different candi-date physical systems, the focus has progressively shiftedtoward the fulfilment of scalability desiderata [6], that is,the realization of quantum gates with very high fidelities,in the range 0.999 – 0.9999.Gate errors in a real implementation of a given quan-tum gate scheme can be reduced by different means.Some errors arise from (or are increased by) experimentalimprecisions of a technical nature and can be controlledby careful alignment, stabilization etc. of the experimen-tal apparatus. Other errors stem from unaviodable in-teractions with the environment and can be reduced bysimply completing the gate in as short a time as possi-ble. Typically, a gate scheme can be made faster by sim-ple scaling to higher intensities, shorter distances etc. Ifsuch simple optimizations of the gate prove insufficient,one needs to consider changes to the scheme itself andtrade simplicity for improved performance. This is ex-actly the goal of quantum optimal control techniques [7]which allow for a precise tailoring of the system’s evo-lution by time-dependent tuning of some external pa-rameters. With sufficient control over these parameters,a given target state can often be reached with minimalerrors even over short gate operation times. The appli-cation of these methods to quantum information systemsrequires in turn a very accurate simulation of the dy- namics and a careful understanding of the targeted errorsources. This is precisely the aim of this paper, in thespecific case of the two-qubit ion gate first proposed in[2] and subsequently analyzed in [8, 9]. In this “push-ing gate” the qubits are encoded in the internal states oftwo ions. Each ion is held in a separate microtrap andstate-selective “push” potentials are applied in order tomodify the distance and thus the Coulomb interactionbetween ions, see Sec. III below. We shall first pointout a series of issues that arise when the assumption ofspatial homogeneity of the ion-driving force is dropped,and subsequently develop a way to correct each of theseissues, exploting a range of ideas, including in a crucialway optimal control methods.It should be noted from the outset that in the presentpaper, we analyze the pushing gate without what is calledthe π -pulse method in Refs. [8, 9], where it was shownto dramatically reduce some types of errors. The π -pulsemethod is a spin-echo technique and requires the gate tobe repeated with the internal state of the ions flipped.Typically, single particle operations like flipping the in-ternal states can be done with high fidelity, and it isreasonable to expect that eventually the π -pulse methodwill be used. However, internal state control is at least inprinciple a separate issue from the pushing gate opera-tion itself. Keeping the design process modular in spirit,it is relevant to optimize the gate without this additionaltrick and thus pave an alternative road to high fidelities.As will become clear below, we extract quite general noisereduction methods from the automated numerical opti-mization results and it is an interesting topic for futureresearch to combine these with the π -pulse method.The remainder of the paper is organized as follows. InSec. II we introduce the general setting of conditionaldynamics gate schemes, a useful approximation for sim-ulating such a gate, and a measure of the gate errors.In Sec. III we specialize to the pushing gate. The un-optimized performance of the gate is reported in Sec. IVand in Sec. V we show how this performance can be signif-icantly improved by a combination of “manual” changesand numerical optimal control methods. Finally we con-clude in Sec. VI. The Appendix contains a number ofmore technical results and derivations. II. CONDITIONAL DYNAMICS
The basic idea of the so-called pushing gate is one of conditional dynamics , i.e. we apply potentials depend-ing on the internal state of the two ions. The inter-nal states themselves are not changed during the gate,or, in the case of the π -pulse method, are changed ona much shorter time-scale than the external dynamics.This means that the analysis of the problem splits intofour separate evolutions for the external state, one foreach of the logical (internal) states, 00, 01, 10, and 11.In the following, these four evolutions will be denoted“branches” and will be indexed by β ∈ { , , , } .The complete Hamiltonian can be written in the form ofa sum of internal ⊗ external factorized termsˆ H tot ( t ) = X β | β ih β | ⊗ ˆ H ex β ( t ) (1)and it results in an evolution operator of a similar formˆ U tot ( t ) = X β | β ih β | ⊗ ˆ U ex β ( t ) . (2)Ideally, when t = T at the end of the gate, the U ex β should differ from each other by at most a phase factormultiplying a common unitary operator ˆ U excom ˆ U ex β ( T ) = e iθ β ˆ U excom (3)so that ˆ U tot itself can be factorized:ˆ U tot ( T ) = X β | β ih β | e iφ β ⊗ ˆ U excom . (4)The internal evolution is then that of a phase gate, whilethe external evolution can in principle be undone usinginternal-state independent potentials.The requirement (3) is very hard to achieve and wouldmake the gate completely independent on the initial ex-ternal state. However, we can typically assume to havesome degree of control over the initial external state, e.g.by cooling the particles before the gate. This meansthat (3) need only hold when restricted to a subset of thecomplete Hilbert space, typically the states of relativelylow energy. In Appendix B we show how to evaluate theperformance of the gate in general. For now, we notethat since only the low energy part of ˆ U ex β ( T ) will be im-portant, we can focus on getting a good approximationfor this part when trying to simulate the gate dynamics. A low initial energy means particles localized near thepotential minimum and this suggests using a harmonicapproximation to the real potential [22]. The simplestchoice is to Taylor expand around a fixed point, which isnot changed during the gate operation. The next level ofrefinement is to expand around the instantaneous poten-tial minimum. This works very well if the gate operationis nearly adiabatic so the particles stay near the (moving)minimum at all times. However, it may be desirable tomake fast and substantial changes to the potential duringthe gate and that may induce pronounced non-adiabaticdynamics. In that case, the harmonic approximation canstill be a good one provided it is done around the classicaltrajectories of the particles. Typically these trajectoriescannot be computed analytically, but for any moderatenumber of particles it is a numerically simple task to findthem. Below we will use this method and show how itleads to a relatively simple characterization of U ex β . A. Harmonic approximation
In this section we focus on a single branch of the evolu-tion and thus suppress the β index. Let us denote by x ( t )the classical trajectory, which is found by solving classi-cal equations of motion. The time-dependent, secondorder Taylor expansion of the potential V ( t, x ) around x ( t ) reads simply V so ( t, x ) = V (0) ( t )+∆ x T V (1) ( t )+ 12 ∆ x T V (2) ( t )∆ x , (5)with ∆ x = [ x − x ( t )] and V (0) ( t ) = V ( t, x ) ,V (1) i ( t ) = ∂V∂x i ( t, x ) ,V (2) ij ( t ) = ∂ V∂x i ∂x j ( t, x ) . (6)Note that we will still use x as our coordinate, i.e. weare not changing to a coordinate system moving with x ( t ), we simply use a potential that approximates thereal potential close to x ( t ). Collecting terms of equalorder in x leads to the alternative form V so ( t, x ) = E ( t ) − x T F ( t ) + 12 x T K ( t ) x , (7)with E ( t ) = V (0) ( t ) − x T ( t ) V (1) ( t ) + 12 x T ( t ) V (2) ( t ) x ( t ) , F ( t ) = − V (1) ( t ) + V (2) ( t ) x ( t ) ,K ( t ) = V (2) ( t ) . (8) B. Gaussian evolution
The big advantage of choosing a second order approx-imation to the real potential is that this restricts thecorresponding approximate ˆ U ex ( t ) to be Gaussian for all t . Let us introduce the compact notation q = ( x , p ) anddefine the matrix J by J = (cid:20) I n − I n (cid:21) , (9)where n is the number of degrees of freedom. Then theusual canonical commutation relations can be written as[ q i , q j ] = iJ ij . (10)With the potential of Eq. (7) and a matrix of parti-cle masses M = diag( m , . . . , m n ) the time-dependentHamiltonian becomesˆ H so = 12 ˆ p T M − ˆ p + V so ( t, ˆ x ) . (11)In Appendix A we show that such a Hamiltonian leadsto an evolution operator of the formˆ U ( t ) = e − iφ ( t ) ˆ D c ( t ) ˆ W b ( t ) , (12)where ˆ D c is a displacement operator and ˆ W b a squeezingoperator ˆ D c = e − i c T J ˆ q ˆ W b = e − i ˆ q T b ˆ q . (13)The scalar φ , the vector c , and the matrix S = exp( Jb )should satisfy the following equations of motion: ∂∂t φ = E − F T x ,∂∂t c = Jh c + (cid:20) F (cid:21) ,∂∂t S = JhS, (14)where the 2 n × n matrix h is defined by h ( t ) = (cid:20) K ( t ) 00 M − (cid:21) . (15)The form of solution (12)–(14) holds for any secondorder Hamiltonian. In the particular case where V so is aTaylor expansion of a real potential around the classicaltrajectory x ( t ), the equation of motion for c reduces tothe exact equation of motion for q = ( x , p ) where p is theclassical momentum. In the following, we will thereforewrite q instead of c . It is then important to rememberthat the right-hand sides of Eqs. (14) are in general non-linear functions of x . C. Fidelity
We can quantify the performance of the gate by calcu-lating the average fidelity F avg (see Appendix B) betweenthe obtained output state and the ideal one when the in-put state is varied. One can then separate out three kindsof contributions to the deviation of F avg from 1 (the per-fect gate) 1 − F avg = E θ + E q + E S (16)The three types of errors each have their physical inter-pretation. The most straightforward one pertains to the sloshing errors E q which correspond to a residual mo-tion of the ions after the gate has been completed andthe micro-traps are again at rest. The phase errors E θ are errors in the gate phase. Finally, the breathing er-rors E S are induced by differences in the harmonic ap-proximation parameters around the classical trajectoryfor different internal states. For example, in the case weconsider, when the particles are pushed closer together,the second order term in the Coulomb repulsion becomeslarger, cf. Eq. (25) below.In our model we assume that systematic, local phaseerrors can be undone. Then an explicit calculation inAppendix B, shows that the phase errors are given byE θ = 120 ( θ − θ − θ + θ − π ) (17)with θ β = − φ β +Tr[ b β γ ], where γ is the covariance matrixof the external state, see Section B 3 of the Appendix.Note the inclusion of the Tr[ b β γ ] terms in the definition of θ β : These terms correspond to phase contributions fromaverage excitation energy in the traps and are thereforetemperature dependent through γ . For our parametersthey are small.The sloshing errors are given byE q = 120 X α<β ( q α − q β ) T γ ( q α − q β ) , (18)and the breathing errors areE S = 140 X α<β Tr [( b α − b β ) γ ( b α − b β ) γ ]+ 1160 X α<β Tr [( b α − b β ) J ( b α − b β ) J ] . (19)Here φ β , q β , and b β are defined in Eqs. (12,13), and arethe variables describing the Gaussian approximation inthe β branch of the evolution. Again γ introduce tem-perature dependence. III. THE PUSHING GATE
Let us now focus on the particular case of the pushinggate. Here we have two ions, each in a separate micro-trap [23]. To further simplify the discussion, we concen-trate on just one spatial dimension, i.e. there are twodegrees of freedom, n = 2. The ions are assumed to be ofidentical mass, m = m = m , and thus M = m I . Thepotential energy consists of a micro-trap for each ion,time- and internal state-dependent pushing potentials,and the Coulomb interaction. The pushing potentialscan be realised as optical dipole potentials generated byfocused laser beams. The time-dependence of these po-tentials is most easily achieved by controlling the inten-sity of the laser and the state-selectivity by polarizationselection rules [9]. We assume that the form for the in-ternal state labeled by β is V ( β ) ( t, x ) = X i =1 , mω x i + e πǫ | d + x − x | + X i =1 , f ( β ) i ( t ) ~ ωa (cid:20) − ( − i x i + Ga x i (cid:21) . (20)The trapping potentials are assumed to be perfectly har-monic. The state-dependent pushing amplitudes f ( β ) i aresuch that the ions are only pushed if they are in the in-ternal state 1, f ( β ) i ( t ) = δ β i , f ( t ). Note that the ( − i factor means that ion 1 is pushed to the left and ion 2to the right. Non-linear contributions to the pushing po-tentials are included via the constant G . The harmonicoscillator ground state size is a = p ~ /mω . The twocoordinates x and x are taken to have origin in therespective trap centers, a distance d apart. Experimen-tally, the parameters in (20) can be varied quite a lot(see e.g. [9]). Trap distances d from 100 µm all the waydown to 1 µm are within technological reach. The trapfrequency ω can be chosen in the range 2 π × − Hzwhich for e.g. Ca + ions will mean an oscillator length a from 150nm down to 5nm.It should be noted that the use of optical dipole po-tentials to generate the state selective pushing forces willin general introduce large single qubit phases due to ac-Stark shifts. This is not a problem as such, but it meansthat even small fluctuations in laser intensities will leadto loss of gate fidelity. For the particular case of thepushing gate, the ac-Stark shifts can be balanced againstthe Coulomb energy as discussed in Ref. [9]. Obviouslyac-Stark shifts are common to many gate proposals thatuses optical potentials. An experimentally demanding,but quite general solution is to compensate the shiftsalong the lines of Refs. [10, 11]. In the present work, wefocus on errors that are more directly related to the mo-tion of the ions and assume that the push potentials areeffectively non-fluctuating.. A. Dimensionless Hamiltonian
The relative strength of the Coulomb interaction to thetrapping potentials turns out to be conveniently quanti- fied by ǫ = e πǫ d mω d = a d e πǫ d ~ ω , (21)which is the ratio of the energy scale of the Coulomb andtrap potential energies at the equilibrium positions ofthe ions. In Ref. [9] it was found that ǫ ≪ β readsˆ H ( β )push ( t ) = X i =1 , (cid:2) ˆ p i + ˆ x i (cid:3) + ǫ d a (cid:12)(cid:12) a d (ˆ x − ˆ x ) (cid:12)(cid:12) − X i =1 , f ( β ) i ( t ) (cid:2) ˆ x i + G ˆ x i (cid:3) . (22) B. Harmonic approximation
When Eqs. (8) are specialized to the pushing gate, weget the following: E ( β ) ( t ) = 14 ǫ a d x ( β )2 − x ( β )1 a d (cid:16) x ( β )2 − x ( β )1 (cid:17) (23) F β ) i ( t ) = − ( − i f ( β ) i ( t ) ∓ ǫ da a d (cid:16) x ( β )2 − x ( β )1 (cid:17)h a d (cid:16) x ( β )2 − x ( β )1 (cid:17)i (24) K ( β ) ii ( t ) = 1 + 2 f ( β ) i ( t ) G + 12 ǫ h a d (cid:16) x ( β )2 − x ( β )1 (cid:17)i (25) K ( β )12 ( t ) = K ( β )21 ( t ) = − ǫ h a d (cid:16) x ( β )2 − x ( β )1 (cid:17)i . (26)When these expressions are inserted into Eqs. (14) and(15) we are ready to simulate the gate. IV. RESULTS OF SIMULATIONA. Choice of parameters
Even with the simplifications we have introduced,there are still a lot of parameters in the problem. Theoptimal “working point” will always be dependent on ex-perimental considerations beyond the simplified modeltreated here. For a discussion of parameters and designdecisions, see Ref. [9]. For concreteness we have cho-sen to focus on a limited set of parameters. We first ofall assume the individual ion traps to be very well sep-arated and let a/d = 0 .
001 in all calculations. Likewise,we assume a reasonably low value for ǫ of 0 .
04. Suchparameters would result from e.g. Ca ions placed inmicro-traps with trapping frequencies of ω ∼ π × ∼ µ m. For a the trav-eling wave configuration with beam waist w consideredin Ref. [9], G = 4( a/w )( w/ x − x /w ) where x is theinitial position of the ion relative to the beam center.For realistic focusing of the push beam, w ∼ µ m thissuggests to vary the non-linearity coefficient G between0 and 3 × − . As the initial temporal shape of the pushpulse we choose a Gaussian f ( t ) = ξ exp( − t /τ ), wherethe amplitude ξ should be chosen to give a gate phaseof π . A simple estimate (for G = 0) suggests that wechoose [8] ξ = π p π/ ǫ p ǫ/ ωτ . (27)The temporal width of the pulse, τ , should be within anorder of magnitude from the trap period if we want a fastgate. We will mainly look at τ in the range 1 to 10 trapperiods, which for the parameters quoted above resultsa maximum excursion due to the push in the range from12 a down to 3 a . B. Phase errors
The choice of push amplitude expressed by Eq. (27) isnot optimal. This can be seen in Fig. 1 where we plotE θ as a function of G . Even for G = 0 the gate phaseis not exactly π . For ωτ = 3 . G can improve the gate phase. This is not surprizing,but also not very useful as we shall see below that E θ is in general easy to reduce. In Fig. 1 results for threedifferent temperatures are plotted, but the dependenceon temperature is completely negligible. C. Sloshing errors
Let us now turn to the errors described by the E q termin Eq. (C8), the sloshing errors. Fig. 2 shows how theseerrors are strongly dependent on G , the strength of thenon-linearity of the pushing potential. A series of minimaof E q as a function G can be seen. The optimal values of G depends on the chosen duration of the pulse, τ . Eachminima is associated with the ions performing an integernumber of “non-adiabatic oscillations” during the pushpulse. This is illustrated in Fig. 3 where the trajectoryof ion 1 with respect to its trap minimum is plotted forvalues of G that are below, at and above the one thatleads to the lowest E q .In contrast to E θ above, E q depends noticeably onwhether T is 0 . × ~ ω/k B . Higher temperaturesalways increase the sloshing errors and for k B T / ~ ω ≫ q scales as T since γ does [see Eqs. (B21)and (B22)].Curves for three different values of τ are plotted inFig. 2 and it is immediately clear that one can dramat- −8 −6 −4 −2 PSfrag replacements G E θ ωτ = 3 . ωτ = 5 . ωτ = 7 . FIG. 1: Phase errors with non-uniform pushing forces. Thephysical parameters are: ǫ = 0 . a/d = 0 . ωτ = 3 . . . T = 0 . , × ~ ω/k B are plotted for each value of τ and are indistinguishable from each other. −10 −8 −6 −4 −2 PSfrag replacements G E q ωτ = 7 . ωτ = 5 . ωτ = 3 . FIG. 2: Sloshing errors with non-uniform pushing forces. Thephysical parameters are: ǫ = 0 . a/d = 0 . ωτ = 3 . . . q as a functionof G . As explained in the text, there exist nonzero valuesof G where the sloshing is strongly suppressed. Results fortemperatures of T = 0 . × ~ ω/k B are plotted asdashed, full and dotted lines, respectively. −0.4 −0.2 0 0.2−0.3−0.2−0.100.10.20.30.4 PSfrag replacements x cl , − x min , p c l , − p m i n , PSfrag replacements x cl , − x min , p cl , − p min , FIG. 3: Phase-space trajectories for ion 1 for a push durationof ωτ = 5 .
5. The coordinates are relative to the potential min-imum in which ion 1 is trapped. Perfect adiabatic evolutionwould correspond to the ion simply following this minimumand thus the trajectory would be the single point (0 ,
0) in thisplot. Since the push is not infinitely slow, the ion will first lagbehind the moving minimum and then oscillate in the movingpotential. When the potential minimum again approaches itsoriginal position, the ion may happen to have just the rightposition and speed in order to end up at rest. Whether thisis the case depends (for fixed push pulse) on G : the higher G , the more the confinement is increased during the push. Inthe figure, the thick line corresponds to G = 0 . G = 0 , . , . , . G givesthe outermost curve over the main part of the “loop” in thefigure. ically decrease sloshing errors by make the gate slowerand thus more adiabatic. The suppression of E q is expo-nential and this is thus in general an efficient strategy. D. Breathing errors
We now turn to the E S term of Eq. (16). These“breathing” errors come from the different changes in theeffective quadratic Hamiltonian for the different branchesof the evolution. In Fig. 4 we plot E S for different valuesof G and τ and for different temperatures of the exter-nal motion. Results for temperatures of T =0.125, 1 and8 × ~ ω/k B are shown and it is first of all clear that E S depends strongly on T . We also see that E S is nearlyproportional to G τ . That larger G leads to larger errorsis not surprising, but that larger τ does is rather counter-intuitive: larger τ means a more smooth and thus moreadiabatic push. It also means a smaller amplitude for thepush since the ions will have more time to pick up the −8 −6 −4 −2 PSfrag replacements G E S ωτ = 3 . ωτ = 5 . ωτ = 7 . FIG. 4: Breathing errors with non-uniform pushing forces.The physical parameters are: ǫ = 0 . a/d = 0 . ωτ =3 . . . T = 0 . × ~ ω/k B are plotted as dashed, full anddotted lines, respectively. For all three temperatures, larger τ leads to larger E S . E S is increasing approximately linearlywith G τ , i.e. for fixed G , E S will be proportional to τ , thepulse duration! This is very different from the behavior of E q above, which is rapidly decreased by increasing τ and therebymaking the push more adiabatic. gate phase, cf. Eq. (27). Let us discuss the explanationfor this behaviour in more detail.For exponential suppression of errors to be valid, theevolution should be well into the adiabatic regime. Atfirst sight, the relevant timescale is ω − , the oscilla-tion period of the micro traps. Since the traps areassumed to be far apart ( a ≪ d ) the parameter ǫ issmall and the normal modes of the system have peri-ods shifted little from this value: the CM mode is infact unaffected by the Coulomb interaction and has fre-quency ω while the relative motion mode oscillates at √ ǫ ω . With ωτ ≫ transfer excitations betweenthe modes, as the adiabatic timescale for this process is( √ ǫ ω − ω ) − ∼ ǫ − ω − /
2. Such a transfer will beinduced by mixing of the CM and relative motion duringthe gate operation. A linear push potential will not mixthe two, but a non-linear one will.Another effect to remember is that when the instanta-neous oscillation frequencies change during the push, theexternal motion will pick up different phases dependingon the number of excitation quanta. Like the transfer ofexcitations, this effect of course disappears if the systemis cooled to the ground state. A perturbative calculationto lowest order in a/d , G and ǫ gives the result E S = π G ξ ω τ ×× ( ( ~ ω/ k B T ) (cid:20) (cid:18) − ǫ ω τ (cid:19)(cid:21) + 2 (cid:20) ( ~ ω/ k B T ) + 2 (cid:21) exp (cid:0) − ω τ (cid:1)) . (28)The prefactor gives the scaling behaviour both in the“na¨ıve” nonadiabatic limit ωτ < < ωτ < ǫ − :E S ∝ ξ G ω τ ∝ G ωτǫ , (29)In fact, even in the adiabatic limit ǫωτ ≫ τ . This un-suppressed term is stemming from the above mentionedeffect of time-varying instantaneous mode frequencies.From Eq. (28) we can also understand the strongtemperature dependence of the breathing errors. For k B T / ~ ω ≫
1, the breathing errors will scale approxi-mately like T . However, as seen in Figure 4, high tem-peratures require very low values for G . For low tem-peratures, note that one term in Eq. (28) is not sup-pressed even at T = 0. This term stems from changes inthe ground state widths of the two instantaneous normalmodes and is adiabatically suppressed when ωτ ≫
1. Tofind the dominant term for very low temperatures andshort pulses one should do a higher order perturbativecalculation.
V. OPTIMIZING GATE PERFORMANCE
From the simulations in Sec. IV we learn that withoutimprovement high fidelities require either very low valuesof G or cooling the external motion almost to the groundstate. In this section we shall see how a better perfor-mance can be achieved by modifying the temporal shapeof the push-pulses. A. Correcting the phase
Our first step will be to correct the gate phase by asimple scaling of the push pulse-shape. From Fig. 1 weknow that typical errors can be well above the percentlevel. Our strategy is based on the observation that thesimplest estimate of the gate-phase suggest that it scalesas the square of the push-amplitude, ξ . [A general gate-phase replaces the π in the numerator of Eq. (27).] Wetherefore divide ξ by the square-root of the ratio of theobserved gate-phase and the ideal gate-phase ( π ) and re-peat the propagation. In Fig. 5 we show the results of ap- −8 −6 −4 −2 PSfrag replacements G E θ n o o p t i m i z a t i o n i t e r a t i o n i t e r a t i o n s FIG. 5: Improving E θ by iteratively adjusting the push-amplitude ξ . The pulse-duration is τ = 7 . T = 0 .
125 (dashed), 1 (full),and 8 × ~ ω/k B (dotted). The uppermost set (looks like a singlecurve) is for the un-optimized pulse-amplitude and identicalto the τ = 7 . θ below 10 − for all the considered G < . plying this algorithm to the τ = 7 . θ is rapidly reduced and can be brought be-low e.g. 10 − in a very modest number of iterations. Forsimplicity we ignore the temperature dependent Tr[ bγ ]contribution to E θ when rescaling the pulse. This is thereason for the k B T = 8 ~ ω (dotted curves) departing fromthe k B T = 0 . ~ ω and k B T = 1 ~ ω curves especially atlow G . B. Fast gate: eliminating sloshing in x
A big advantage of the simple harmonic approximationis that it becomes feasible to solve the equations of mo-tion many times with different temporal shapes of f ( t )in order to optimize the performance of the gate. Ratherthan simple trial-and-error we will apply the global con-trol algorithm of Krotov, which is guaranteed to improvethe performance at each iteration [12–14]. The relevantequations for our case are given in Appendix C.In general, it is desirable to complete the gate in asshort a time as possible. This will limit many undesiredeffects and it will ultimately enable faster quantum com-putations. A fast gate, however, means that the pushingforce will deliver a rather abrupt impulse. This can lead −10 −8 −6 −4 −2 PSfrag replacements iterations E q ωτ = 3 . ωτ = 5 . ωτ = 7 . FIG. 6: Optimization of pulse-shape to eliminate sloshing.We plot E q , the contribution of sloshing to the total infidelity,as a function of the number of Krotov iterations performed.The parameters in this example are: ǫ = 0 . a/d = 0 . G = 2 × − . Each curve corresponds to a separate valueof the pulse duration, ωτ =3.5, 5.5, 7.5, with larger τ alwaysgiving a lower value of E q . to excitations of the external motion being left after thecompletion of the gate, limiting the fidelity. In this sec-tion we show how such “sloshing” effects can be avoidedby using optimal control.We start from an initial Gaussian temporal shape ofthe push. The overall amplitude is first optimized itera-tively to get the desired gate phase as described above.We then run the Krotov algorithm to get a better shapeof the pulse. We assume a non-uniform pushing force, G = 2 × − . The result is plotted in Fig. 6. As can beseen, the influence of sloshing motion can be decreasedby a couple of orders of magnitude in a modest numberof iterations.To investigate the physical mechanism behind the re-duction of the sloshing error, we plot in the lower panel ofFig. 7 the difference between the optimized pulse f opt ( t )and the original Gaussian pulse f ( t ) = ξ exp( − t /τ ).This difference looks a lot like a simple cosine-wavewith a period close to 2 πω − multiplied by a Gaus-sian of the same width as f . Thus the optimizedpulse is approximately of the form f opt ( t ) ∼ f ( t ) + A cos( ωt ) exp( − t /τ ). An approximative calculation ofthe sloshing excitation for the simplest G = 0 case re-veals that for such a pulse, the non-resonant contribu-tion of the bare Gaussian pulse is cancelled by a resonant contribution from the cosine-modulated pulse. Since theresonant response is much stronger, only a small, nega-tive A is needed for this cancellation. More precisely, theoptimal A from first order pertubation theory is given −20 0 20−1−0.500.51 −20 0 20−10−50510−20 −10 0 10 20051015 PSfrag replacements ωt f o p t − f ξ e x p ( − ω τ / ) G = 0 G = 2 × − ωt ωtωt f FIG. 7: Upper panel: pulse shape both before (full line) andafter (dashed line) Krotov optimization. Curves for three dif-ferent pulse duration are shown, ωτ = 3 .
5, 5.5 and 7.5, otherparameters are as in Fig. 6. Only for the shortest pulse is theoptimized curve distinguishable from the original Gaussian.Lower panel: Difference between optimized pulse and initialGaussian pulse (after adjusting overall amplitude to reducephase errors). The left plot is for G = 0, the right plot isfor G = 2 × − . Each curve has been normalized to theprediction of a G = 0 perturbative calculation, which is seento describe well the G = 0 case as all curves have a maximalexcursion of approximately -1, while the G = 0 case is onlyqualitatively similar. by − ξ exp( − ω τ /
4) and Fig. 7 shows that this is alsowhat the Krotov algorithm converges to for G = 0. The“strategy” of the Krotov algorithm in this case seemstherefore to be well understood. For G = 0, it is moredifficult to predict the value of A , but nonetheless theKrotov algorithm seems to be highly efficient. C. Minimizing breathing errors
We now know that phase errors and sloshing errorscan be controlled and we turn to the breathing errorsof Fig. 4. Without optimization, these errors put ratherstringent limits on the parameters. In order to keep E S atan acceptable level, either very low temperature or verysmall G is required. For very low temperature k B T =0 . ~ ω , we need just G . − to get E S below 10 − ,but if we assume a more modest cooling to k B T = 1 ~ ω the same error-level require G . × − . Note that evenat G = 0, breathing errors persist and that for k B T =8 ~ ω they never get below the 10 − level. These errorsstem from the high order terms in the Coulomb potentialwhich have been ignored in Eq. (28).As described in Sec. IV D, breathing errors cannot beeliminated by simply increasing the push duration τ :First of all the adiabatic time-scale is ∼ ǫ − ω − / π -pulse method,see Ref. [9]. Allowing negative push amplitudes andputting “by hand” an optimized cos( ǫωt/
2) modulatedcontribution, we have been able to e.g. reduce E S below10 − for ωτ = 7 . G = 2 × − , and T = 1 × ~ ω/k B .Compared to the results reported in Fig. 4, this is a re-duction by more than 3 orders of magnitude. Unfortu-nately, the strongly modified pulse now gives rise to largesloshing errors. To obtain an overall safisfactory fidelitywe use a combined strategy: we first put the breathingerror reduction by hand, then iteratively reduce phaseerrors and finally use the Krotov algorithm to reduce thesloshing errors. In Fig. 8 we show results of this strategystarting from ωτ = 7 . VI. CONCLUSION
In this paper we have shown how a time-dependent,quadratic approximation to the Hamiltonian can be auseful tool when analyzing quantum gates based on con-ditional external dynamics. The resulting equationsare much more manageable than the original two-bodySchr¨odinger equations. This is especially true if one wereto include more spatial dimensions than the one consid-ered here: A full time-dependent, three dimensional, two-body wavefunction calculation is an extremely demand-ing numerical task, whereas the corresponding quadraticapproximation will be much more manageable.We used the developed method to show how to im-prove on a na¨ıve design of the pushing gate. This was −8 −6 −4 −2 PSfrag replacementsiterations E θ Krotov
FIG. 8: Combined strategy for reducing infidelity for ǫ = 0 . G = 2 × − . The three contributions E θ , E q , and E S are labeled by ▽ , ◦ , and (cid:3) , respectively. Their sum is la-beled by × . The leftmost column are the results for a simpleGaussian with ωτ = 7 .
5. The breathing errors dominate. Inthe next column, breathing errors have been reduced (butsloshing increased) by using a cosine-modulated pulse basedon Eq. (28). The overall amplitude is then iteratively opti-mized to reduce phase errors as described in Sec. V A. Ascan be seen, 3 iterations are more than sufficient to render E θ completely insignificant. The dominating error type is nowsloshing and in the third column, the Krotov algorithm isapplied to finally reduced the total infidelity below 10 − . done including a non-uniform contribution to the pushingforce. An important lesson of our analysis and simula-tions, is to pay attention to changes in the symmetryof the Hamiltonian during the gate operation. In thepresent case, a nonlinear push potential invalidates theseparation of the dynamics into CM and relative motion.This opens up another type of non-adiabaticity, namelytransfer of excitations between the two normal modes.The adiabaticity parameter for this type of error is ǫωτ and for small ǫ , adiabaticity will require gate times muchlarger than the charateristic time of the micro-traps. Anefficient counter-measure is to decrease temperature sothat there are in fact no excitations to transfer betweenmodes. Failing that, one should increase ǫ as much aspossible and, somewhat counter-intuitively, do the gateas fast possible. Sloshing errors puts a lower limit onthe gate-time, but as we show an optimized choice of thetemporal shape of the push can dramatically reduce thisproblem.By analyzing the way that the Krotov optimized pulsereduces sloshing errors, we identified the basic mechanismas a destructive interference between the non-resonant,non-adiabatic contribution from the finite push pulse du-ration and a resonant contribution from a small ampli-tude superposed oscillation of the push force. Generaliz-0ing this idea to deal with breathing errors, we were ableto reduce them by several orders of magnitude. However,since breathing errors are not significantly adiabaticallysuppressed for the considered pulse durations, the de-structive interference required sign-changes in the pushamplitude, which introduce experimental complications.The suppression of breathing errors also came at the priceof increased sloshing errors, but we showed that the Kro-tov algorithm once again was able to improve the pulseshape.One may ask to what extent is the final pulse shapeoptimal for the given overall gate-time? This is an in-teresting question in general and in this study we sawexamples both where the Krotov algorithm seemed toexhaust the potential in its “strategy” (the G = 0 casein Fig. 7) and where it was not able to find an optimiza-tion. In the latter case we could improve the pulse byhand (eliminate the breathing errors by destructive in-terference). The problem of optimality is related to thequestion of a quantum speed limit (QSL) [15] and thisconnection has been studied in Ref. [16]. Note, how-ever, that in our case the hamiltonian is time-dependentand what we want is in fact to leave the external motionunaffected after the pulse. It would be interesting to in-vestigate such a general adiabaticity problem along thesame lines as the work on the QSL. Acknowledgments
This research was supported in part by the NationalScience Foundation under Grant No. PHY05-51164, andin part by the EC projects AQUTE and EMALI.
Appendix A: Evolution under quadratic Hamiltonian
In this section we show that a second order Hamilto-nian leads to an evolution operator that can be writtenlike ˆ U ( t ) Eq. (12). For alternative parameterizations, seeRefs.[17, 18]. We will do a direct calculation showingthat ˆ U ( t ) fulfills the Schr¨odinger equation i ∂∂t ˆ U = ˆ H so ˆ U . (A1)The demanding part of the the calculation involves dif-ferentiating expontials of time dependent operators. Auseful formular can be found in e.g. [19] and involves in-tegration over an auxillary variable η . It results in i ∂∂t ˆ D c = Z ˆ D η c ˙ c T J ˆ q ˆ D † η c dη ˆ D c = Z ˙ c T J (ˆ q − η c ) dη ˆ D c = ˙ c T J (cid:18) ˆ q − c (cid:19) ˆ D c (A2) for the displacement operator and i ∂∂t ˆ W b = Z ˆ W ηb
12 ˆ q T ˙ b ˆ q ˆ W † ηb dη ˆ W b = 12 Z (cid:0) e − ηJb ˆ q (cid:1) T ˙ b (cid:0) e − ηJb ˆ q (cid:1) dη ˆ W b = 12 ˆ q T J T Z e ηJb J ˙ be − ηJb dη ˆ q ˆ W b = 12 ˆ q T J T (cid:18) ∂∂t e Jb (cid:19) e − Jb ˆ q ˆ W b (A3)for the squeezing operator. It is now easy to show thatthe equations of motion (14) for φ , c and S = exp( Jb )leads to ˆ U = exp( − iφ ) ˆ D c ˆ W b fulfilling Eq. (A1) [24]. Appendix B: Fidelity for Gaussian evolutions
In this section we derive expressions for the fidelityas a function of the variables used to characterize theevolution in the harmonic approximation, φ β , q β , and S β , β = 00 , , , ρ ⊗ σ , we get the following statefor the internal degrees of freedom after the applicationof ˆ U tot of Eq. (2): ρ ′ =Tr ex (cid:20) ˆ U tot ρ ⊗ σ (cid:16) ˆ U tot (cid:17) † (cid:21) = X αβ | α ih β | [ ρ ◦ R ] αβ . (B1)Here “ ◦ ” denotes the element-wise matrix product (theHadamard product) and R is the matrix given by:[ R ] αβ = Tr (cid:20) ˆ U ex α σ (cid:16) ˆ U ex β (cid:17) † (cid:21) (B2)It is easy to see that R er Hermitian and that all itsdiagonal elements are 1. In particular, Tr R = 4. Slightlyless obvious is it that R is positive semi-definite: Let c ∈ C . Then: c † Rc = Tr (cid:20)(cid:8)X c ∗ α ˆ U ex α (cid:9) σ (cid:8)X c β (cid:16) ˆ U ex β (cid:17) † (cid:9)(cid:21) = Tr h(cid:8)X c ∗ β ˆ U ex β (cid:9) † (cid:8)X c ∗ α ˆ U ex α (cid:9) σ i ≥ R = P h w h w † h , we get instead a Krauss operator sum form: ρ ′ = X h =1 K h ρK † h (B4)1with K h = diag( w h ).There are different ways to define the fidelity of thegate. The one used in Refs. [9] and [8] is as the minimumfidelity of the obtained final state w.r.t. the wanted finalstate when the input state is varied. This means that F min = min ψ h ψ | U † ρ ′ ψ U | ψ i (B5)with ρ ψ = | ψ ih ψ | and U the gate operator we aim for.In the case of a phase gate, U is diagonal in the logicalstate basis: ˆ U = X β | β ih β | e iθ β (B6)and we get the simpler minimization problem: F min = min { p i } p T ˜ Rp (B7)where the p ∈ R and P p i = 1. The matrix ˜ R = U RU † is Hermitian, but since p is confined to be real, only itsreal symmetric part contributes. The minimization inEq. (B7) is a so-called quadratic programming problemand very efficient numerical methods for its solution ex-ist. Given ˜ R it is therefore simple to calculate F min onthe computer. However, a more direct evaluation is pos-sible if fidelity is instead defined as an average over inputstates as F avg = Z S n − h ψ | U † ρ ′ ψ U | ψ i dV. (B8)Here S n − denotes the normalized states (unit sphere) in C n and the volume element dV is such that R S n − dV =1. For F avg a compact formula exist [20] and using it inthe present case leads to F avg = 14(4 + 1) Tr ˜ R + X αβ ˜ R αβ (B9)or when using the properties of ˜ R :1 − F avg = 110 X α<β (cid:16) − Re h ˜ R αβ i(cid:17) . (B10)
1. General small errors
Typically, we will be mostly interested in situationswhere the four logical states leads to almost identicalevolutions for the external states. It is then useful towrite ˆ U ex β = exp( i ˆ D β ) ˆ U excom = exp (cid:16) i h ˆ D β i (cid:17) exp( i ∆ ˆ D β ) ˆ U excom , (B11) where h ˆ D i = Tr h ˆ D ˆ U excom σ ex ( ˆ U excom ) † i and ∆ ˆ D = ˆ D − h ˆ D i .Calculating ˜ R to second order in the ∆ ˆ D ’s, we get:˜ R αβ = exp ( i ∆ θ α ) × (cid:26) − D ∆ ˆ D α + ∆ ˆ D β −
2∆ ˆ D α ∆ ˆ D β E(cid:27) × exp ( − i ∆ θ β ) . (B12)This form is useful, as it separates the infidelity into sys-tematic phase-errors (the exp( ± i ∆ θ α ) factors) and “de-coherence” (factor in curly brackets). The phase errors∆ θ β = h ˆ D β i − θ β can be made small by tuning the av-erage of laser powers etc. and this can usually be donevery well. The challenge will therefore most often be tosuppress the fluctuations, i.e., the terms in curly bracketsin Eq. (B12).Assuming also the ∆ θ ’s to be small, Eq. (B10) be-comes:1 − F avg = 120 X α<β (cid:18) [∆ θ α − ∆ θ β ] + (cid:28)(cid:16) ∆ ˆ D α − ∆ ˆ D β (cid:17) (cid:29)(cid:19) . (B13)At first sight, this form might seem dubious since only differences in the ∆ θ ’s and ∆ ˆ D ’s enter. However, oneshould remember that any common evolution on the fourbranches can be absorbed into ˆ U excom in Eq. (B11). Thisemphasizes that for the implementation of a single gateon the logical state, the external motion must not neces-sarily be returned to its initial state as long as the finalstate is common to all logical input states. Typically,the further requirement that energy is not pumped intothe external degrees of freedom by repeated applicationof the gate must be made. In the particular case of thepushing gate, this requirement is in fact already hidden inEq. (B13) since the β = 00 branch contains no pushing.In other cases, one could apply cooling to the externalstate between gate operations.
2. The non-local part of the phase
We are seeking to implement the phase-gate (B6). Inmany cases, the θ β ’s are not so important individuallysince single-particle operations are easy to perform andonly the truely non-local phase is interesting. Assumingthat perfect single-particle phase changes can be imple-mented on average , it is straightforward to show that oneshould replace P α<β [∆ θ α − ∆ θ β ] by[∆ θ − ∆ θ − ∆ θ + ∆ θ ] (B14)in Eq. (B13).This simplified view of single-particle phase-changesshould of course be revisited in a more complete anal-ysis of any given proposal for quantum-computing. In2the present work we use the replacement (B14) through-out, but let us emphasize that fluctuations in the single-particle phaserotations are more naturally incorporatedin the ∆ ˆ D terms of Eq. (B13) than in the ∆ θ terms: Onesimply model the fluctations as a consequence of somefluctuating parameter which can be included in σ ex .
3. The Gaussian case
For Gaussian evolutions like (12) and a Gaussian (e.g.thermal) external state σ with covariance matrix γ ij = 12 h q i q j + q j q i i − h q i ih q j i =ReTr [ q i q j σ ] − Tr [ q i σ ] Tr [ q j σ ] (B15)and vanishing means h q i i = Tr [ q i σ ] = 0 (B16)one gets phase contributions h ˆ D β i = − φ β − Tr [ b β γ ] (B17)and decoherence terms D (∆ D α − ∆ D β ) E = ( c α − c β ) T JγJ T ( c α − c β )+ 12 Tr [( b α − b β ) γ ( b α − b β ) γ ]+ 18 Tr [( b α − b β ) J ( b α − b β ) J ] (B18)In general, for a harmonic oscillator in thermal equi-librium at temperature T , the covariance matrix is givenby γ thermal = 12 1tanh ~ ω k B T (cid:20) ~ mω ~ mω (cid:21) (B19)where k B is the Boltzmann constant. In the present case,the CM and the relative motion are separately in ther-mal equilibrium and for the corresponding dimensionlessposition and momentum operators [ x CM = ( x + x ) / γ = γ CM ⊕ γ rel (B20)with γ rel = 12 1tanh (1+ ǫ ) / ~ ω k B T " ǫ ) / (1+ ǫ ) / (B21)and γ CM = 12 1tanh ~ ω k B T (cid:20)
00 2 (cid:21) . (B22)In the limit ǫ ≪
1, a we have approximately γ ∝ I if weuse the set of individual-ion operators ( x , x , p , p ). Appendix C: The Krotov Algorithm
Optimizing the temporal shape of the push pulse isdone using the Krotov algorithm [12]. For an introduc-tion to the method, see e.g. Ref. [13].
1. Auxillary variables
The key ingredient in this approach is a functionΦ( t, { φ β , q β , S β } β =00 , , , ) which allow us to translatethe global goal of improving the final ˆ U tot ( T ) to a local problem of choosing a better f ( t ) for each t . Construct-ing Φ is in general very difficult, but it is relatively simpleto get a linear approximation to it. The coefficients inthis approximation will constitute a set of auxillary vari-ables. For each branch, the equations of motion for theauxillary variables ˜ φ β , ˜q β and ˜ S β are determined by therequirement that they are conjugate to the physical vari-ables φ β , q β , and S β , respectively.Let us focus on a single branch and suppres the β index like in Sec. II A. We then need to construct H ( t ; φ, q , S ; ˜ φ, ˜q , ˜ S ) such that Eqs. (14) can be written ∂∂t φ = ∂∂ ˜ φ H (C1) ∂∂t q i = ∂∂ ˜ q i H (C2) ∂∂t S ij = ∂∂ ˜ S ij H . (C3)This leads simply to H ( t ; φ, q , S ; ˜ φ, ˜q , ˜ S ) = ˜ φ h E ( t, x ) − F T ( t, x ) i + ˜q T h Jh ( t, x ) q − J F ( t, x ) i + Tr h ˜ S T Jh ( t, x ) S i (C4)Then the equations of motion for the auxillary variablesbecome ∂∂t ˜ φ = − ∂∂φ H = 0 (C5) ∂∂t ˜q = −∇ q H (C6) ∂∂t ˜ S = ∇ S H = h ( t, x ) J ˜ S. (C7)The equation of motion for ˜q is rather involved since H depends on x in a complicated manner through h , F and E . It can be rewritten as two coupled time-dependent,forced harmonic oscillators. Note on the other hand that˜ φ is time-independent and that J ˜ S solves the same equa-tion as S .3
2. Objective function
Our ultimate goal is to improve the fidelity of the gate.However, it is somewhat impractical to apply this as theobjective in the Krotov algorithm: Calculating the fi-delity is only simple for small errors and in general itdepends on e.g. the temperature of the external motion.Instead we shall work with a simpler function of the vari-ables φ , q and S for the four branches. The reduction ofthis objective function should tend to increase the fidelityof the gate. Based on the fact that in the pushing gate,the branch β = 00 is not subject to any time-dependentforces, we choose the following: J (cid:0) { φ β , q β , S β } β =00 , , , (cid:1) = J φ + J q + J S = 12 h φ − φ − φ + φ − π i + 12 X β (cid:20)(cid:16) x β − x (cid:17) + (cid:16) p β − p (cid:17) (cid:21) + 12 X β Tr (cid:20)(cid:16) S β − S (cid:17) T (cid:16) S β − S (cid:17)(cid:21) . (C8)The term with φ ’s aim to ensure the correct phase in thephasegate, while the other terms aim at identical evolu-tion for the external motin in the four branches. In thelimit ǫ ≪
3. Terminal conditions
The objective function supplements the auxillary-variable equations of motion (C5–C7) with the following terminal conditions , i.e. boundary conditions at t = T :˜ φ β ( T ) = − ∂∂φ β J (cid:12)(cid:12)(cid:12)(cid:12) T = − ( − β (cid:0) φ − φ − φ + φ − π (cid:1)(cid:12)(cid:12) T (C9) ˜x β ( T ) = −∇ x β J (cid:12)(cid:12) T = − (cid:16) x − x − x − x (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 − (cid:16) x β − x (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 (C10) ˜p β ( T ) = −∇ p β J (cid:12)(cid:12)(cid:12) T = − (cid:16) p − p − p − p (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 − (cid:16) p β − p (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 (C11)˜ S β ( T ) = −∇ S β J (cid:12)(cid:12) T = − (cid:16) S − S − S − S (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 − (cid:16) S β − S (cid:17)(cid:12)(cid:12)(cid:12) T , β = 00 (C12)where ( − β is +1 for β = 00 and 11 and − β =01 and 10. These equations express the values of theauxillary varibles at time t = T in terms of the physicalvaribles also at t = T and give the input to the backwardspropagation of the auxillary variables, cf. Ref. [12]. [1] J. I. Cirac and P. Zoller, Phys. Rev. Lett. , 4091 (1995).[2] J. I. Cirac and P. Zoller, Nature , 579 (2000).[3] D. Jaksch, H.-J. Briegel, J. I. Cirac, C. W. Gardiner, andP. Zoller, Phys. Rev. Lett. , 1975 (1999).[4] T. Calarco, E. A. Hinds, D. Jaksch, J. Schmiedmayer,J. I. Cirac, and P. Zoller, Phys. Rev. A , 022304 (2000).[5] G. Burkard, D. Loss, and D. P. DiVincenzo, Phys. Rev.B , 2070 (1999).[6] D. P. DiVincenzo, Fortschr. Phys. , 771 (2000).[7] D. D’Alessandro, Introduction to quantum control anddynamics , Chapman & Hall/CRC applied mathematicsand nonlinear science series (CRC Press, Boca Raton,FL, 2008).[8] T. Calarco, J. I. Cirac, and P. Zoller, Phys. Rev. A ,062304 (2001).[9] M. ˇSaˇsura and A. M. Steane, Phys. Rev. A , 062318(2003).[10] H. H¨affner, S. Gulde, M. Riebe, G. Lancaster, C. Becher,J. Eschner, F. Schmidt-Kaler, and R. Blatt, Phys. Rev.Lett. , 143602 (2003). [11] A. Kaplan, M.F. Andersen, and N. Davidson, Phys. Rev.A , 045401 (2002).[12] V. F. Krotov, Global methods in optimal control theory ,no. 195 in Monographs in pure and applied mathematics(Marcel Dekker, Inc., New York, 1995).[13] S. E. Sklarz and D. J. Tannor, Phys. Rev. A , 053619(2002).[14] D. Tannor, V. Kazakov, and V. Orlov, in Time Depen-dent Quantum Molecular Dynamics , edited by J. Broeck-hove and L. Lathouwers (Plenum Press, New York, 1992),vol. 299 of