Pulse based Variational Quantum Optimal Control for hybrid quantum computing

This work studies pulse based variational quantum algorithms (VQAs), which are designed to determine the ground state of a quantum mechanical system by combining classical and quantum hardware. In contrast to more standard gate based methods, pulse based methods aim to directly optimize the laser pulses interacting with the qubits, instead of using some parametrized gate based circuit. Using the mathematical formalism of optimal control, these laser pulses are optimized. This method has been used in quantum computing to optimize pulses for quantum gate implementations, but has only recently been proposed for full optimization in VQAs. Pulse based methods have several advantages over gate based methods such as faster state preparation, simpler implementation and more freedom in moving through the state space. Based on these ideas, we present the development of a novel adjoint based variational method. This method can be tailored towards and applied in neutral atom quantum computers. This method of pulse based variational quantum optimal control is able to approximate molecular ground states of simple molecules up to chemical accuracy and is able to compete with the gate based variational quantum eigensolver in terms of total number of quantum evaluations. The total evolution time $T$ and the form of the control Hamiltonian $H_c$ are important factors in the convergence behavior to the ground state energy, both having influence on the quantum speed limit and the controllability of the system.


I. INTRODUCTION
Quantum computing is presently in the noisy intermediate-scale quantum (NISQ) era [4], where the available quantum computers cannot outperform their classical counterparts. Nevertheless, even in the NISQ era, quantum computers can be used for specific, well-designed cases [5]. One such application is in solving optimization tasks. In this field, Variational Quantum Algorithms (VQAs) have emerged as the leading optimization technique to obtain quantum advantage on NISQ devices. These algorithms are hybrid, i.e. the algorithm exploits both the classical and quantum computer. Certain VQAs have shown proof of concept for small dimensional problems with various designs of qubits [6][7][8][9][10]. The common aim of these algorithms is to determine the lowest eigenvalue of a quantum Hamiltonian H mol , which can be equated to finding the ground state energy of a molecule. One example of a VQA is the variational quantum eigensolver (VQE), which evolves an initial state through a gate-based quantum logic circuit to prepare an approximate ground state |ψ f . In recent years, VQE has been implemented on many different qubit systems and has become a highly active area of research (cf. [11,12] for recent comprehensive overviews).
Generally, in VQAs, the final qubit states |ψ f are prepared from an initial state |ψ 0 by a quantum circuit of parametrized gates. These gates can either manipulate a single qubit or multiple qubits at once, and particularly * Corresponding author: r.j.p.t.d.keijzer@tue.nl on a neutral atom qubit system, are implemented using specific laser pulses. The pulses necessary to implement such a gate can be optimized in terms of fidelity, robustness, etc. using the mathematical formalism of optimal control [13][14][15]. This theory finds its use in quantum algorithms such as GRAPE [16] and CRAB [17]. The entire gate sequence can in principle be seen as a discretization of one continuous laser pulse, as thoroughly explored by Magann et al. [1]. From this stems our idea of Variational Quantum Optimal Control (VQOC) as an alternative to VQE. Instead of optimizing a sequence of gate parameters, a laser pulse is optimized to construct a state |ψ f close to the ground state of a Hamiltonian H mol . This idea has been principally investigated by Meitei et al. [2] and Choquette et al. [3]. The aim of our work is to explore this option more deeply and systematically by developing a theory for VQOC, and to compare our results with the more established VQE method. The main advantages of our pulse-based method are that the pulses provide more flexibility in the evolution of the initial state, resulting in a broader exploration of the state space and allowing for the total evolution time T to be minimized. Both of these advantages mitigate unwanted effects of qubit lifetimes and decoherence, and are therefore especially important in the NISQ era. We further demonstrate how to realistically implement VQOC on a Rydberg atom-based qubit system. From this analysis, we see that VQOC can approximate ground state energies of small molecules beyond chemical accuracy [18] on a Rydberg system. We note that our proposed VQOC algorithm bares similarity to the ctrl-VQE algorithm by Meitei et al. [2], e.g. pulse-based controls with piecewise constant pulses. However, contrary to ctrl-VQE, the adjoint based ap-proach used in VQOC results in an exact expression of a continuous time gradient for the cost functional with respect to the pulses, thus allowing to extend the method beyond square pulses, e.g. sum of Gaussians. Furthermore, this approach allows us to prove minimizer existence results, solidifying the well-posedness of the control problem.
The layout of this paper is as follows. Sec. II describes VQAs in the NISQ era. Sec. III explains the standard VQE method and its shortcomings. In Sec. IV our novel VQOC method is rigorously introduced. In Sec. V the computational cost on the quantum processing unit of both algorithms is described, as well as our methods for comparing the algorithms. Finally, Sec. VII shows initial results of VQOC, and compares the performance of VQOC to that of VQE.

