Generalized filtering of laser fields in optimal control theory: application to symmetry filtering of quantum gate operations

We present a modified version of a previously published algorithm (Gollub et al 2008 Phys. Rev. Lett.101 073002) for obtaining an optimized laser field with more general restrictions on the search space of the optimal field. The modification leads to enforcement of the constraints on the optimal field while maintaining good convergence behaviour in most cases. We demonstrate the general applicability of the algorithm by imposing constraints on the temporal symmetry of the optimal fields. The temporal symmetry is used to reduce the number of transitions that have to be optimized for quantum gate operations that involve inversion (NOT gate) or partial inversion (Hadamard gate) of the qubits in a three-dimensional model of ammonia.


Introduction
Optimal control theory (OCT) opens up the ability to manipulate the dynamics of quantum systems, for example by using tailored laser pulses to influence the system dynamics in a desired way [1]- [12]. Processes to be controlled include bond-selective chemistry, optimization of product channels in chemical reactions, or increasing/decreasing strengths of transitions in spectroscopic studies [13]. Most relevant for the present work is the proposal by Tesch and de Vivie-Riedle [14] for the implementation of quantum gate operations using tailored laser pulses and the use of vibrational modes of molecules as qubits-a topic of growing interest [15]- [27].
To achieve the general goal of controlling quantum dynamics, several numerical methods have been deployed such as iterative techniques [1,2,9,12] and genetic algorithms [28]- [30]. While a variety of algorithms are available for theoretical studies, genetic algorithms are primarily used in experiments. The main difficulty within theory, however, arises from the fact most algorithms require repeated and explicit solution of the time-dependent Schrödinger equation for the system. This process can be very time and memory demanding, especially for multi-dimensional systems. Therefore, we use the multi-configuration time-dependent Hartree (MCTDH) method [31]- [34] to perform the propagations. OCT has been previously implemented within the Heidelberg MCTDH package [35]-the so-called OCT-MCTDH approach [26,36,37]-and here we report further extension of this implementation.
Control fields to achieve the desired objective or objectives can readily be found for most problems [38]. However, often the fields obtained with standard OCT using iterative methods exhibit rather complicated temporal structures. In most cases, the complex field arises due to multiple different frequency components within the field 3 . As these complicated fields are not generally accessible within the control experiments, efforts have been made to restrict the search space for the optimal fields, for example by projecting out parts of the wavefunction that cause unwanted frequencies [39], by applying filters [40,41], or utilizing a filtered reference field [42]. Algorithms applying other restrictions on the field, such as constraints on the fluence or on 3 the pulse envelope, have also been applied (see, e.g. [41]). However, many of these iterative algorithms with constraints are not monotonically convergent, which restricts their general utility.
Recently, Gollup et al [40] published a monotonically convergent algorithm based on the standard Krotov algorithm that is able to explicitly suppress unwanted frequencies. In the present contribution, we show that their algorithm actually holds for a much larger class of restrictions on the field. By defining a general projector on a wanted (or unwanted) subspace of all possible fields, we derive a slightly modified algorithm than the one in [40]. The modification accounts for unwanted contributions that are generated because of an approximation during the first few iterations of the algorithm. We demonstrate the general applicability of the algorithm by imposing constraints on the temporal symmetry of the fields when determining optimal pulses for the NOT and Hadamard (HAD) quantum gates. Both, the NOT and especially the HAD gate are used in many quantum computing algorithms as they provide basic logical operations [43]. Time-symmetric optimal pulses have been obtained earlier for NOT gates [44,45]-without any explicit restrictions on the symmetry. This was not completely by chance. Both the NOT and HAD gates are so-called multi-target optimizations on a two-qubit system, i.e. the timedependent Schrödinger equation has to be solved for both states involved. In section 3.2, we show that time-symmetry can (under certain conditions on the states) be utilized to avoid the propagation for one of the states.
The paper is organized as follows. In section 2.1, we review the OCT algorithm while introducing generalized constraints through projection operators. The implementation of the algorithm within the Heidelberg MCTDH software package [35] is discussed in section 2.2. In section 3, we demonstrate the applicability of the algorithm by imposing symmetry constraints on the pulse to optimize quantum gate operations in three-dimensional (3D) ammonia. The 3D ammonia model, including the states defining the qubit basis, is introduced in section 3.1. Section 3.2 provides details on the time-symmetry constraints for the NOT and HAD quantum gates. The results for the NOT and HAD gates are discussed in sections 3.3 and 3.4, respectively. A short summary of the general utility of the method and its specific application to timesymmetry constraints is given in section 4.

