Simulating quantum search algorithm using vibronic states of I2 manipulated by optimally designed gate pulses

In this paper, molecular quantum computation is numerically studied with the quantum search algorithm (Grover's algorithm) by means of optimal control simulation. Qubits are implemented in the vibronic states of I2, while gate operations are realized by optimally designed laser pulses. The methodological aspects of the simulation are discussed in detail. We show that the algorithm for solving a gate pulse-design problem has the same mathematical form as a state-to-state control problem in the density matrix formalism, which provides monotonically convergent algorithms as an alternative to the Krotov method. The sequential irradiation of separately designed gate pulses leads to the population distribution predicted by Grover's algorithm. The computational accuracy is reduced by the imperfect quality of the pulse design and by the electronic decoherence processes that are modeled by the non-Markovian master equation. However, as long as we focus on the population distribution of the vibronic qubits, we can search a target state with high probability without introducing error-correction processes during the computation. A generalized gate pulse-design scheme to explicitly include decoherence effects is outlined, in which we propose a new objective functional together with its solution algorithm that guarantees monotonic convergence.


Introduction
Coherent control of molecular dynamics is initiated to manipulate quantum interferences induced by (shaped) laser pulses, thereby improving selectivity and energetic efficiency of chemical reactions [1]- [7]. Its principles are common to those of other research areas in which coherent control is essential, such as quantum information processing [8]. As there are huge numbers of states in molecules, it would be natural to consider the possibility of using them as quantum information resources (qubits) [9]. For example, if we consider q vibrational degrees of freedom in a molecule, they can be regarded as q qubits. With sophisticated-design laser pulses that achieve gate operations, the molecular states satisfy, at least in principle, the conditions required to be qubits, e.g. those discussed by DiVincenzo and Loss [10]. From the point of view of mathematical isomorphism, the direct product space spanned by the q vibrational qubits can be alternatively realized by the 2 q states in a specified vibrational mode, which is usually called a unary representation of states. Although the unary representation scheme is known to require exponential physical cost [11] (i.e. not scalable in the sense of [10]), it could provide a practical way especially for small-scale molecular computation, because each vibrational degree of freedom contains a number of eigenstates. In the present paper, we therefore study the unary representation approach through a case study of vibronic states of I 2 .
There are several experimental and theoretical proposals that demonstrate the usefulness of molecular qubits [12]- [38]. Those proposals can be classified into two categories: weakfield and strong-field schemes. The weak-field schemes [12]- [20] are based on the perturbative interactions between laser pulses and molecules, which result in the limited number of gate operations. Qubits are defined by using molecular eigenstates [12]- [17] as well as by their superposition states [18], and gate operations are realized by the free propagation of wave packets and/or optical transitions controlled by shaped laser pulses. The main advantage is that all the basic ingredients required for computation, i.e. the initial preparation, gate operation and detection step, can be achieved by experimentally feasible spectroscopic methods, such as fourwave mixing (FWM) with different pulse sequences [13]- [16], multi-photon ionization [17] and quantum interferometry [12], [18]- [20].
In the strong-field schemes [21]- [38], laser pulses are designed in a non-perturbative manner so that they act as gate operators. At present, because of the challenges posed by actual experiments, many of the schemes are investigated theoretically and numerically, often by 3 using optimal control simulation [3,4,39]. Some of these studies adopt parameter-optimization approaches [21], whereas others employ full optimization approaches with different criteria to evaluate gate fidelities [22]- [26], [28]- [38]. The performance of the designed gate pulses is often evaluated through case studies of the simulation of the Deutsch-Jozsa algorithm [29,30,37], quantum Fourier transform [25,26], etc.
In the present paper, we study the quantum search algorithm (Grover's algorithm) [40] by combining vibronic qubits (unary representation) and optimally designed gate pulses. So far, the quantum search algorithm has not attracted much attention to discuss the possibility of molecular quantum computation although it only needs the superposition principle [41] and therefore we can use an ensemble of molecules. In the present numerical study, we show that sequential irradiation of separately designed gate pulses leads to the population distribution predicted by the quantum search algorithm. In section 2, we briefly summarize Grover's algorithm to explain our settings, and examine the optimal control procedure to design gate pulses and their solution algorithms. We discuss the methodological aspects of a gate pulsedesign problem [25,26], in which we show that it has the same mathematical structure as a stateto-state control problem in the density-matrix formalism. We thus derive alternative algorithms for solving the gate pulse-design equations according to our previous study [42]. The numerical results are shown in section 3, in which we examine how the lack of extremely high precision in the design of gate pulses and the presence of decoherence result in reduced accuracy of computational performance. Finally, we outline a generalized gate pulse-design procedure to explicitly include decoherence effects. [8,40]: setting the stage We briefly summarize the quantum search problem considered here. Our purpose is to find an unknown integer number, a, between 0 and N − 1 with N = 2 q , in which each integer x ∈ [0, N − 1] is associated with a computational basis state (q-qubit input register), |x . We introduce a quantum oracle, V , which can recognize solutions to the search problem without explicitly specifying its structure. The oracle is defined by a unitary operator that solely changes the phase of the target state, |a , associated with the unknown (target) number:

Quantum search algorithm
In other words, our purpose is to find the special state involved in the oracle. If we execute the search with the oracle on a classical computer, we can find the unknown number with 50% probability after half of N trials. However, the quantum search algorithm (Grover's algorithm) can find the unknown number with virtually 100% probability after √ N trials [8,40]. As usual, we initially transform the q-qubit 'standard' state, |0 , into the equal superposition state, where H is the Hadamard transform. Then, we apply sequentially the oracle, the Hadamard transform, the conditional phase shift and the Hadamard transform. The set of these four 4 operations is called the Grover iteration, which we denote by . The last three operations are often brought together to form one unitary operation, W , called the inversion about mean operation, which is defined by After the jth Grover iteration, the initially prepared state is transformed into where the state |a ⊥ = 1 x=0(x =a) |x is orthogonal to the state |a . The angle, θ, in equation (4) is given by sin θ the satisfaction of the condition, ( j + 1 2 )θ = π/2, leads to the solution that requires only O( In the present study, we consider a two-qubit search space as the minimal model, in which the angle is given by θ = π/3. After the jth Grover iteration, we find the unknown number a = 0, 1, 2 or 3 with probability P ( j) a = | a| j |ψ | 2 , which has a value of 1 = P (1) a = P (4) a and that of 0. 25 As the essential operations in the Grover iteration are the oracle and the inversion about mean, we regard them as the basic operations and design two optimal pulses that achieve those operations. Then, we apply the independently designed operation pulses to the qubits implemented in the vibronic states of I 2 to evaluate how they work.

Optimal gate pulse design
The Hamiltonian of a molecule interacting with a time-dependent electric field, E(t), is given by where H 0 and µ are the field-free molecular Hamiltonian and the electric dipole moment operator, respectively. The time evolution operator is expressed as where the arrow denotes the time ordering. If the laser pulse realizes a unitary operation, ( = V, W ), at a specified time, t f , then we have for an arbitrary set of basis functions {|n } that span the Hilbert space, H q , associated with qubits. In fact, equation (7) guarantees that U (t f , 0)|ψ = |ψ for an arbitrary qubit state, |ψ = n∈H q |n C n . When we design a laser pulse that realizes the gate operator, , by means of optimal control simulation, we introduce a quantitative criterion called an objective functional to evaluate the difference between an actual state, U (t f , 0)|n , and an ideal state, |n . One possible evaluation is to measure the norm of the vector of the difference between these states, which leads to the functional 5 Then the optimal pulse can be defined by the pulse that minimizes the objective functional. If we remove constant terms that do not contribute to the minimization, the problem ends up with the maximization of the following functional: which was first used by Palao and Kosloff [25,26]. In equation (9), P is the projector onto the subspace, H q . Another possible evaluation is to measure the absolute square of the overlap between the states, U (t f , 0)|n and |n , which defines the functional This was used by Rangan and Bucksbaum [22] and Tesch and de Vivie-Riedle [23,24] and is called a (simultaneous) state-to-state control scheme depending on the number of basis functions. When adopting the functional in equation (10), it is known that we have to be careful in choosing the basis functions. For example, suppose that the designed gate pulse introduces a state-dependent (erroneous) phase, ϕ n , in the operation such that Although the extra phase does not change the value of the functional in equation (10), it is obvious that U (t f , 0)|ψ = |ψ for an arbitrary qubit state, |ψ . This feature has been pointed out by Palao and Kosloff [26] in a different and more general manner and has been numerically shown by Zhao and Babikov [35] through a case study. We therefore adopt the functional in equation (9) to design gate pulses.
Here, we would like to point out that the gate pulse-design problem in equation (9) has the same mathematical form as a state-to-state control problem within the density-matrix formalism [42]. To show this explicitly, we introduce the double-space notation [42]- [44], whereby the functional in equation (9) is expressed as The scalar product of two operators, O 1 and O 2 , is defined by . The timeevolution operator obeys the following equation of motion: The initial condition is given by |U (0, 0)P = |P = n∈H q |nn , in which the doublespace vector |mn corresponds to |m n|. The commutators of H t , H 0 and µ are expressed as where the cap and tilde denote the right-hand side and left-hand side acting operators. For example, If we introduce the penalty term to minimize pulse fluence, the optimal gate pulses are designed to maximize the following objective functional: where a positive function, A(t), weighs the physical significance of the penalty. This objective functional has the same form as a state-to-state control within the density-matrix formalism [42]. 6 Applying calculus of variations to equation (14) under the constraint of the equation of motion in equation (13), we obtain the gate pulse-design equations. The optimal gate pulse is expressed as where | (t) is the Lagrange multiplier associated with the constraint of equation (13) and obeys the equation of motion with a final condition, | (t f ) = | .

Solution algorithms
As shown in the previous subsection, the double-space notation shows that the optimal gate pulse-design problem in equation (9) has the same mathematical form as a state-to-state control problem. We can solve the gate pulse-design equations according to the previously developed solution algorithm in the density-matrix formalism, which guarantees the monotonic convergence [42], [44]- [46]. The kth iteration step is summarized as follows: with the final condition, | (k) (t f ) = | , and with the initial condition, |U (k) (t = 0, 0)P = |P , in which the electric fields are given bȳ and If we ignore the intermediate improvement of the electric field in equation (19), the above algorithm is reduced to the Krotov method [26,47]. In appendix A, we give a matrix representation of the solution algorithm.

Results and discussion
We consider the vibrational states in the X-electronic state, {|Xu X }, and those in the B state, {|Bv B }, of I 2 , in which the four vibronic states, |Bv B = 28 , |Bv B = 29 , |Bv B = 30 and |Bv B = 31 , are assigned to the logical states |00 , |01 , |10 and |11 (two qubits), respectively. The other vibronic states are regarded as auxiliary states. The gate operations are realized through the optical transitions between the two electronic states. The potential curves are obtained by the interpolation of the Rydberg-Klein-Rees (RKR) points given in [48] and [49]. The vibrational states in each electronic state are calculated by a conventional diagonalization 7 procedure. The wave function is expanded in terms of these vibronic states, and its time evolution is calculated by the Runge-Kutta method, in which optical transition matrix elements are obtained by using the electronic transition moment function given in [50]. According to the optimal control procedure in section 2.3, we numerically design the gate pulses with a temporal width of t f = 1000 fs. These pulses achieve the gate operations through resonant optical transitions between the X-and B-electronic states. As the transition frequencies are much greater than the Rabi frequencies induced by the gate pulses, we can safely assume the rotating-wave approximation (RWA). In this case, we can freely choose the value of the electronic excitation energy, say 7000 cm −1 , without changing the dynamics. The optimal gate pulses are designed by using 27 states, {|Xu X , u X = 0-14} and {|Bv B , v B = 24-35}. The calculated pulses are substituted into the equation of motion with a larger basis to check the convergence with respect to the number of vibronic states. Except for the optimal pulses themselves, all the numerical results are obtained by the calculations with the larger vibronic states. In the penalty term in equation (14), the positive function is assumed to have the form of A(t) = A 0 φ(t), in which φ(t) is given by a sine function with a period of 4 T (T = 200 fs) for 0 t T and t f − T t t f , whereas φ(t) = 1 for T t t f − T . First, we design the gate pulses assuming a small value of A 0 . Then, we restart the calculations by employing a larger value of A 0 , in which the optimal pulse obtained in the previous stage is used as the initial trial field. We repeat the procedure until the gate fidelity virtually does not depend on the choice of the value of A 0 .
As we have mentioned in the last paragraph of section 2.1, we optimally design the gate pulses associated with the oracle and the inversion about mean operations. Figure 1 shows the sets of these designed gate pulses that realize the Grover iteration. Each pulse has a temporal width of t f = 1000 fs. The first halves of the pulses, t ∈ [0, 1000 fs], correspond to (panel (a)) the oracle pulse associated with the target state, |00 , and (panel (b)) that with |10 . The second halves of the pulses, t ∈ [1000 fs, 2000 fs], correspond to the inversion about mean operation that is independent of the target states. If we evaluate the fidelity based on the control achievement, F PK = 2Re tr{ † U (t f , 0)P}, we will have higher than 99% fidelity for the oracle operations and 94% for the inversion about mean operation. The difference in fidelity originates from the different structures of the gate operations. The oracle operations are expressed by diagonal matrices with three +1 elements and one −1 element. The element, −1, specifies the target number. On the other hand, a matrix that represents the inversion about mean operation does not have such a simple structure. This results in difficulty in designing the pulse, leading to a decrease in fidelity. The fidelity can be improved by assuming a longer control time and/or could be done by introducing intermediate targets to avoid undesirable transitions [51,52]. In the following analyses, however, we use the gate pulses designed here and examine how the lack of extreme precision affects computation accuracy. Figure 2 shows the power spectra of these gate pulses. As we have assumed an artificial value for the electronic transition frequency, the abscissa shows a relative frequency. As can be seen from figure 2(c), the pulse associated with the inversion about mean operation has a much more complicated structure than those of oracle operations, which also explains its low fidelity (94%).
To see the performance of these gate pulses, we solve the time-dependent Schrödinger equation under the sequential irradiation of the gate pulses assuming the initial condition, |ψ(t = 0) = 1 at 4000, 6000, 10 000 and 12 000 fs. As shown in figure 3, the results are in good agreement with the prediction. In the present example, the fidelity of the gate pulse associated with the inversion about mean operation is 94%. The lack of extremely high precision could lead to poor performance in finding the target state. We show such an example in figure 4, in which |11 is chosen as a target number. As shown in figure 4(a), the overall behavior of the time-dependent population of the target state appears similar to that in figure 3; that is, the target population reaches a maximum value at 2000 and 8000 fs; however, those values are ∼60% although they should be 100% according to Grover's algorithm. We attribute the reduction in probability of obtaining the target to the reduction in the total population within the qubit space, not to the distortion of the wave packet within the qubit space. In fact, the thin (green) line shows the time-dependent population of the excited electronic state 'outside' the qubit space. It always has a value typically larger than 40-50%. Our explanation is confirmed by the result in figure 4(b) that shows the timedependent population normalized with respect to the total qubit population. The normalized population well reproduces the predicted time-dependent behavior at each Grover iteration. That is, we can find the target state with quite a high probability as long as we focus on the vibronic states associated with the qubits without introducing error-correction processes during the computation. This suggests that a reduction in the target population does not necessarily lead to a reduction in the probability of finding the target state. Next we show that a similar result is numerically observed for the quantum search in the presence of electronic decoherence. For this purpose, we use the set of gate pulses in figure 1(a), which are designed in the absence of decoherence. We apply it to the qubit state that is initially in the state, |ψ(t = 0) = 1 2 31 v B =28 |Bv B , and calculate the target population after the irradiation of the set of gate pulses in figure 1(a). To systematically examine the decoherence effects on the quantum search, the time evolution of the vibronic system is assumed to be described by the non-Markovian master equation (in the double-space notation) [43,53]: where |ρ(t) is the density operator of the vibronic system, and L t corresponds to the commutator, [H t , · · ·], that describes the unitary time evolution. In equation (21), the electronic decoherence is characterized by the coupling strength, γ , and the reservoir correlation time, τ . In the following examples, we assume artificially large values for γ to clearly see the decoherence effects within our short simulation time. Figure 5(a) shows the memory kernel, exp (−t/τ ), as a function of time for several values of τ . Figure 5 to the results in the cases of τ = ∞ (no decoherence), τ = 1000, 300, 100 and 0 fs (Markovian decoherence) from the top to the bottom. As expected, the target population after the first Grover iteration decreases with 1/τ . Figure 6 summarizes the target population as a function of γ (cm −1 ) and τ (fs). As the electronic decoherence suppresses desirable optical transitions, it leads to a reduction in the accuracy of the gate operations. In fact, the smaller coupling parameter and the longer reservoir correlation time yield a higher probability of finding the target state. Figure 7 shows the target populations normalized with respect to the populations within the qubit space. Even in the Markovian-decoherence case with a large coupling parameter, γ = 9 cm −1 , the normalized population has a value of 0.76. This means that if we focus on the populations within the qubit space, we can find the target state with high probability even in the presence of electronic decoherence.
If the designed pulses achieve the gate operations with 100% probability, the change in the relative phase between them does not change the computational results as long as the RWA is valid. However, because of the lack of 100% fidelity, here we examine how the computational results are affected by the relative phases. For this purpose, we calculate the population of the target state, |00 , after the first Grover iteration as a function of the relative phase, which is shown in figure 8. We see that both the population and the normalized population are virtually constant with the relative phases, being independent of the presence/absence of  electronic decoherence. This indicates that the computation is robust against the relative phases between gate pulses even in the presence of imperfection in the pulse design and/or electronic decoherence.
Finally, we outline a generalized optimal pulse-design procedure to explicitly include (non-Markovian) decoherence effects. For convenience, we rewrite equation (21) in terms of where In order to explicitly include decoherence effects in the optimal gate pulse design, we propose a new functional where the operator, ⊗ =ˆ ˜ , represents the double-space gate operator. According to this procedure, an optimal gate pulse is defined by the pulse that maximizes the functional in equation (24). As shown in appendix B, there exists a monotonically convergent algorithm to solve the optimal control problem in equation (24). If we consider the special case in which the decoherence is neglected in equation (22), then the time-evolution operator is expressed as the product of the cap-and tilde-space timeevolution operators associated with the right-hand side and left-hand side acting Hamiltonians, 0). Ignoring the constant terms in equation (24), the problem is reduced to the maximization of the following functional: where tr ⊗ {· · ·} = m,n mn|{· · ·}|mn and P ⊗ =PP are the double-space trace and projector, respectively. In this special case of no decoherence, the functional in equation (24) is reduced to that introduced by Palao and Kosloff [26]. However, we would like to emphasize that in the presence of decoherence, we have to explicitly deal with the functional in equation (24). Its numerical applications will be shown elsewhere [54].

Summary
By means of optimal control simulation, we have numerically studied the possibility of molecular quantum computation through the case study of Grover's algorithm. Qubits are implemented in the vibronic states of I 2 (unary representation), and basic quantum operations are realized by optimally designed laser pulses. We have discussed how the computational accuracy is reduced by the imperfect quality of designed gate pulses and by the electronic decoherence processes that are described by the non-Markovian master equation with a wide range of values of coupling strength and memory time. In all the cases considered here, the computational accuracy is considerably improved by normalizing the superposition states of the qubits with respect to the qubit population without introducing error-correction processes. That is, as long as we focus on the population distribution of the vibronic qubits, we can search a target state with high probability. This indicates that there could be some degree of robustness in the 'shape' of a qubit-wave packet against undesirable processes. We have discussed the methodological details of gate pulse design. In the absence of decoherence, as the pulse-design procedure used by Palao and Kosloff [25,26] has the same mathematical structure as a state-to-state control problem in the density matrix formalism, we have derived the alternative solution algorithm according to our previous procedure [42]. To explicitly include decoherence effects, we have outlined a generalized gate pulse-design scheme, in which we proposed a new objective functional together with its solution algorithm that guarantees monotonically convergent behavior.

Appendix A. Matrix representation of the solution algorithm
We derive the matrix representation of equations (17)- (20). The formal solutions to equations (17) and (18) are given by and respectively, where superscripts,k and k, indicate the time evolution under the electric fields, E (k) (t) and E (k) (t), respectively. Using these formal solutions, equation (20) can be rewritten as We can derive a similar expression forĒ (k) (t) given in equation (19). If we introduce the state vectors such that then the solution algorithm in the matrix form is summarized as follows: ih ∂ ∂t |λ Following the general conclusions [42,[44][45][46], it is apparent that the algorithm guarantees monotonic convergence. and the Lagrange multiplier, | mn (t) , which represents the constraint due to equation (B.2), then the solution algorithm at the kth iteration step is summarized as follows: with a final condition, | (k) mn (t f ) = |ρ (k−1) mn (t f ) + ⊗ |mn , and ih ∂ ∂t |ρ To prove the monotonically convergent behavior of the algorithm, it is convenient to introduce another state vector, |σ (k) mn (t) = |ρ (k) mn (t) + ⊗ |mn . We consider the difference in objective functional between the kth and (k−1)th adjacent iteration steps, δ J (k, k−1) where |δσ (k) mn (t) = |σ (k) mn (t) − |σ (k−1) mn (t) . We then introduce a function defined by S(t) = 2Re m,n∈H q mn (t f )|δσ (k,k−1) which proves the monotonic convergence behavior.