II. VQAS IN THE NISQ ERA
The aim of a variational quantum algorithm is to evolve an initial qubit state |ψ 0 ∈ H m := span({|0 , |1 }) ⊗m C 2 m in a given total evolution time T , to the ground state |ψ g of a problem Hamiltonian H mol ∈ L(H m ). Here, H m is the Hilbert space of m-qubit states, and L(H m ) is the Hilbert space of operators on m-qubits. The Frobenius inner product on L(H m ) is given by The state evolution, as with every quantum system, is mediated through the Schrödinger equation, given in propagator formalism as where U (t) ∈ L(H m ) is the unitary propagator which describes the evolution of the qubit state as |ψ(t) = U (t)|ψ 0 [19]. The goal of a VQA is thus to manipulate H(t) in a way that makes U (T ) map the initial state |ψ 0 to the ground state |ψ g . To do so, the Hamiltonian H(t) is parametrized, and the variational principle [20] is employed to approximate the ground state. Two important factors in this process are the possible parametrization of the Hamiltonian H(t) and the total evolution time T . Figure 1 illustrates the Hilbert space of qubit states.
The possible Hamiltonians H(t) determine a subset of the Hilbert space that is reachable, e.g. if the initial state is unentangled and there are no interaction terms in H(t), then all entangled states are unreachable. An even smaller subset is the set of states reachable in a finite time T (given a bound on the infinity norm of H(t)). The concept of a minimal time necessary to reach a state |ψ f from an initial |ψ 0 is called the quantum speed limit (QSL) [21][22][23]. Several factors are important in creating a NISQ-friendly VQA: (A) the total evolution time T should be as small as possible to suppress decoherence effects, (B) a large portion of Hilbert space should be expressible, so the algorithm can be used to find the ground states of many H mol , and (C) the evolution should be straightforward to implement and control on the quantum computing system. We will argue in Secs. III and IV that our VQOC method adheres to these criteria better than VQE does.