OCT
To measure the quality of a given electric field ε(t) in terms of reaching the optimization goal, one introduces the so-called control or cost functional as (h is set to unity throughout the paper): The first term, J 1 , in equation (1) measures the population of the target state that is to be optimized, while the second term, J 2 , punishes for the difference of the optimal field ε from a reference field ε ref .
Here α is the so-called penalty-factor and S(t) later serves as a pulse envelope. The fourth term of equation (1), J 4 , guarantees that the time-dependent Schrödinger equation is fulfilled. The difference to the usual control functional lies in the third term, J 3 , which punishes for contributions the field may contain after operation of a linear operator Q with γ being a Lagrange multiplier. Here we adopted (analogous to the notation used for Hilbert spaces) the use of round brackets for indicating a scalar product in a real linear vector space of functions in the interval [0, T ].
In the following, we will assume that Q is a real valued and symmetric (i.e. Q(t, τ ) = Q(τ, t)) projector onto an 'unwanted' subspace of the complete search space for the optimal field. That is, the field may not contain contributions within the subspace projected onto by Q. This could, for example, be a certain frequency range as in [40], but other subspaces are also possible. If Q is a projector, its orthogonal complement P is given by with the usual relations for projectors: The functional derivative of J with respect to the initial state, the target and the field leads to an expression for the optimal field and equations of motion for the states. The optimal field is obtained as where σ = sign (γ | Q |ε) and The side condition δ J/δγ = 0 can be satisfied by demanding Qε = 0. Then the Lagrange multiplier can be determined by operating Q on both sides of equation (4).
It follows immediately that In other words, the choice of the pulse envelope S(t) influences the calculation of the Lagrange multiplier γ by its transfer properties from the wanted (P) into the unwanted (Q) subspace. To simplify equation (7), one can choose the pulse envelope such that operating S(t) on a given field within the wanted subspace does not lead to contributions outside the wanted subspace, i.e. QSP = 0. This leads to One possible solution of equation (8) can be obtained as Inserting this result into equation (4) leads to The action of Q on C cannot usually be obtained during the current iteration step unless Q is diagonal in the time domain (but in this case the action of Q could have been absorbed in S) and has to be approximated. Following [40], we propose the following iterative scheme to obtain the optimal field during the (n + 1)th forward iteration step: where the approximation is obtained by using the result of C(t) during the previous backward propagation with Equation (11) is a slightly modified version of the algorithm obtained in [40]. The difference lies in operating P on the reference field ε n . This ensures that no unwanted contributions from previous iterations are carried through to later iteration steps. That is, the restrictions on the field are enforced by removing undesired components of the field in the subsequent iteration.
(Alternatively, the reference field could be defined as ε ref = Pε n in which case equation (9) would simplify without affecting the resulting fields.) In fact, Gollub et al showed that the functional, equation (1), is monotonically convergent if γ is obtained exactly. In this case, the action of P leaves the reference field ε n unchanged and equation (11) reduces to the scheme proposed in [40]-provided that the initial guess-field lies completely within the required subspace. It should be noted that the unwanted contributions are typically only generated during the first few iteration steps when the approximation equation (12) is poor and the constraint Qε n+1 = 0 is not fulfilled.

Implementation
We have extended the existing implementation of OCT algorithms within the Heidelberg MCTDH Software Package [26,36,37]. The usual Krotov algorithm [10,11] has been implemented along with optional constraints on the field as discussed in section 2.1. In our implementation, the constraints on the field can be either applied with or without filtering of the reference field. At present, we have implemented projections within the frequency domain, as discussed in [40], and on the temporal symmetry of the field, as discussed in section 3.2. Implementing new constraints into the code is straight forward. In addition, a new program for calculating superpositions of wavefunctions in MCTDH-form is made available. We hope that our code will be made available in one of the future releases of the Heidelberg MCTDH Package.

