Laser control in a bifurcating region
1 Laser control in a bifurcating region
D. Sugny*, C. Kontz*, M. Ndong + , Y. Justum + , G. Dive ** and M. Desouter-Lecomte ‡1 *Laboratoire de Physique de l’Université de Bourgogne, Unité Mixte de Recherches 5027 CNRS et Université de Bourgogne, BP 47870, 21078 Dijon, France. + Laboratoire de Chimie Physique, Unité Mixte de Recherches 8000, CNRS et Université de Paris-Sud-11, 91405 Orsay Cedex, France. **Centre d’Ingéniérie des Protéines, Université de Liège, Sart Tilman B6, B-4000 Liège, Belgium. ‡ Laboratoire de Chimie Physique, Unité Mixte de Recherches 8000, CNRS et Université de Paris-Sud-11, 91405 Orsay Cedex, France and Département de Chimie, Université de Liège, Institut de Chimie B6, Sart-Tilman, B-4000, Liège 1, Belgium. Corresponding author : [email protected] CO which is a molecule-benchmark for such dynamics. We show the feasibility of different processes and we investigate their robustness against variations of laser field. We discuss the conditions under which each method of control gives the best results. We also point out the relation between optimal control theory and local control. 3
I. Introduction.
Control of physico-chemical processes by ultra-short laser pulses remains nowadays an attractive and challenging domain. The aim of this kind of control is to design a laser pulse which drives the system from an initial state to a specific target state or even better, to find laser fields able to perform unitary transformations on molecular qubits. By this way, shaped laser pulses have become new reagents for chemical reactions. Some of the most important experimental contributions to this field have been reviewed recently [1]. On the other hand, different control schemes have been proposed, among others we can cite the Brumer-Shapiro coherent control [2 3], the Tannor-Rice-Kosloff local control approach [4 5], the Rabitz optimum control theory (OCT) based on learning algorithms or closed-loop control procedures [6 7 8 9] or the simulated Raman adiabatic passage (STIRAP) scheme [10 11 12]. This article is devoted to a theoretical analysis of different scenarios based on STIRAP (or extension of this process as f-STIRAP [13]) and OCT by working only in the Infra-Red domain, i.e., without transitions via excited electronic states. We consider a two-dimensional model of a bifurcating region in the ground potential energy surface of a polyatomic system. Such a region connects three non equivalent wells. A deep reactant well is connected to a symmetric double-well. One passes from the reactant well to the double basin with a large amplitude bending mode of a migrating hydrogen atom around a given bond. The double-well corresponds to an internal rotation of this atom around the axis defined by the particular bond. This three-well bifurcating region is an interesting pinball topography which suggests different processes of control : - transformation of a delocalized state into a localized state in the double-well potential [14] - transformation of a localized state of the double-well potential into the other. This has been already proposed in the spirit of Cope rearrangement [15] or enantiomer selection [16] - isomerization from the reactant well to a given basin of the surface like in hydrogen transfer in organic molecules [17]. In our case, this reaction involves a break of symmetry. 4- realization of one or two qubits systems [18 19 20 21 22 23]. The double-well region offers different possibilities for the choice of the quantum numbers which allow to define the qubits (parity or excitation) [24] We address different control issues: the efficiency of various strategies which depend on the shape of the dipolar surface or, equivalently on the structure of the dipolar matrix and the robustness of the control with respect to the process used and the duration of the pulse. In each case, we also analyse the different pathways which are enforced by the laser field. Finally we briefly discuss the relationship between local control and OCT in the particular case where the objective is to maximize the average value of the projector on a superposed state.
II. Model
We consider a model recently proposed which reproduces the main features of a bifurcating region connecting three potential energy wells [25]. Isoenergy contours are presented in Fig. 1. The model is calibrated on an ab initio computation at the QCISD level [25] of the isomerization of the methoxy radical into hydroxymethyl which is a molecule benchmark for such energy landscapes [26]. We should emphasize that we do not intend to control dynamics of this particular radical moiety but, we are rather interested in this particular topography for which a convenient analytical expression has been proposed. Figure 1 near here 5FIG. 1. Isoenergy contours (in eV) in the model potential energy surface of the isomerization H CO → H COH as a function of two active angular coordinates (see Fig.2). The zero of energy is at the bottom of the product well ( P or P ’); R =0.181 eV, TS1 = 1.854eV and
TS2 = 0.195eV. The model describes the rotation of the hydrogen atom around a polar bond connecting two atoms (here CO) in different chemical environments (see Fig.2). The two active coordinates [ ] θ π∈ and [ ] , φ π π∈ − are the spherical angles of the migrating hydrogen atom with respect to the center of the bond. In C S geometry ( φ = 0), the first active bending coordinate θ connects the reactant well R to a second well trough a first transition state TS1 (the barrier height from the reactant is 1.673 eV). This second well is a transition state
TS2 according to a second symmetry breaking active coordinate φ . TS2 is the top of the small barrier (0.195eV) of the double well corresponding to rotational conformers P and P ’. Between TS1 and
TS2 lies a valley ridge inflexion point (VRI). Mathematical definitions of a VRI point can be found in different works [27 28 29 30]. Roughly speaking, it is a point where a valley corresponding to a particular internal mode becomes an unstable ridge. We point out that field-free dynamics has already been carried out in such bifurcating regions by assuming that initial wave packets can be prepared in the valley uphill from the VRI [31 32]. In the quasi-harmonic regime, one can introduce vibrational quantum numbers for θ and φ oscillators. The ground state of the R well is denoted by
0, 0 R . The delocalized states of the double well are of parity even and odd and are thus noted , n m + and , n m − . The splitting of the first level + , − is 4.3 10 -5 eV. This corresponds to a rather long tunneling time of about 95 ps much longer than the duration of the pulses used in the control. The first localized states coming from the in phase and out of phase superposition are ( ) L,0 , 0 , 0 / 2 n n n = + + − and ( )
R,0 , 0 , 0 / 2 n n n = + − − . They are associated to m = 0 for the θ vibrator. We recall that the notations R and L do not refer to enantiomers in this example (R is P and L is P ’). In the dipolar approximation, the reduced 2D Hamiltonian takes the form 6 ˆ ˆ ˆ ( ) k kk H H E t µ= − ∑ (1) where ˆ ˆ ˆ H T V = + is the field free Hamiltonian and k denotes the polarization direction. The exact constrained 2D kinetic energy operator can be numerically computed by the TNUM algorithm [33] by freezing the inactive coordinates at the TS1 geometry. We extract an approximate kinetic energy operator by fitting the standard angular momentum expression in spherical coordinates. In Euclidian normalization convention, ˆ T is then equal to
1ˆ ( cotan )2 2 sin
Eucl
T I I θ φ θθ θ θ φ∂ ∂ ∂= − + −∂ ∂ ∂ (cid:61) (cid:61) where constant inertia moments I θ = 6160 a.u. and I φ = 4430 a.u. are estimated from the TNUM grids. In the Wilson normalization convention in which the volume element is d d θ φ , ˆ T becomes: ˆ ( )2 2 sin Wil
T vI I θ φ θθ θ φ ∂ ∂= − − +∂ ∂ (cid:61) (cid:61) where )( θ v is an extra-potential term. This analytical expression is particularly suited to the use of the split operator algorithm [34] which is needed to propagate the wave packets. This point is due to the fact that the coefficient of a given differential operator / k q ∂ ∂ does not depend on k q but only on the other coordinates. We assume that the molecules are aligned in the laboratory frame with the polar bond oriented along the z e (cid:71) axis (see Fig.2). This could be obviously an important constraint [35]. We consider linear polarizations with directions x e (cid:71) in the C S plane and y e (cid:71) perpendicular to the C S plane. Figure 2 near here 7 FIG. 2. Active coordinates θ and φ for the isomerization H CO → H COH and polarization directions for aligned molecules. To pursue the construction of the model, we also propose a simple form for dipolar surfaces based on a chemical analysis of the molecule. Indeed, the system can be roughly described as the rotation of a charged particle around a polar bond. The dipolar components ( , ) x µ θφ and ( , ) y µ θφ are larger on the P , P ’ side ( θ < π /2) when the particle is close to the most electronegative atom. They decrease quickly for θ > π /2, when the particle enters a region near the weakly electronegative atom. The analytical model is given in the appendix. µ x ( θ , φ ) and µ y ( θ , φ ) are respectively symmetric and antisymmetric with respect to φ . Some cuts are given in Fig.3. Figure 3 near here C OH θ HH φ zxy C OH θ HH φ zxy
012 0 60 120 180 isomerization coordinate θ (degree) d i po l a r m o m en t u m ( D eb y e ) φ = 0° φ = 75° φ = 25° φ = 50° φ = 50° φ = 25° . Cuts in the model dipolar momentum surface (see Appendix) for different values of the torsion angle φ . Full lines : ( , ) x µ θφ ( symmetric), open circles: ( , ) y µ θφ ( antisymmetric). III Control methodologies A. Local and optimal methods
The control algorithms are usually classified as local [4 36 37 38] or global depending on whether the field is determined from the instantaneous dynamical properties by maximizing a performance index or from the variational calculus of a cost functional. The objective functional can be defined in different manners [7 8] which are strongly connected [9]. The procedure to maximize the cost functional under the constraint of satisfying the time dependent Schrödinger equation is described in details in the literature [39]. The Zhu, Botina, Rabitz formulation [7] leads to three coupled equations: the Schrödinger equation for ( ) t ψ with an initial condition ( 0) i i t ψ φ = = (forward propagation), the Schrödinger equation for the Lagrange multiplier ( ) f t ψ with a final target condition ( ) f f T ψ φ = (backward propagation) and an equation for the optimum field ( ) ( ) 1/ ( ) ( ) ( ) ( ) j i f f j i E t m t t t t α ψ ψ ψ µ ψ = − ℑ (cid:61) (2) where α is a positive penalty factor chosen to weight the significance of the laser fluence. An experimental switching function ( ) ( ) sin / s t t T π = is usually introduced [39], α is then replaced by / ( ) s t α α → . The equations are solved by an iterative formulation [7] adapted to a discrete implementation based on a second order Split Operator scheme [34]. We have used the improvement proposed in Ref. [40]. At each iteration, the field is given by ( ) ( 1) ( ) k k kj j j E E E − = + ∆ where ( ) kj E ∆ is calculated by [Eq.(2)]. It is worthy to note that the local approach is strongly related to the Zhu, Botina and Rabitz approach when the performance index involves a projection on a non stationary state. The local control methodology 9is overviewed in Ref. [38]. The field is chosen in order to maximize the rate of variation of a performance index ( ) ˆ( ) ( ) j y t y O t = which is a function of expectation values ˆ ˆ( ) ( ) ( ) ( ) j i j i O t t O t t ψ ψ = of Hermitian operators with j = 1, N . In the case where the target operator is a projector on a non stationary wave packet at a final time T : ˆ ( ) f f O T φ φ = , the rate depends on a single expectation value ˆ( ) / ( ) / dy t dt d O t dt = . If the time dependence of the operator is fixed by the field free Hamiltonian ( )/ ( )/ ˆ ( ) ( ) ( ) iH t T iH t Tf f f f O t e e t t φ φ φ φ − − − = = (cid:61) (cid:61) , (3) in other words, when the operator projects down to the wave packet which freely evolves towards the target state at time T , then one obtains [38] ( ) ˆˆ( ) / 2 ( ) . ( ) dy t dt m O t E t µ = − ℑ (cid:71)(cid:71) . (4) The local control field giving a monotonous increase of the performance index is obtained by setting, for a polarization direction, ( ) ˆ ˆ( ) ( ) j j j E t m O t λ µ = − ℑ . By inserting expression (3) into this last equation, one gets an expression corresponding to the first step (without zero order field) of the iterative optimum control [Eq.(2)], i.e. when ( ) f t ψ evolves with the field free Hamiltonian ( ) ˆ( ) ( ) ( ) ( ) ( ) j j i f f j i E t m t t t t λ ψ φ φ µ ψ = − ℑ (5) The method focuses on the λ j coefficient. With few trials, it is possible to find values of λ j providing an acceptable field. The latter is then used as an initial-guess field to continue the iterative optimum control procedure. This speeds up the rate of convergence of the algorithm by finally choosing the best α . The optimum field able to steer a set of initial states to a set of target states, i.e. to apply a unitary transformation to the 2 N states of N qubits 10 ˆ: : N N f igatef i U φ φφ φ = can be obtained by the multitarget generalization of OCT [18 21]. We have to propagate simultaneously a set of 2 N wave packets forward in time ( 0) n ni i t ψ φ = = with 1,...2 N n = and a set of 2 N Lagrange multipliers wave packets backwards ( 0) n nf f t ψ φ = = with 1,...2 N n = . The optimum field is given by a sum of contributions from each state ( )
20 1 ( ) 1/ ( ) ( ) ( ) ( ) N n n n nj i f f j in E t m t t t t α ψ ψ ψ µ ψ = = − ℑ ∑ (cid:61) (6) A constraint on the phase of the quantum gate could be added [41]. The fidelity of the quantum gate is measured by ˆ ˆ( ) / 2 Ngate control
F tr U U + = (7) B. STIRAP and adiabatic processes
The second strategy for the control is based on adiabatic passage (for a recent overview, see [11 12] and references therein). Such processes are widely used in a variety of fields, extending from nuclear magnetic resonance and quantum information to atomic and molecular excitations. Adiabatic methods are usually achieved by using a series of intense pulses which can be frequency-chirped, the frequencies and the chirping being adapted to the structure of the energy levels. However, the modification of the shape of the pulse envelope and the chirping rate must be sufficiently slow so as to fulfill adiabatic conditions. One of the most well-known adiabatic processes is the STIRAP excitation which involves a counterintuitive sequence of two pulses in a three-level system, in which the field of the Stokes pulse precedes and overlaps the field of the pump pulse. These adiabatic techniques allow a complete population transfer from an initial state to a target state which can be either a stationary state i.e. an eigenstate of the field free Hamiltonian or a coherent superposition of such states. They are also robust in the sense that they are not sensitive to small 11variations of laser parameters. Due to these remarkable properties, such processes seem to be particularly suitable for the control of chemical reactions. For instance, they have been applied with success for controlling the isomerization of HCN [42 43]. However, the relevance of adiabatic techniques in a complex system can be questioned. We stress that 100% efficiency of the control is generally ensured only for a subset of levels with particular couplings such as the tripod system. If the molecular system is rich in the energy range considered, the effect of coupling to background states can deteriorate noticeably the population transfer in particular if the background states are resonant or almost resonant with laser fields. Another major drawback of these methods is the duration of the pulses which is longer than the time needed by optimal or local control to reach their objective. This point can be problematic if other concurrent chemical processes with time-scale of the same order occur during the control. To avoid the preceding problems, we combine in this article adiabatic processes which allow determining a simple form for the overall field and optimization of some parameters of the pulse, leading to a shorter (of the order of few picoseconds) and efficient control. The strategy can be summarized as follows. We first select a subset of levels and we determine an adiabatic process in order to achieve the objective of the control. These levels have to be carefully chosen, as otherwise the value of the electric field is too large. More precisely, we recall that the Rabi frequency )( tE µ =Ω between the states 1 and 2 ( µ being the matrix element of the dipole moment) must be sufficiently large so as to fulfilled adiabatic conditions. For instance, a standard condition is >>Ω T where T is a characteristic duration of the pulse, which is the full width half maximum for a Gaussian pulse. In addition, in order to avoid other unwanted chemical processes such as ionization, the intensity of the electric field has to be limited to /10 cmW which roughly leads to a minimum of the order of 0.1 a.u. for matrix elements of the dipole moment. In a second step, considering all the levels of the system, we decrease the pulse duration to few picoseconds and we optimize both the intensities and the delay between the different pulses to keep efficient control. We now describe the computational details of the method. We have used Gaussian pulses, the pulses being polarized in the x e (cid:71) or the y e (cid:71) direction. The field E ( t ) is equal to the sum of terms of the following form : 12)cos()2/)(exp( kkk tttE ϕωγ +−− (8) where γ , k ω , E and k ϕ are respectively the width, the frequency, the amplitude and the phase of the pulse. To simplify even more the overall field, we assume that the width and the amplitude are the same for all the pulses (except for the quantum gates). The delay is defined by the difference between the times k t . IV Wave packet c ontrol A. Double well scenarios
We consider two control schemes in the double well product region ( P and P’, see Fig.1): the localization of the ground delocalized state into one localized state, in the spirit of the previous control on H POSH [44] and the transformation from a localized state of one well ( P ) to a localized state of the other well ( P ’) [15, 16]. We schematize these processes as follows + → ( ) + + − = (9) → ( ) + − − = (10) (i) Description of the adiabatic processes We first analyze in details the processes L + → or R + → which consist in preparing one of the conformer from a delocalized state. The presentation of the results follows the different steps of the strategy. We begin by selecting the first three levels of the system, that is + , − and + . The method for determining the adiabatic process consists in using the particular symmetry of the dipole moment. For instance, we recall that x µ only couples the levels + and + , the transition − to + being forbidden. We consider a f-STIRAP scheme which, as the STIRAP technique, only uses two pulses, the pump and the Stokes fields. We choose to fix the frequencies k ω and the phases k ϕ of each pulse as follows: 13 kk E E E ωϕ + + − = − += (11) where E + is, for instance, the energy of the level + . In the three-state basis + , − and + , the total Hamiltonian can be written as cos( ) cos( )cos( ) cos( )cos( ) cos( ) S PS SP S
E t tt E tt t E α ω ωα ω ωω ω + − + Ω Ω Ω Ω Ω Ω (12) where P Ω and S Ω are respectively the Rabi frequencies of the pump and the Stokes pulses for the transitions + → + and − → + . The α parameter is equal to / x y µ µ + − + + . The pump pulse is polarized along the x e (cid:71) direction whereas the Stokes pulse is polarized along the y e (cid:71) one. Note that the Rabi frequencies are chosen real without loss of generality. Using the RWA approximation [12], the Hamiltonian I H reads in the interaction representation ΩΩ ΩΩ=
000 00
SP SPI H (13) where the small detunings E E ω + + − − and E E ω + − − − are neglected. The idea is then to use a f-STIRAP technique. f-STIRAP is an extension of STIRAP which allows the creation of coherent superpositions of states [11]. f-STIRAP is now a well-known process which has already been used in a variety of systems to implement qubit gates or to generate superposed states [45 46]. The process can be schematized by the following diagram 14 |0+,0 〉 |0-,0 〉 |1+,0 〉 |1-,0 〉 |2+,0 〉 |2-,0 〉 pump (E x ) Stokes (E y ) The extension of the STIRAP technique consists in the fact that the amplitudes of the two pulses are required to have a constant ratio at the end of the pulse. More precisely, if the eigenvector ψ of eigenvalue 0 of I H writes S PP S ψ = Ω + − Ω −Ω + Ω (14) then the following conditions have to be fulfilled by the two Rabi frequencies: =ΩΩ −∞→
SPt and ε=ΩΩ +∞→
SPt lim (15) where ±=ε . One deduces for the two limit cases that: ( ) 0 , 01( ) [ 0 , 0 0 , 0 ]2 ψψ ε−∞ = ++∞ = + − − (16) It is then clear that for an adiabatic evolution, the localized state can be obtained from f-STIRAP with +=ε for and −=ε for . Another equivalent scheme can be constructed by replacing the intermediate state + of the f-STIRAP process with + : 15 |0+,0 〉 |0-,0 〉 |1+,0 〉 |1-,0 〉 |2+,0 〉 |2-,0 〉 pump (E x ) Stokes (E y ) One can also imagine other mechanisms of the same kind using other intermediate states and the particular symmetry of the dipole moment. Finally, the following points can be noticed. A more complex superposed state can be obtained with f-STIRAP if the ratio SP ΩΩ / is different from 1 or -1 when +∞→ t . Moreover, a process using only one linear polarized laser field but with a frequency ω and its second harmonic ω can also be built to control the tunneling [47]. For the transformation → from a localized state to the other, we have slightly modified the previous scheme. The limits of [Eq. (15)] become ε−=ΩΩ −∞→ SPt lim and ε=ΩΩ +∞→
SPt lim (17) leading thus to the following limit states
1( ) [ 0 , 0 0 , 0 ]21( ) [ 0 , 0 0 , 0 ]2 ψ εψ ε−∞ = + + −+∞ = + − − (18) which correspond either to or according to the value of the parameter ε . Having determined an adiabatic process able to control the specified reaction, we are now in a position to examine the conditions (choice of Rabi frequencies and delay between the pulses) under which 16the mentioned scheme of control continues to work for a shorter duration of the pulse, i.e. not in the adiabatic limit. The optimized laser field has not been constructed using optimal algorithms in order to preserve as much as possible the robustness of the solution which, as stated above, is one of the most important features of adiabatic processes. For that purpose, we have considered a 2-dimensional grid (Rabi frequencies, delay) and we have calculated for each point of the grid, i.e. for particular values of Rabi frequencies and delay, the corresponding time evolution. The Rabi frequency (the same for each pulse) varies from − a.u. to − a.u. whereas the limits of the delay are CO, this amounts to a pulse duration of about 20 ps and a field amplitude of Vm − that corresponds to a very weak field. Figure 4 near here P opu l a t i on (a) | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | Lo c a li z a t i on (b) | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | -10 -5 0 5 1000.20.40.60.81 x 10 -4 R ab i f r equen c i e s time(ps) (c) FIG. 4. Dynamics controlled by f-STIRAP strategy for the preparation of the superposed state through the intermediate state + . Panels (a) and (b) show respectively the evolution of populations in the Hamiltonian eigenbasis and in the superposed states and . Populations of other vibrational states remain small during the process. The Rabi frequencies of the different pulses are displayed on panel (c). Rabi frequencies are in atomic units. The solid line corresponds to the Stokes pulse and the dashed line to the pump pulse. The total duration of the pulse is of the order of 20 ps. Figure 4 illustrates the results of applying the f-STIRAP strategy for a total duration of 20 ps and for the intermediate state + . In the adiabatic limit, only states + and − are expected to be populated. 18A different behavior is obtained for the process. This is due to coupling to background states and to the fact that adiabatic conditions are not rigorously fulfilled. For instance, the product Rabi frequencies*duration of the pulse is of the order of 10. However, this deviation from the theoretical description does not decrease its efficiency. The same behavior can be obtained for different intermediate states and different total durations. Figure 5 near here -5 Delay R ab i f r equen c i e s -5 Delay R ab i f r equen c i e s FIG. 5. Robustness of the f-STIRAP process as a function of the Rabi frequencies and the delay between the pulses for a total duration of 20 ps (upper part) and 4.5 ps (lower part) of the overall field. Rabi frequencies and delay are in atomic units. The intermediate state is + . 19The robustness of the strategies has been checked against two parameters: the time delay between successive pulses and the Rabi frequency of each pulse. Figure 5 shows the robustness against these two variables for the f-STIRAP process and for two total durations of about 20 ps and 4.5 ps (see Fig. 4). Remarkable robustness is achieved, especially for the longer pulse duration, advocating for a possible experimental feasibility of the control scheme. Moreover, it can be clearly seen that the strategy is more robust for longer pulses. This point can be explained by the fact that larger the duration of the pulse is, the more the adiabatic conditions are fulfilled and consequently the more the process is robust [48]. The same study with similar results can be done for the transformation from a localized state to the other. However, we notice that the overall dynamics is generally much more oscillatory with more complex structures as compared to the preceding reaction. (ii) Optimal control Methods based on local or global control allow finding optimum fields with a smaller duration T than fields obtained by the adiabatic approach. Fig. 9 displays the evolution of populations in the eigenbasis of ˆ H [Eq.(1)], in particular for the states + and − and the average value of the operator φ showing the localization in the P ’ well ( φ = - 75°) at the end of the process, after about 4.5 ps. This time is chosen because it is the shortest time ensuring a good performance index in the STIRAP approach. The objective is reached at 99.99% in 10 iterations with the Rabitz algorithm [Eq.(2)] improved by the correction proposed in Ref.[40] using α = 1.2 without any zero order field. Focusing on the first step of the procedure and using local control [Eq.(5)], we obtain a zero order field with λ x = 8 and λ y = 1.2 leading to a performance index of 91%. The Rabitz algorithm then converges at 99.99% in 3 iterations. Figure 6 near here 20 FIG. 6. Dynamics controlled by the optimum field for the preparation of . Left axis: evolution of populations in eigenvectors + and − , the excited states 3, 4, 5 nearly correspond to + , − and + respectively; right axis: evolution of the average φ position. One observes that the populations of the + and − eigenstates become equal very early but the average φ position shows that the equality of populations does not involve the correct phase to form the localized superposition ( ) + + − = . Transient excitations help at reaching the target superposition. The OCT field is shown in Fig. 7. A part of the structure of the pulse can be understood as follows. The pulse is composed of two sub-pulses, one along the x e (cid:71) direction and the other along the y e (cid:71) direction. We consider this latter part. This sub-pulse can be viewed as a half-cycle pulse (HCP) [49], i.e. only one half of an optical field cycle. HCPs have already been used in different applications; we can cite the control of molecular alignment or orientation [49 50] or the control of tunneling in a double-well system [51]. This pulse, being of short duration with respect to the tunneling time but long in comparison with the period associated to the transitions + → + or − → + , produces a superposition of states popu l a t i on -80-400 ♦ ( deg r ee ) ( ) ( ) t t ψ φψ t ψ− t ψ+ + and − . Using the sudden approximation [52], the evolution operator HCP U for the HCP can be written as follows in the basis defined by + and − : )exp( xHCPHCP iAU σ= (19) where HCP A is the area of the pulse times the corresponding matrix element of the dipole moment and σ x the Pauli matrix. Starting from a delocalized state, it can then be shown that a HCP of area π and a free evolution of a quarter of the tunneling time lead to a completely localized state [51]. A numerical calculation shows that the area of the optimal pulse along the y e (cid:71) direction is very close to π . Notice that the condition max Ω << ∆ where max Ω is the peak Rabi frequency and ∆ the detuning has to be fulfilled to avoid the appearance of other resonances and the transfer of population to excited states. ∆ can be roughly estimated by E E + + ∆ = − which leads to the following larger possible value of the electric field E − = × a.u. The second part of the pulse along the x e (cid:71) direction is a more complex field which cannot be explained so simply as it does not respect the condition max Ω << ∆ . This field is used by OCT to decrease the duration of the control from 20 ps to 4.5 ps. Figure 7 near here f i e l d E x ( V c m - ) -0.2-0.10 f i e l d E y ( V c m - ) E x E y + → [Eq. (9)]. We have also tested the robustness of this process against the area of the different pulses. As can be seen in [Eq. (19)], the area is the main feature of such short pulses. We have observed that the process is robust (of the order of 10%) with respect to such inaccuracies affecting the area of the pulse. It seems that this feature can be attributed to the simple form of the optimal field. B. Bifurcation scenario
We now investigate the possibility of steering the ground state
0, 0 R of the R well towards the P or P ’ product basin through the VRI region (see Fig.1). This control is summarized by the following schematic diagram:
0, 0 R → or
0, 0 R → . (20) This scenario involves a break of symmetry after a passage over a high barrier (1.67 eV). From a dynamical point of view, it can be wondered which path the controlled wave packet will follow, i.e. whether the bifurcation occurs early (near the VRI) or not (near TS2 ). It should be noted that this example illustrates the extreme sensibility of the control to parameters of the model. This scenario looks like isomerization processes between two wells which have already discussed in the literature [53 37]. However, the topography between
TS1 and
TS2 with a change of curvature along φ leads to delocalized eigenvectors strongly coupled by the dipolar momentum. This is unfavourable to the STIRAP scheme which needs intermediate states well decoupled from all the others. In OCT, we do not succeed in finding a satisfactory optimum field for the current model inspired from the QCISD ab initio level with a TS1 barrier of 1.67eV. The algorithm finds a path involving a too high excitation ˆ( ) ( ) t H t Ψ Ψ of several eV up to 7eV which is completely unrealistic. We adopt another potential energy surface inspired from other ab initio calculations 23(MP2) with a smaller energy barrier from the reactant (1.56 eV) at
TS1 ( R = 0.451eV, TS1 = 1.911eV and
TS2 = 0.216eV) and a slightly different profile along φ for φ > 80°. We use the same dipolar momentum model. In this case, the OCT gives a reasonable field leading to an average unperturbed energy of the order of E TS1 as shown in Fig.8. For a pulse duration of 4.5 ps, the target is reached with 95.3% in 240 iterations starting with two zero order fields ( ) cos( ) j j
E t E t ω= where E = a.u. and the ω j are the harmonic frequencies of the two θ and φ vibrators in the R well. Fig.8 shows the average value of the two active coordinates during the process. The controlled break of symmetry occurs in a sequential manner. The θ angle first reaches the TS2 value before the break of symmetry and the cooling occurs in the double well region. Figure 8 near here ene r g y ( e V ) -80-1060130 ang l e ( deg r ee ) ( ) ( ) t t θΨ Ψ ( ) ( ) t t φΨ Ψ ˆ( ) ( ) t H t Ψ Ψ
FIG. 8. OCT for the isomerization R → P ’. Black line: average unperturbed energy, gray lines: average value of the two active coordinates. Figures 9 and 10 give the optimal fields and the corresponding Gabor Transforms ( , ) ( , ) ( ) i s F t H s t E s e ds ω ω τ +∞−∞ = − ∫ (21) where H ( s, τ ) is the Blackman window [54] 24 H s s s sH s π π ττ τ ττ = + + ≤= and τ is the time-resolution. Here we have fixed τ = 0.2 ps. The Gabor transforms contain the zero order frequencies (1715 cm -1 for E x and 1578 cm -1 for E y ). A lot of frequencies are used during the time interval [0.8, 1.5] ps. They permit to increase the unperturbed energy above TS1 . The intermediate states playing a significant role in this heating are the low excitations of the θ vibrator (nearly R ,
0, 2 R and
0, 3 R ) and the first excitation in the φ vibrator
1, 0 R . Between 1.5 ps and 3 ps, a lot of delocalized states are populated with a weight smaller than 5%. At this point, the OCT path involves a large number of intermediate states. Note that this is the extreme opposite of the situation favorable for applying the STIRAP technique. The cooling occurs after 3 ps and mainly involves two states of the double well region (nearly + and − ). In this example, cooling is easier than heating probably because the dipolar momentum is very different in the two regions. Figure 9 near here FIG. 9. OCT fields for the izomerization R → P ’. Black line: E x ( t ), gray lines with open circles: E y ( t ) -1.2-0.8-0.400.40.8 0 1.5 3 4.5time (ps) f i e l d s ( V c m - ) E x E y R → P ’. We have also tried without any success the local approach [37] which consists in heating and cooling the wave packet according to the average θ position. Besides the difficulty of the expected break of symmetry, the average θ position does not reach the TS1 value during the heating but remains in the corresponding well excluding an efficient cooling.
V. Logical gates . We examine the possibility of realizing logical gates on one or two qubits. We recall that a quantum computation can be described as a sequence of logical gates which determine a unitary transformation gate U [55]. Different molecular systems such as vibrationally excited molecules [18 19 20 21 22 23 24] have already been proposed for the implementation of one and two-qubits gates and several control schemes using either pi-pulses [22] or optimal control theory [18 19 20 21 22] have been constructed. In the present case, the low lying states can be thought of as a qubit = + and = − . The previous transformation + → is obviously related to the well known Hadamard transformation 26 + + = = − −− . Following this idea, it can be shown that arbitrary unitary operations can be performed on the preceding qubit. For that purpose, we can use a universal set of one-qubit gates composed of the rotation gate and the phase gate which is defined by : i e ϕ + + = − − The basic transformation on a two-qubit system is the controlled-not gate which permutes the state of the second qubit only if the first qubit is in state 1:
00 001 0 0 001 010 1 0 0CNOT 10 100 0 0 111 110 0 1 0 =
There are different ways of defining a two-qubit system in our example. According to the proposal of Sola et al [22], we can choose excitation-parity or parity-excitation of the φ vibrator. This gives, respectively, the following definitions 00 0 , 001 0 , 010 1 , 011 1 , 0 + − = + − or 00 0 , 001 1 , 010 0 , 011 1 , 0 + + = − − We explore also the usual realization of a two-qubit system using states of two φ and θ vibrators in the double well 00 0 , 001 0 ,110 0 , 011 0 ,1 + + = − − . In the next section, we will give examples of control which aim at implementing the Hadamard gate, the phase gate and the C-NOT gate. 27 A. One-qubit gate (i) Adiabatic process
More complex methods than f-STIRAP strategy have to be used for implementing qubit gates. For one-qubit gates, we follow schemes proposed in [56 57]. We only consider the phase gate. The Hadamard gate can be derived by using a similar strategy. The procedure is composed of two STIRAP processes which are aimed at transferring the population between states − and + via state + which is not populated in the adiabatic limit. Only three pulses can be used because the second one serves as a pump field for the first STIRAP excitation and as a Stokes field for the second STIRAP process. The first part of the scheme can be viewed as follows : |0+,0 〉 |0-,0 〉 |1+,0 〉 |1-,0 〉 |2+,0 〉 |2-,0 〉 pump (E x ) Stokes (E y ) The frequencies are chosen so that the different excitations are resonant, i.e. we have xy E EE E ωω + ++ − = −= − where x ω and y ω are respectively the frequencies of the fields along the x e (cid:71) and y e (cid:71) directions. The phase of the first pulse is fixed to ϕ , the phase of the phase gate, whereas other phases are chosen to be zero. Note that other procedures with different phases can also be considered [56 57]. The eigenvector ψ of eigenvalue 0 of I H in the basis − , + and + reads as follows during the first STIRAP excitation : 28 iS PP S e ϕ ψ = Ω − − Ω +Ω + Ω and the second STIRAP ( 1 , 0 0 , 0 ) i S PP S e ϕ ψ −= Ω + − Ω −Ω + Ω In these expressions, the Rabi frequencies are assumed to be real and the dependence on the phase ϕ has been explicitly written in order to clarify the proof. Figure 11 near here P opu l a t i on (a) | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | -8 -4 0 4 800.20.40.60.81 F i de li t y time(ps) (b) -8 -4 0 4 800.20.40.60.81 F i de li t y time(ps) (c) -8 -4 0 4 802468 x 10 -5 R ab i f r equen c i e s time(ps) (d) FIG. 11. The phase gate. Panel (a) displays the evolution of populations in the Hamiltonian eigenbasis during the phase gate transformation, the initial state is − . Panels (b) and (c) represent respectively the evolution of the fidelity for πϕ= and 2 πϕ= . Panel (d) shows the Rabi frequencies of the different pulses. The solid and dashed lines correspond respectively to the field x E and the field y E (see text). It has been more difficult to control and to optimize one qubit gates than processes involved in the double-well scenarios. This point is basically due to the fact that both populations and relative phases have to be controlled in a quantum gate. This difficulty is particularly relevant in this case because the levels of the qubit are not degenerate. The correct unitary transformation is therefore achieved by the adiabatic process only in the interaction representation and not in the bare state basis. The optimization allows to set up the relative phases. Very good results have nevertheless been obtained. Moreover, one of the advantages of 30adiabatic processes is that the same form of the overall field can be used to realize different quantum gates. This point is illustrated for the phase gate in Fig. 11. Modifying only the phase of the first pulse of the first STIRAP excitation and keeping constant other parameters, two phase gates for 4 πϕ= and 2 πϕ= have been built by our strategy. (ii) Optimal control OCT confirms its efficiency in order to find fields of smaller duration. We have derived a field for the Hadamard gate on the + and − states with T = 4.5 ps. The field is completely similar to the one used to realize the scheme (9). One obtains exactly the same behavior as shown in Fig.6. This field is very simple and has been discussed in Sec. IV.A. We have also checked that OCT can be used to implement the phase gate. B. Two-qubit gate (i) Adiabatic process
We only consider the first C-NOT gate involving the states + , − , + and − . Similar processes can be constructed for other choices of qubits. The control scheme that can be used for the C-NOT gate is a strategy similar in its spirit to the preceding process. We make a step further by considering now a superposition of states. The scheme can be represented as follows : 31 |0+,0 〉 |0-,0 〉 |1+,0 〉 |1-,0 〉 |2+,0 〉 |2-,0 〉 pump (E x ,E y ) Stokes (E y ) where three different pulses have been considered. The values of the frequencies for the pump pulses and the Stokes field are chosen resonant with the corresponding transition. We consider the subset of levels + , − , + and − . In this basis, the total Hamiltonian I H can be written as ΩΩΩ ΩΩΩ −+ −+
11 11
S S where the Rabi frequencies (with straightforward notations) are assumed to be real. Diagonalizing the matrix C-NOT, we determine the corresponding eigenvectors involving the states + and − . These eigenvectors denoted + h and − h of eigenvalues 1 and -1 can be defined as follows : −−+= −++= −+ )0,10,1(21 )0,10,1(21 hh In the basis + h , − h , + and − , I H is given by ΩΩΩ ΩΩΩ −+ −+ S S Ω−Ω=Ω Ω+Ω=Ω −+− −++ )(21 )(21
11 11
The idea is then to decouple the eigenvector + h from other states of the basis. For instance, if we choose + Ω = Ω and − Ω = − Ω , one obtains for I H : ΩΩ ΩΩ
00 000 000 0000
S S
The last step consists in applying the scheme of the phase gate described above for a phase equal to π . In the adiabatic limit, + h and − h will be respectively transformed into + h and − − h which corresponds to the transformation of the C-NOT gate. Figure 12 near here P opu l a t i on (a) | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | | 〈 ψ (t) 〉 | -6 -2 2 6 1000.20.40.60.81 F i de li t y time(ps) (b) -6 -2 2 6 10-5-4-3-2-1012345 x 10 -5 R ab i f r equen c i e s time(ps) (c) FIG. 12. Same as Fig. 11 but for the C-NOT gate involving the states + , − , + and − . Figure 12 shows the results of this strategy. We have obtained a fidelity close to 0.95. The fact that the levels of the two qubits are not degenerate implies a quick loss of the fidelity of the order of 1 ps which seems problematic in view of experimental applications. We emphasize that this behavior can be observed in most of quantum gates constructed from vibrationally excited states. It is a disadvantage of this kind of systems in comparison of other schemes such as optical cavity [56] where all states are degenerate or almost degenerate. (ii) Optimal control We present only the gate C-NOT on the two qubits using fundamental and first excited states of two φ and θ vibrators in the double well. We impose the pulse duration T = 4.5 ps. Fig. 13 displays the time evolution of 34the population when each initial state
00 0 , 0 = + ,
01 0 ,1 = + ,
10 0 , 0 = − and
11 0 ,1 = − is driven by the optimum field which has been obtained with 21 iterations. Panels (a) and (b) show the inversion of population of the states of the second qubit. Gray lines display intermediate populations of the different eigenstates. Panels (c) and (d) show the population of the first qubit states. The final value is again equal to one at the end of the process even if intermediate depopulation occurs. Figure 13 near here
Time (ps) P opu l a t i on
10 0 ,0 = −
11 0 ,1 = −
Time (ps) P opu l a t i on
00 0 ,0 = +
01 0 ,1 = +
Time (ps) P opu l a t i on
46 8
11 0 ,1 = −
10 0 ,0 = −
Time (ps) P opu l a t i on
01 0 ,1 = +
00 0 ,0 = + (a) (d)(c) (b)
FIG. 13. C-NOT gate on the + − + − states. Black lines: population of the states of the two-qubit system, gray lines: intermediate transitions towards other eigenstates which are denoted from 3 to 10 and nearly correspond to excitation of the even and odd states of the φ vibrator only. The θ vibrator remains in its ground state with no node along θ . The optimal field obtained for this C-NOT gate is given in Figure 14. Only the E x component is used by the OCT. The maximum of the weak E y component is of the order of 1.5 10 -3 Vcm -1 . The field is again very 35simple. The Gabor transform [Eq.(21)] shows that a main frequency 1916 cm -1 corresponding to the − → − transition acts during the whole process. Figure 14 near here FIG. 14. OCT field for the C-NOT gate on the + − + − states. VI Concluding remarks
This article has focused on the application of OCT and adiabatic processes to various situations that can be encountered when a potential energy surface presents a bifurcating region connecting three potential wells (isomerization, tunneling and implementation of one or two qubits quantum gates). In the present case, the symmetric double well region is the most favorable to realize control scenarios due to the shape of the dipolar momentum surface. The results are expected to be transposable to other molecules such as H POSH presenting the same kind of double well region. We have also investigated the advantages and limits of the different methods. We recall that the goal of a control is to reach a defined objective, the field solution being subject to some physical constraints: on -3-113 0 1.5 3 4.5time (ps) f i e l d ( V c m - ) E x Appendix
This appendix gives the analytical expression of the dipole moment that has been used in our calculations. We first define the function CS µ as an approximation of the x µ in C S plane ( φ = ): ( ) cos ( ) kcs kk a µ θ θ = = ∑ where the parameters are given by a a a a a = = = = − = − . The two active coordinates of the dipole moment are then equal to : ( ) [ ]
21 23 ( , ) ( ) ( ) cos( ) ( ) 0.25 cos ( ) 1.75 cos( ) 1( , ) ( ) ( ) sin( ) x csy cs f ff µ θφ µ θ θ φ θ φ φµ θφ µ θ θ φ = + + − = where
1( ) 3 0.92 21( ) 3 0.92 21( ) 3 2.42 2 f arctgf arctgf arctg πθ θ πθ θ πθ θ = − − + = − + = − +
Acknowledgments
We thank Prof. M. Persico, Prof. H. Jauslin, Dr. S. Guérin, Dr C. Koch, Dr D. Lauvergnat and Dr F.Remacle for helpful discussions. Dr. G. Dive is research associate of the FNRS of Belgium. The computing facilities of IDRIS (Project numbers 061247 and 2006 0811429) as well the financial support of the FNRS in the University of Liège SGI Nic project are gratefully acknowledged. 38 [1] M. Dantus and V. V. Lozovoy, Chem. Rev. , 1813 (2004). [2] M. Shapiro and P. Brumer, Principles of Quantum Control of Molecular Processes (Wiley, New York, 2003). [3] P. Brumer and M. Shapiro, Annu. Rev. Phys. Chem. , 257 (1992); M. Shapiro and P. Brumer, Rep. Prog. Phys. , 859 (2003). [4] D. Tannor and S. A. Rice, J. Chem. Phys. , 5013 (1985); D. Tannor, R. Kosloff, and S. A. Rice, J. Chem. Phys. , 5805 (1986). [5] V. S. Malinovsky, C. Meier and D. J. Tannor, Chem. Phys. , 67 (1997). [6] R. S. Judson and H. Rabitz, Phys. Rev. Lett. , 1500 (1992). [7] W. Zhu, J. Botina, and H. Rabitz, J. Chem. Phys. , 1953 (1998). [8] W. Zhu and H. Rabitz, J. Chem. Phys. , 385 (1998). [9] Y. Ohtsuki, G. Turinici, and H. Rabitz, J. Chem. Phys. , 5509 (2004). [10] U. Gaubatz,P. Rudecki, S. Schiemann, and K. Bergmann, J. Chem. Phys. , 5363 (1990). [11] N.V. Vitanov, T. Halfmann, B. Shore and K. Bergmann, Annu. Rev. Phys. Chem. , 763 (2001). [12] S. Guérin and H. R. Jauslin, Adv. Chem. Phys. , 147 (2003). [13] N. V. Vitanov, K. A. Suominen and B. W.Shore, J. Phys. B , 4535 (1999). [14] K. Hoki, Y. Ohtsuki, and Y. Fujimura, J. Chem. Phys. , 1575 (2001). [15] M. V. Korolkov, J. Manz, and G. K. Paramonov, J. Chem. Phys. , 13927 (1996). [16] H. Umeda, M. Takagi, S. Yamada, S. Koseki and Y. Fujimura, J. Am. Chem. Soc. , 9265 (2002). [17] Y. Zhao and O. Kühn, Chem. Phys. Lett. , 7 (1999). [18] C. M. Tesch, L. Kurtz and R. de Vivie-Riedle, Chem. Phys. Lett. , 633 (2001). 39 [19] C. M. Tesch and R. de Vivie-Riedle, Phys. Rev. Lett. , 157901 (2002). [20] B. M. R. Korff, U. Troppmann, K. L. Kompa and R. de Vivie-Riedle, J. Chem. Phys. , 244509 (2005). [21] D. Babikov, J. Chem. Phys. , 7577 (2004). [22] U. Troppmann and R. de Vivie-Riedle, J. Chem. Phys. , 154105 (2005). [23] T. Cheng and A. Brown, J. Chem. Phys. , 034111 (2006). [24] I. R. Sola, V. S. Malinovsky, J. Santamaria, J. Chem. Phys. , 10955 (2004). [25] B. Lasorne, G. Dive, D. Lauvergnat, and M. Desouter-Lecomte, J. Chem. Phys. , 5831 (2003). [26] T. Taketsugu, N. Tajima, and K. Hirao, J. Chem. Phys. , 1933 (1996); T. Taketsugu and Y. Kumeda, J. Chem. Phys. , 6973 (2001). [27] M. V. Basilevsky, Chem. Phys. , 81 (1977); M. V. Basilevsky, Theor. Chim. Acta , 63 (1987). [28] J. Baker and P. M. W. Gill, J. Comput. Chem. , 465 (1988). [29] M. Hirsch, W. Quapp, and D. Heidrich, Phys. Chem. Chem. Phys. , 5291 (1999); W. Quapp, and V. Melnikov, Phys. Chem. Chem. Phys. , 2735 (2001); W. Quapp, M. Hirsch and D. Heidrich, Theor. Chem. Acc. , 40 (2004). [30] R. Palmeiro, L. M. Frutos and O. Castaño, Int. J. Quantum Chem. , 422 (2002). [31] T. Taketsugu and Y. Kumeda, J. Chem. Phys. , 6973 (2001). [32] B. Lasorne, G. Dive and M. Desouter-Lecomte, J. Chem. Phys. , 184304 (2005). [33] D. Lauvergnat, A. Nauts, J. Chem. Phys. , 8560 (2002). [34] M. D. Feit, J. A. Fleck, and A. Steiger, J. Comput. Phys. , 412 (1982). [35] P. Van Leuven and M. Persico, J. Chem. Phys. , 054319 (2006). [36] F. L. Yip, D. A. Mazziotti, and H. Rabitz, J. Phys. Chem. A , 7264 (2003). [37] S. Gräfe, C. Meier, and V. Engel, J. Chem. Phys. , 184103 (2005). 40 [38] M. Sugawara, J. Chem. Phys. , 6784 (2003). [39] K. Sundermann and R. de Vivie-Riedle, J. Chem. Phys. , 1896 (1999). [40] J. P. Palao and R. Kosloff, Phys. Rev. Lett. , 188301 (2002). [41] C. M. Tesch and R. de Vivie-Riedle, J. Chem. Phys. 121, 12158 (2004). [42] V. Kurkal and S. A. Rice, Chem. Phys. Lett. , 125 (2001). [43] I. Vrabel and W. Jakubetz, J. Chem. Phys. , 7366 (2003). [44] Y Fujimura, L. González, K. Hoki, J. Manz and Y. Ohtsuki, Chem. Phys. Lett. , 1 (1999). [45] M. Amniat-Talab, S. Guérin, N. Sangouard and H. R. Jauslin, Phys. Rev. A , 023805 (2005). [46] T. Wang, M. Kostrun and S. F. Yelin, Phys. Rev. A , 053822 (2004). [47] N. Sangouard, S. Guérin, M. Amniat-Talab and H. R. Jauslin, Phys. Rev. Lett. , 223602 (2004). [48] F. Ticozzi, A. Ferrante and M. Pavon, IEEE Transaction on Automatic Control , 1742 (2004). [49] C. M. Dion, A. Keller and O. Atabek, Eur. Phys. J. D , 249 (2001). [50] D. Sugny, A. Keller, O. Atabek, D. Daems, C. M. Dion, S. Guérin and H. R. Jauslin, Phys. Rev. A , 063402 (2005). [51] A. Matos-Abiague and J. Berakdar, Appl. Phys. Lett. , 2346 (2004). [52] D. Sugny, A. Keller, O. Atabek, D. Daems, S. Guérin and H. R. Jauslin, Phys. Rev. A , 043407 (2004). [53] N. Doslic, O. Kühn, J. Manz and K. Sundermann, J. Phys. Chem A , 9645 (1998). [54] M. Sugawara and Y. Fujimura, J. Chem. Phys. , 5646 (1994). [55] M. A. Nielsen and I. L. Chuang, Quantum computation and Quantum information , Cambridge University Press, Cambridge, U. K. 2000. 41 [56] Z. Kis and F. Renzoni, Phys. Rev. A , 032318 (2002). [57] H. Goto and K. Ichimura, Phys. Rev. A , 012305 (2004). [56] T. Pellizari, S. A. Gardiner, J. I. Cirac and P. Zoller, Phys. Rev. Lett.75