III. VQE
The variational quantum eigensolver, first proposed in a paper by Peruzzo and McClean in 2014 [24], is the most widely implemented VQA to date. Using a parametrized quantum logic circuit, a sequence of quantum gates, it has proven able to closely approximate the ground state of several small molecules [7,8,11]. See also cf. [25] for a recent overview of VQE.
In the VQE algorithm, the evolution of the initial qubit state is done via qubit gates. These gates can be described as matrix elements L(H n ), where 1 ≤ n ≤ m is the number of qubits the gate acts on. Examples of these gates are R X (θ x ) and R Z (θ z ) gates, which respectively rotate a qubit around the x-axis by an angle θ x and z-axis by an angle θ z on the Bloch sphere (see Fig. 3). Other examples include entangling gates such as the canonical CNOT gate and the SWAP gate. A sequence of parametrized and unparametrized gates together with a set of parameters θ determines the evolution.
One commonly used template for such a gate sequence alternates between blocks of parametrized single qubit gates U i,j , with i indicating the qubit and j the block, and an unparametrized entangling gate U ent (which itself could consist of multiple small gates, such as CNOT gates). The final propagator will take the form where the depth d of a state preparation is defined as the number of blocks in the gate sequence. This method is called the hardware-efficient ansatz [26] and is especially relevant in NISQ machines where single qubit gates U q,d can be implemented with high fidelity, and the passive interaction between the qubits can be used to create entanglement gates U ent [4]. If such an interaction is not present and the qubit system possesses some form of control in entanglement operations, then there is more freedom in choosing what U ent should look like, and this choice is highly influential on the performance of the algorithm [27].
After state preparation, the energy of the prepared state on H mol is determined, which can be done efficiently on a quantum computer for Hamiltonians with sparse decompositions [28]. Based on these measurement outcomes, the set of parameters is updated in order to create a state with lower energy. In VQE, stochastic approaches are often implemented to update the parameters. One commonly used example is the simultaneous perturbation stochastic approximation (SPSA) method [6,29,30]. This method samples two points in the parameter space close to θ and measures their energies. Based on these results, a step is taken against the direction of an estimated gradient of the energy w.r.t. the parameters. This work will employ the SPSA method for VQE as, similar to our own VQOC method, it is a first order method. The total number of quantum evaluations (QE) per parameter update is thus given by where #shots = O( −2 ) is the number of times a measurement is to be repeated in order to get an O( ) error in the expectation values. The product of the number of iterations and the number of quantum evaluations per iteration is thus the total number of times a quantum circuit is run on the system. In the NISQ era, where quantum resources are limited, it is important to keep this number as small as possible while desiring a given accuracy. For SPSA-based VQE, two expectation values are necessary per parameter update, hence the factor 2 (see [31] for more detail). In this work, expectation values are calculated as products between matrices and vectors instead of estimated via simulated measurements.
When looking at the NISQ-friendly criteria of Sec. II we see that for VQE, the following issues exist: (A) The total evolution time T is fully determined by the gate sequence and can not be decreased unless the gate sequence is altered (see Fig. 3).
(B) Even with enough control available, the gate sequence might inhibit the space of reachable states. For example, if the qubits can be fully rotated on the Bloch sphere by the quantum computing system; if only Z-gates were included in the gate sequence, a large part of the Hilbert space becomes inaccessible.
(C) The necessary control to create certain gates might be demanding. For instance, in a Rydberg system, complex and time-consuming laser pulses might be necessary to create specific single qubit or entanglement gates. Implementation is, therefore, not straightforward.
Based on these shortcomings in VQE, we believe that VQOC can serve as a more NISQ-friendly VQA. Figure 2: Schematic diagram of VQE illustrating the interplay between quantum space search and classical optimization. An initial state evolves through a hardwareefficient ansatz gate sequence of parametrized gates U i,j and entanglement gates U ent . Afterwards, the state is post-rotated (PR) [32] and measured. A classical optimization algorithm then updates the parameters of the U i,j gates. The northern hemisphere state (red dot) needs to be mapped to the southern hemisphere state (red dot). By taking a R Z gate followed by a R Y gate, the green path is taken. In terms of technical control, the blue path could have been followed, resulting in faster evolution.