Model system
For the numerical studies we use a reduced dimensional model of ammonia, where we restrict ourselves to the internal, vibrational degrees of freedom, i.e. global rotations are not taken into account. We use a set of coordinates consisting of the three N-H bond lengths, two bond angles and one dihedral angle as depicted in figure 1. This set of coordinates has been successfully utilized to describe vibrational states for a number of tetra-atomic molecules such as HONO, HFCO, H 2 CS and NH 3 [46]- [48]. The kinetic energy operator for a tetra-atomic using this coordinate system is given in [46], where the inverse-mass matrix has to be replaced by Here m N and m H denote the masses of the nitrogen and hydrogen atoms, respectively. For our calculations, we used accurate potential energy and dipole surfaces that have been published by Yurchenko et al [49,50]. The 3D model used in the following was obtained by setting the three N-H bond lengths to their equilibrium values of 1.01 Å, and then neglecting all terms in the kinetic energy operator, i.e. equation (2) of [46], that contain a derivative with respect to R i . Within this approximation, the three stretching modes present in a full 6D model are removed from the model system. Note that the exact Hamiltonian for a reduced model can be obtained using a procedure discussed by Gatti and Iung [51] but the model Hamiltonian used here provides sufficiently accurate results [36]. Table 1 shows the eigenstates of the model Hamiltonian selected as qubit states. The states were calculated using the improved relaxation [46] option of the MCTDH software, where we used the same setup for the calculations (i.e. grids, basis states and single particle functions) as detailed in [36]. The quantum number ν p 2 labels excitations of the inversion mode, where p indicates the parity. ν l 4 4 labels excitations in the two-fold degenerate bending mode, with l 4 labelling the vibrational angular momentum. In a full 6D model, additional quantum numbers for the stretching modes ν 1 (symmetric stretching) and ν l 3 3 (two-fold degenerate asymmetric stretching) would have to be introduced. Since the dipole matrix element between the qubit states is zero, i.e. 0|µ|1 = 0, due to selection rules, there is no direct transition between |0 and |1 . Therefore, the gate operations will require population transfer through at least one  [54] state outside the qubit system. In such a case, the optimal pulse is expected to contain multiple frequencies. This choice of qubits is not optimal if the goal is solely optimization of quantum gate operations, where one would normally choose a qubit basis with states connected by direct one-photon transitions, but it is useful to test the time-symmetry constraints for a nontrivial control problem.

Temporally symmetric Hamiltonians
In the following, we will demonstrate the general applicability of the algorithm by imposing constraints on the temporal symmetry of the optimal electric field over the time interval [0, T ], i.e. the field is demanded to be of the form Then the projector on the unwanted space can be constructed as which is a projector on all fields with ε(t) = −ε(T − t). In order to fulfil the condition QSP = 0, the pulse envelope is chosen as Temporally symmetric Hamiltonians exhibit a number of interesting features. In particular, one can show by time-slicing the time evolution operator U , that, provided the initial and target states are real (e.g. eigenstates of the system Hamiltonian), applying the time evolution operator generated by a Hamiltonian that is symmetric in time leads to That is, a pulse optimized to transfer population from to automatically performs the inverse with the same amount of population transfer if operated on . This can be interpreted as a quantum NOT gate if for example | labels a qubit |0 and | labels a qubit |1 , where the NOT gate is defined as Moreover, both the initial and the target state will have the same global phase φ at final time T after operation of U , provided that they also had initially the same phase. As mentioned in the introduction, time-symmetric pulses have been observed earlier when optimizing NOT gate operations using standard iterative techniques [44,45]. If there is a direct transition between the two qubits, such that the population transfer is very symmetric, time-symmetric pulses can be favoured as they symmetrically optimize both transitions. However, it must be emphasized that in these previous reports, no explicit time-symmetry constraints were imposed on the fields-the symmetry arose naturally from the optimization process. A similar relation can be found for the HAD gate which is defined as Provided that | +|U |0 | 2 = 1, then Using U 2 |0 = e 2iφ |0 , this yields That is, also for the HAD gate it is sufficient to optimize only one of the operations of equation (20) under the constraint of temporal symmetry of H in the time interval [0, T ]. As in the previous case, automatic phase alignment is guaranteed at the final time T , provided the initial phases are equal and the population transfer is 100% (we observed that the phase alignment can be rapidly destroyed if the population transfer is not very close to 100%). Numerically, the use of the time-symmetry constraint means a reduction in the effort to solve the time-dependent Schrödinger equation as only one of the two states needs to be propagated. However, experimentally, constructing a temporally symmetric Hamiltonian could be a challenging task.

NOT gate
In figure 2, the control functional equation (1) for the optimization of a NOT gate with the transition |0 → |1 is depicted as a function of the number of iterations. The field is constrained to be symmetric as discussed in section 3.2. The iterations were started with the initial guess field where ε 0 = 257.1 MV m −1 , T = 1 ps and ω ≡ 775 cm −1 , which is equivalent to half of the energy difference between the states |0 and |1 . It can clearly be seen that the complete functional J (blue line) is monotonically convergent. However, this is not true for its single components. The curve for the population part of the functional, J 1 (black line), decreases during the first two iterations. Also the filter part, J 3 (green line), increases temporarily between iterations 5 and 10 before converging to zero. After 130 iterations, the change of the components of J has decreased to the order of 10 −5 (after 50 iterations about 10 −3 ) so that the functional can be considered converged. The resulting electric field is depicted in figure 3. It clearly fulfils the constraint condition equation (15) even though it contains multiple frequency components and has a complicated temporal structure. The field contains mainly three frequency components: the frequency ω that was used in the  initial guess field and two additional frequencies at ≈1020 cm −1 and at ≈510 cm −1 . The latter two frequencies correspond approximately to the transition frequencies between the states |0 and |2 − , 0 0 and |2 − , 0 0 and |1 , respectively, where the quantum numbers are used as in table 1. It should be noted that the optimization process started with the same initial conditions but without symmetry constraint or without filtering the reference function does not lead to a symmetric field (not shown here). The resulting population dynamics of both the optimized transition |0 → |1 as well as the transition |1 → |0 , that was not optimized explicitly, are shown in figure 4. The population of the respective target states differs from each other for times smaller than the final time T since the pathways are not symmetric. At the final time, however, both target populations reach approximately 99.9%. The fidelity of the gate operation including both targets-using explicit propagations for both initial states-is obtained as  Here it should be stressed again that this fidelity can only be reached using the method described above if both initial states are real, except for an (identical) initial phase.

HAD gate
In the following, the results for the optimization of the HAD gate using only the transition |0 → |+ under the constraint of a symmetric field are presented. Here we used the same initial guess field equation (23) as for the NOT gate. The control functional is shown in figure 5 as a function of the number of iterations. In the present case, the single components of J show monotonically convergent behaviour, except for the constraint term J 3 (green line), which slightly increases during iteration 2. Compared to the NOT gate, the convergence is faster. This is probably due to the less complicated field as depicted in figure 6. The optimized field mainly contains a single frequency and has an almost Gaussian pulse envelope such that it is not very different from the initial guess field.
The population dynamics of both target states obtained with the optimized field after 130 iterations and using explicit propagations of both initial states is depicted in figure 7. Both target state populations are initially 50% and start oscillating after about 100 fs. The fast oscillations of the target populations result from the different energies of the states |0 and |1 from which the superpositions |+ and |− are constructed. It can be seen that the population dynamics of the two target states are approximately the same after about 800 fs and reach in both cases nearly 100% after 1000 fs. At 1000 fs, the fidelity of the HAD gate operation including both targets-again using explicit propagations for both initial states-has been obtained as

Summary
We have derived a slightly different version of an optimal control algorithm previously proposed by Gollub et al [40], where we have introduced general restrictions on the optimal field using projection operators on desired and undesired contributions. By introducing a different solution of the side condition imposed by a filter, we were able to obtain more rigorously optimized fields with respect to the given constraints. The algorithm proposed in this work accounts for the inaccuracy of the prediction of the Lagrange multiplier during the first iterations where undesired components of the field can be generated and carried through to subsequent iterations.
Here we would like to stress that this is a small, but useful, amendment to the original algorithm proposed by Gollub et al [40]. Importantly, the algorithm has been implemented in the broadly applicable MCTDH software package [35].
We have demonstrated the applicability of the algorithm by imposing symmetry constraints on the field. Precisely constrained fields are especially crucial for optimizing quantum gate operations using symmetry conditions in the time domain. We have used the algorithm to optimize both targets for the NOT and HAD quantum gates while solving the time-dependent Schrödinger equation for only one of the two involved states. The propagation of the second state has been replaced by imposing constraints on the field and demanding the initial and target states to be real (except for an identical global phase), which automatically guaranteed that also the second target is optimized as well.