IV. OPTIMAL CONTROL
In this section, the derivation of the pulse-based variational quantum optimal control algorithm (VQOC) as a VQA is given. The motivation for a pulse-based method stems from the idea that every gate is implemented as a laser pulse. Therefore, the entire gate sequence can be seen as a discretization of one continuous laser pulse, as principally explored by Magann et al. [1]. On the other hand, the origin of gate-based quantum logic circuits stems from an analogy to digital computing. In principle, however, there is nothing digital about the outset of VQAs. Therefore, there is also no a priori reason as to why gate-based algorithms would be preferable.
Quantum optimal control theory has been connected to variational state preparation in several works [1,2,33] over the last few years. Recent work has shown that pulse-based methods can indeed prepare desired states in the minimum evolution time set by the quantum speed limit [33]. Furthermore, constructing specific control pulses for gate designs focusing on high-fidelity [34][35][36] and robust [37,38] gates have been a lively topic of research.
Our VQOC algorithm can be used to optimize various types of control functions in the qubit system. In particular, this algorithm can optimize the controllable laser pulses that interact with individual qubits. Because of this, we refer to the control functions as pulses. The novelty in our approach lies in the use of an adjoint equation in computing the gradient of the energy with respect to the pulses, which provides a descent direction in the parameter space. A similar approach has been taken by Zhu et al. [39] to specifically optimize electric fields for atom state preparation. Furthermore, we propose a way to compare the results of pulse-based algorithms to gatebased algorithms in terms of energy convergence and time spent calculating on the quantum system, as described in Sec. VI.
To find the ground state energy of a quantum Hamiltonian H mol with the optimal control formalism, we will need to define the set of admissible pulses. Let The space of admissible pulses Z ad ensures that the control pulses remain physically realizable. The inner prod-uct in Z is given by with z being the complex conjugate of z.
We then consider a parametrization of the Hamiltonian where H d is a drift Hamiltonian representing the passive evolution of the system (see Eq. (7) below), and The goal of our optimal control problem setup is then to minimize a given real-valued cost functional J = J(U, z), where U = U (t) are time-dependent unitary propagators satisfying the Schrödinger equation (1) with prescribed pulses z ∈ Z ad . For an extensive description of our optimal control theory, we refer to the Supplementary Material. For the ground state energy problem, we consider the cost functional J(U, z) = J 1 (U ) + J 2 (z), with J 1 and J 2 given by where |ψ 0 is a fixed initial qubit state. The functional J 1 , therefore, describes the energy of the state generated by the unitary U at time t = T , while the J 2 term punishes for high amplitude pulses.
Altogether, we have a well-defined ground state energy minimization problem, which is shown to admit a minimizer z * ∈ Z ad in the Supplementary Material. In practice, the parameter λ > 0 in J 2 is chosen to be large, in which case the constraint on the pulses does not take effect. In this case, the optimality conditions for the minimization problem read as follows (the derivations are provided in the Supplementary Material): It can be verified that the solution P ∈ X to the adjoint equation is given explicitly in terms of U as which is fully dependent on U and the initial parameters. This allows for the second term in the control equation to be expressed as where Γ(t, s) := U (t)U † (s) describes the evolution by the pulses from s to t, [A, B] = AB − BA is the usual commutator, and V k,l are unitaries decomposing Q l as where K ∈ N is the necessary number of unitaries. For practical purposes, where the control terms work on only 1 (or occasionally 2) qubits K = O(1). These terms can be efficiently determined on a quantum computer by first applying the pulse until time t, performing a single gate operation based on the parameter gradient algorithm [40] (see Supplementary Material), then applying the rest of the pulse up to T and finally measuring the expectation of H mol . Experimentally, this is feasible as long as the gate operation can be done quickly with respect to the drift of the system (e.g. τ V τ g , see Sec. V). Moreover, the implementation of the pulses will be exactly similar to the previous state preparation except for the one added gate, thus the experimental architecture remains largely unchanged.
In order to optimize for z ∈ Z ad , we employ the gradient descent algorithm with learning rate α k > 0 chosen by Armijo step rule [41]: where the z (0) l are initialized as low absolute valued constant functions on the interval [0, T ]. In order to implement the algorithm, the pulses are discretized as equidistant piecewise constant functions with N = T /τ ∈ N steps. For this discretization, the propagator can be expressed in closed form. Let t n := nτ , n = 0, . . . , N − 1, and set z l (t) := z l,n for t ∈ [t n , t n+1 ), l = 1, . . . , L. Then for t ∈ [t n , t n+1 ), the propagator U takes the form The minimization procedure is then executed by the algorithm below.
Algorithm 1: Discrete pulses VQOC algorithm The total number of quantum evaluations necessary in one parameter update is given by

V. RYDBERG PHYSICS
This section introduces basic Rydberg physics to identify what drift and control Hamiltonians as in Eq. (4) can look like on a Rydberg atom quantum computing system. This will also yield candidates for the control operators Q l , as in Eq. (5).
The following analysis focuses on a ground-Rydberg qubit system [11] (or Rydberg-Rydberg in the case of dipole-dipole interactions). Because the qubits never leave the H m manifold, only these states are taken into account (assuming that the system evolution time T is well below the lifetime of the states).
The Rydberg states have a passive 'always-on' interaction, which is described using a drift Hamiltonian H d , depending on the choice of qubits scheme [42], as a Van der Waals interaction [43] or a dipole interaction [44] where R ij is the interatomic distance. In this work, the qubit are arranged in a line with equal nearest neighbour distance, such that R ij = R|i − j|, with R being the nearest neighbour distance.
The coupling strength and the interaction strength are respectively characterized by the Rabi frequency Ω R and the interaction strength V = C i /R i , i ∈ {3, 6}. This leads to two characteristic timescales τ g = 1/Ω R , characterizing single qubit manipulation time, and τ V = 1/V , characterizing entanglement time.
To perform single qubit manipulations on qubit j, a nearly monochromatic laser interacts with the atom to realize the Hamiltonian [42,45,46] H 01 where 'h.c.' refers to the hermitian conjugate of the term that precedes it. Here, Ω j (t) [Hz] denotes the coupling strength, ϕ j (t) the phase of the laser coupled to atom j, and ∆ j (t) = ω j (t) − ω 0 [Hz] the detuning of the laser frequency ω j (t) from the energy level difference ω 0 . In the current state-of-the-art systems, the control over the laser coupling Ω j (t) and the detuning ∆ j (t) is much higher than that of the phase ϕ j (t) [47]. Therefore, often ϕ j (t) = ϕ j ∈ [0, 2π) is taken constant. It is however not unimaginable that laser phase control will improve significantly in the coming years.
The control Hamiltonian H c can contain several terms depending on which of the mentioned laser parameters can be controlled. In general, the control Hamiltonian will take the form (cf. Section IV) where z l (t) ∈ C [Hz], and Q l is a m-qubit operator. The choice of z l (t) being a complex number stems from the fact that it represents both the coupling strength and the phase in Eq. (8). However, as indicated before, the phase is often not controllable. In that case, z l (t) ∈ R is required. This would result in where Q l + Q † l is Hermitian. However, choosing Q l Hermitian from the start gives where z l + z l ∈ R. Thus, if there is no control of the phase of the complex number, the formulation of Eq. (9) remains valid by choosing Q l Hermitian.
In this part all the imaginable and realistic controls on a Rydberg system are listed, and at every future simulation it is indicated which terms are chosen to be included. respectively. Notice that the coupling Hamiltonian H coup c , already allows for full control on the Bloch sphere of each individual qubit. This is why this term is also referred to as rotational control. The entanglement control for respectively Van der Waals interactions and Dipole-Dipole interaction is given by (10) An entanglement control could for instance be realized by tuneable Rydberg dressing [45]. In many cases, the entanglement controls will be taken equal for all pairwise interactions, such that all z lk are equal for every combination of l and k.

VI. COMPARING VQOC AND VQE
This section outlines our devised method for comparing and contrasting VQOC and VQE. We compare both the expressibility of the algorithms and their performance in variational problems.

Equivalent evolution processes
In order to fairly compare VQOC and VQE, several algorithmic parameters need to be considered to create an equivalent evolution processes, meaning that the system evolves under the most similar circumstances for both algorithms. Concretely, this means that the algorithms are processed on the same hardware, with the same Hamiltonian, but with different control of the time-dependent interaction parameters to either do pulse-based or gatebased interactions.
First, the drift and control Hamiltonians in VQOC should be linked to the ansatz of VQE. In our simulations, both VQE and VQOC are considered for a grqubit [42] Rydberg system, where |0 is a ground state and |1 a Rydberg state. The passive interaction between |1 states is given by either the Van der Waals or dipole interaction Hamiltonian H d , as in Eq. (7) [11]. The qubits are considered on a straight line, equidistant from one another. For VQOC only rotational control is considered (L = m). The hardware-efficient ansatz is considered in VQE with alternating layers of single qubit rotations and m-qubit entanglement gates U ent as in Eq. (2). The assumption τ V τ g is made so that the hardware-efficient ansatz can physically be realized by taking U ent = exp(−iH d τ V ). Without this assumption, the system would passively evolve considerably while the single qubit manipulations are being executed, thus preventing them to be treated as independent. Each state on the Bloch sphere of a single qubit can be reached with a ZXZ-rotation. Such a rotation on a qubit q at a depth i can be written as where θ has three elements for every qubit and depth pair. For VQE, these ZXZ rotations are considered on the individual qubits in order to get analogous control to the rotational control for VQOC. We assume a Rabi frequency of Ω R = 1 kHz. This gives a gate time of τ g = 1 ms. Adhering to the condition τ V τ g can always be realistically done by varying R in V = C 6 /R 6 . Here we set V = 0.1 kHz. For the VQE algorithm, T VQE = d(τ g +τ V ) with d the depth of the hardware-efficient ansatz. To get similar evolutions, we take T VQOC = T VQE , which completes the construction of the equivalent evolution processes, as illustrated in Fig. 4.

Expressibility analysis
To compare the expressibility of VQE and VQOC, we follow an approach employed by Sim et al. [48]. In their work, the parameters of the evolution process are chosen according to some probability distribution. This results in a distribution over the unitaries. Expressibility of the evolution process is then equated to correspondence with the Haar measure, i.e. the uniform measure on unitaries [49]. To determine this correspondence, unitaries are drawn from the distributions and applied to an initial state |ψ 0 . The fidelity of the prepared state with respect to the initial state F = | ψ f |ψ 0 | 2 is determined, resulting in a distribution over fidelities. The distribution of fidelities of the evolution process is then compared to that of the Haar measure in terms of the Kullback-Leibler (KL) divergence D KL , which is a measure of the similarity between two distributions P and Q, given by where X is a probability space. The distribution of fidelities for Haar random states is given by where N is the dimensionality of the unitaries [50]. We consider a VQE hardware efficient ansatz as described above and the VQOC equivalent evolution process. Initial gate parameters/pulses are then chosen randomly for both VQE and VQOC. In VQE, the rotation angles are chosen uniformly in the interval [0, 2π]. For VQOC the piecewise constant values are chosen uniformly at random between [−1, 1] kHz. In the results, it will be shown that even for a low choice of N , thus limiting the degrees of freedom in VQOC, the expressibility, in terms of correspondence to the Haar measure, is higher in VQOC than for the equivalent VQE circuit.

Variational problems
Next, the algorithms are quantitively compared for their performance in solving variational problems. A good quantity used to compare two quantum optimization algorithms is the total number of circuits run on the quantum computer, during the algorithm. This holds especially true in the NISQ era, where the availability of quantum computing power is limited. Based on Eqs. (3) and (6), we find the ratio of number of quantum evaluations for both methods per iteration to be Recall that K is the number of decomposition terms in the control factors, L is the number of controls scaling linearly with m, and finally, N = T /τ , the number of considered time steps in VQOC. For all simulations in this work, we set K = 2, L = m for rotational control, and N = 100. In order to characterize the performance on variational problems, several small molecule ground state energies are approximated using both algorithms for a fixed number of total quantum evaluations.
Obviously, many other ansatzes for VQE, as well as discretizations for VQOC, are possible. The hardwareefficient ansatz and piecewise constant discretization are the most common in the literature, and we foresee no reason why other procedures should yield substantially different results. Nonetheless, different implementations and comparisons might yield interesting results in future research. Figure 4: The equivalent evolution processes used to compare VQE and VQOC in the ground-state-energy problem. VQE is conducted on a hardware-efficient ansatz with d layers of alternating ZXZ-gates, taking time τ g , and passive interaction via the drift Hamiltonian, for a time τ V . VQOC is executed with piecewise constant functions consisting of N steps and L = m rotational controls (giving K = 2). We assume Ω R = 1 kHz, leading to τ g ≈ 1 ms. Adhering to τ V τ g , we choose τ V = 10 ms which in turn gives V = 0.1 kHz and T VQE = T VQOC = 11d ms.

Expressibility analysis
The expressibility analysis, as described in Sec. IV, is performed between VQE and VQOC. Because of finite sampling size, the obtained fidelities are binned together in equally sized bins in the interval [0, 1]. The KLdivergences are calculated over the obtained histograms. Fig. 5 depicts the results for d = 2, V = 0.1 kHz, N = 4 and T VQE = T VQOC = 22 ms. Fig. 5a shows that for 10 4 samples, the VQOC evolution procedure approximates the Haar measure considerably better than VQE with d = 2, as it is less concentrated around fidelity F = 0. Interestingly, this is already significant for the low number of time steps N = 4. Furthermore, Fig. 5b shows the Bloch sphere representation of the partially traced first qubit density matrices ρ 1 = Tr 2 [|ψ f ψ f |]. From these plots, one observes that the density matrices for VQE concentrate more on the poles of the Bloch sphere than for VQOC, indicating less expressibility in terms of correspondence to the Haar measure.

Electronic structure problems
This work considers ground state problems for H 2 , LiH and H 4 molecules. The H 2 and LiH molecules are aligned on the x-axis with varying interatomic distance. For H 2 the 1s orbitals are considered active, higher energy orbitals are ignored [51]. This results, after spin and electron number reduction [27], in a 2-qubit problem. For LiH, the 1s orbital on the Li atom is fixed as fully occupied. The 2s and 2p x orbitals of Li and the 1s orbitals of H are considered active. This results in a 4-qubit problem after spin and electron number reduction. For H 4 , the atoms are aligned on the x-axis with all interatomic distances between neighbouring atoms taken equally. Only the 1s orbitals of all 4 H atoms are considered active. This results, after spin and electron number reduction, in a 6-qubit problem [52].
The choice for these molecules as illustrating use cases is based on the fact that all possess weakly and strongly entangled ground states at the interatomic distances considered. For all molecules considered, the exact H mol ground state energy has been calculated using the FCI method [51]. In certain simulations, random Hamiltonians are considered. These Hamiltonians are constructed by taking a fixed number of Pauli strings at random and picking weights in the interval [−1, 1] uniformly at random. The reason for using these random Hamiltonians is that, by construction, they will possess (with high probability) a highly entangled ground state. This state will also (with high probability) be a linear combination of all basis states, and therefore will have a high QSL w.r.t. all basis states.
We note that the systems considered in this work are small in terms of the number of qubits. In larger systems, for both gate-based and pulse-based control, a trade-off will have to be made between expressibility and optimisability due to the presence of barren plateaus [53].

Variational problems
As a first qualitative illustration of VQOC, the ground states of 2-qubit H 2 Hamiltonians (4-qubit LiH Hamiltonians) are calculated for varying interatomic distances, with the vacuum state (|ψ 0 = |00...0 ) as the initial state. VQOC is run for 100 iterations with rotational control H coup c , L = m, K = 2, Rydberg interactions as in Eq. (7) with C 6 /R 6 = 0.1 kHz, T = 100 ms and N = 100 [54]. The results in Fig. 6 for H 2 show that VQOC is able to find energies below the Hartree-Fock energy and reach chemical accuracy [18] after as few as 30 iterations. The errors reached are comparable to other (pulse-based) VQAs [2,55,56] and well below chemical accuracy. For LiH, all possible basis states (|00...00 , |00...01 , ...|11...11 ) are taken as initial states, after which the mean and best case energy results are reported, as similarly done in Ref. [57] by Rabitz et al. Fig. 7 shows that VQOC is able to get below the Hartree-Fock state energy and reach chemical accuracy [18] for LiH after 50 iterations in many occasions.
As another qualitative analysis, Fig. 8 shows final pulses and density matrix evolution for simulations with rotational and entanglement control as in Eq. (10) (L = m + 1, K = 2), Rydberg interactions with V = 0.1 kHz, T = 100 ms and N = 100. It is shown that the density matrices are closely approximated, and the achieved energy errors, which are around 10 −5 Hartree, are well below chemical accuracy. Moreover, the fidelities F of the prepared state with respect to the ground state are close to 1, as desired.   To compare the VQE and VQOC algorithms as described in Sec. V, LiH and H 4 Hamiltonians are analyzed. Figures 9 and 10, respectively, show VQE vs. VQOC results for LiH with d = 9 and H 4 with d = 4, giving T VQE = T VQOC = 99 ms and T VQE = T VQOC = 44 ms respectively. In both cases, all possible basis states are taken as initial states [18] and the mean and best cases are reported.
The first thing to note is that for both types of entanglement, the best VQOC result generally outperforms the VQE results, while the average results are quite similar, even slightly better for VQE in the H 4 case. Moreover, the VQOC results are way below chemical accuracy, whereas this is not reached for VQE. For LiH in Fig. 9, the dipole interaction results perform well in the intermediate interatomic distances (where the ground state is highly entangled), while the Van der Waals interaction results perform better in the low and high interatomic distances (where the ground state is lowly entangled). For H 4 , similar behaviour is seen in Fig. 10, as VQOC performs better at low interatomic distances, where the ground state is lowly entangled, than at high interatomic stances where the ground state is highly entangled. The reason for this might have to do with the qubit architecture, e.g. for the LiH Hamiltonians, two non-neighbouring qubits (1 and 3) have to be entangled without being entangled to qubit 2, which is challenging given the strong nearest neighbour interactions in a Rydberg system. This indicates the importance of architecture choice in VQOC and will be a future point of investigation, as previously done in VQE [27].   9b shows two convergence behaviour plots. The left plot shows illustrative behaviour for the low interatomic regime where VQOC with dipole interaction performs best and the behaviour of the other methods is quite comparable. The right plot shows illustrative behaviour for the high interatomic regime, where VQOC is generally better than VQE for both interaction types, and the Van der Waals interaction excels. Fig. 10b shows convergence plots for H 4 , where it is seen that VQOC is able to find energy levels below chemical accuracy, whereas VQE stifles.

QSL Results
To highlight the advantage of VQOC opening up a bigger search space at low evolution times T , we consider randomly generated 4-qubit Hamiltonians, generated as described at the beginning of this section. For VQE we again consider the hardware efficient ansatz with varying d, τ g = 1 ms, Van der Waals entanglement gates with τ V = 10 ms and entanglement time T ent ∈ [0, 2πτ V ] resulting in the entanglement gate exp(iH d T ent ). This results in total evolution times T VQE = d · (τ g + T ent ). We then consider the equivalent VQOC circuit, as described in Sec. VI, with T VQOC = T VQE . Note that at time T ent = 2πkτ V , k ∈ N, the Van der Waals entanglement gate is equal to the identity. In VQOC, something similar happens when T VQOC /N = 2πkτ V , k ∈ N as the step size becomes too small to capture the behaviour caused by the drift Hamiltonian. Fig. 11, shows an example of the errors w.r.t. the ground state of a randomly generated Hamiltonian after 64 · 10 4 total quantum evaluations. Here it can be seen that the reached errors and QSL are independent of circuit depth (for d ≥ 3) for VQE and that VQOC reaches lower errors faster than VQE can. This supports the hypothesis that VQOC can open up a larger search space of states than VQE can, given a fixed evolution time. The errors increase again because of the entanglement gates transforming back to identity occurring at the expected times. Figure 11: Errors achieved on a 4-qubit random Hamiltonian for VQOC and VQE at different values of T VQE = T VQOC after 64 · 10 4 total quantum evaluations.

VIII. CONCLUSION
In this work, we introduce a novel pulse-based variational quantum optimal control (VQOC) algorithm for quantum chemistry applications. The method has, through simulations, shown proof of concept for a neutral atom quantum computing system with interactions between the qubits mediated via Rydberg excitations, by approximating ground state energies of simple molecules with chemical accuracy. In this process, hyperparameters like the evolution time T and the terms in the control Hamiltonian H c were found to play important roles.
Comparing VQE and VQOC is not trivial, as both algorithms use entirely different optimization methods. For our criteria of quantum overhead, we have shown that our VQOC method can compete with, and in certain cases outperform VQE, while doing so in a more NISQfriendly manner. Especially when considering equivalent evolution process at small depths (short T ), VQOC can outperform VQE. In particular, the better convergence of VQOC in the low evolution time regime will help stifle the influence of decoherence in NISQ-era quantum computing problems. The reason for the faster convergence lies in the accessibility of the Hilbert space, and thus links to the QSL. In regimes with higher T , depending on the type of drift Hamiltonian introduced and the ground state of the problem Hamiltonian, VQOC either outper-forms VQE or has similar convergence behaviour.
The exact reason for this behaviour should be part of future research. However, we believe that for systems where the interatomic interactions are very strong, giving rise to a strongly-entangled ground state, VQOC is able to find this state faster due to entanglement being accounted for at every time step. Furthermore, there is a difference in convergence behaviour under varying entanglement schemes, which could be a result of a nonoptimal qubit arrangement. An obvious solution to this would be to try and find optimal permutations of the qubits or fit the entanglement scheme to the problem Hamiltonian. A perhaps more interesting solution would be to treat the qubit positions as optimizable parameters (which due to technical constraints would be timeindependent) in the VQOC scheme.
Other points of future research would include the physical implementation of the algorithm. Specifically, the Rydberg version of the VQOC algorithm would be a perfect implementation on a neutral atom optical tweezer platform with Rydberg excitations. Several points to focus on in an experimental implementation are the repeatability of the algorithm, the robustness due to errors within the pulses and the required precision of the pulses, and creating specific controlled gates as required by the parameter gradient algorithm. Another implementation issue is the matter of noise. On the pulse-based level, certain noise terms can be introduced eloquently by switching to Lindbladian evolution instead of Hamiltonian [58]. A last interesting point would be to look into adapting the VQOC algorithm for qutrit or qudit manifolds [59] and its influence on barren plateaus [60,61].