{'fraction_non_alphanumeric': 0.028233440874873148, 'fraction_numerical': 0.007527907786916284, 'mean_word_length': 4.424979457682826, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Counterdiabatic Optimised Local Driving

Adiabatic protocols are employed across a variety of quantum technologies, from implementing state preparation and individual operations that are building blocks of larger devices, to higher-level protocols in quantum annealing and adiabatic quantum computation. The problem of speeding up these processes has garnered a large amount of interest, resulting in a menagerie of approaches, most notably quantum optimal control and shortcuts to adiabaticity. The two approaches are complementary: optimal control manipulates control fields to steer the dynamics in the minimum allowed time while shortcuts to adiabaticity aim to retain the adiabatic condition upon speed-up. We outline a new method which combines the two methodologies and takes advantage of the strengths of each. The new technique improves upon approximate local counterdiabatic driving with the addition of time-dependent control fields. We refer to this new method as counterdiabatic optimised local driving (COLD) and we show that it can result in a substantial improvement when applied to annealing protocols, state preparation schemes, entanglement generation and population transfer on a lattice. We also demonstrate a new approach to the optimisation of control fields which does not require access to the wavefunction or the computation of system dynamics. COLD can be enhanced with existing advanced optimal control methods and we explore this using the chopped randomised basis method and gradient ascent pulse engineering.


I. INTRODUCTION
Time-dependent manipulation of few and manyparticle quantum systems is important across all implementations of quantum computing and simulation. In such processes, decoherence and undesired transitions reducing the state fidelity are relatively ubiquitous. One important example is given by the undesired transitions that can occur between instantaneous eigenstates of the dynamical Hamiltonian upon the application of an external drive. This is why many driving protocols rely on adiabatic dynamics, where the system follows the instantaneous eigenstates and transitions are naturally suppressed. Ideal adiabatic processes are reversible making them -in principle -robust. However, to approach ideal adiabatic processes the dynamics must always be very slow, requiring compromises on the time-scales of competing heating and decoherence processes.
Speeding up adiabatic protocols to enable their completion within the system's coherence time is important for the development of any quantum technologies relying on such protocols [1]. One approach to do this is the implementation of optimal driving protocols, which aim to end up with the system in a desired final state. For example, numerically optimised paths can be employed to avoid points where gaps in the spectrum of the system become small, or additional control fields can be tuned to increase the size of these gaps [2][3][4]. In broad terms, this is the goal of protocols collectively referred to as quantum optimal control. Another option is to design techniques which speed up the adiabatic dynamics, often termed * ieva.cepaite@strath.ac.uk shortcuts to adiabaticity (STA). The primary aim of STA is to entirely remove or suppress diabatic transitions between instantaneous eigenstates of the dynamical Hamiltonian [5,6]. One particularly successful technique is counterdiabatic driving (CD), which was first utilised in physical chemistry by Demirplak and Rice [7,8], and was independently introduced by Berry [9] under the name 'transitionless driving'. CD aims to suppress losses that arise due to fast deformations of the system far from the adiabatic limit by analytically compensating for them in the Hamiltonian. In general, to suppress diabatic losses exactly, the full analytical or numerical solutions of the Schrödinger equation are required. This makes the implementation of CD in complex systems -e.g. for manybody dynamics -difficult and requires the need for new techniques to be introduced.
Links between optimal control and STA have existed throughout the development of both approaches [10][11][12], but there are few examples of their explicit combination in a way that exploits their complementary nature. Some attempts to achieve this have included an emulation of CD through fast oscillations of the original Hamiltonian [13,14] as well as through recent advances in reinforcement learning methods aimed at optimizing quantum protocols [15]. Such methods have been shown to achieve a significant improvement in performance when implemented using concepts borrowed from CD [16]. In this work, we offer a significantly different new approach in combining elements from STA and quantum optimal control which we will call counterdiabatic optimised local driving (COLD).
A key ingredient in the development of COLD is a recent approach designed for implementing CD in the setting of larger, more complex systems: local counterdiabatic driving (LCD) [17][18][19]. LCD offers a method to arXiv:2203.01948v2 [quant-ph] 24 Oct 2022 derive approximate CD protocols, with the aim of suppressing undesired transitions instead of fully eliminating them. This allows it to account for some physical constraints of the system, e.g. locality conditions. However, the approximate nature of the LCD protocol can lead to poor performance, necessitating the introduction of additional non-local, long-range corrections [17,20,21]. If all possible corrections are added, then LCD is equivalent to the normal analytical approaches of CD, but the additional terms are generally difficult to control experimentally. COLD offers an alternative approach, with additional control fields which allow for an optimisation of the dynamical Hamiltonian for a given local form of LCD. The impact of more complex corrections can then be radically reduced, giving a corresponding improvement in the desired protocol.
An important consequence of optimising for a given local order of the LCD is the possibility of bypassing the need to have access to the wave function, dynamics or experimental data of the given system in order to perform the optimisation. LCD is an analytic method and can be calculated using only the coefficients of the Hamiltonian. We find that it is possible to perform numerical optimisation of the path of the system by simply minimising higher order LCD integrals and/or amplitudes, a method that not only bypasses the need to compute system dynamics but is also independent of system size. This makes it an exceptionally useful tool in practical settings.
The structure of this paper is as follows: first, we give a detailed description of the new method, COLD, with a focus on the elements of quantum optimal control and CD required for its implementation. In Sec. III we explore a 2-spin annealing protocol, that showcases the strengths of COLD. Sec. IV describes and analyses the improvements gained with COLD and its combination with other optimal control techniques in the case of state preparation in the Ising model as well as the potential computational advantage of optimising higher order integrals of LCD instead of the final state fidelity. In Sec. V we show the improvement that COLD can achieve on the recently realised example of LCD for state transfer on a synthetic lattice in ultracold atoms and in Sec. VI we demonstrate that, when implemented for second order LCD, COLD can be used to quickly and effectively prepare highly entangled multipartite GHZ states. Finally, in Sec. VII we explore the possibility of minimising both the power and amplitude of higher order LCD drives as a means to efficiently optimise COLD parameters, bypassing the requirement of computing system dynamics. A list of abbreviations used in this work can be found in Table. I for reference. In the context we consider, we employ quantum optimal control to optimise the function f (ψ, β) in the Schrödinger equationψ where ψ is the quantum wave function and β is the set of optimisable control parameters. Optimisation of Eq. (1) in most cases means taking the system from an initial state |ψ 0 to a final target state |ψ T by finding the optimal values of β with respect to some target metric (e.g. the time taken to evolve the system from |ψ 0 to |ψ T ). There is a large variety of techniques available to achieve this goal [3,22]. The success/target metric needs to be defined prior to the optimisation of β. Often this is done by constructing a cost function, which in turn defines the optimisation landscape. In general, we can optimise for any desired property of the final state of the system, with some examples being the entropy, energy, energy fluctuations or some other observable. A commonly used cost function in state preparation is related to the fidelity of the final, post-evolution state |ψ f (β) with respect to the target state: In performing such a numerical optimisation, it is common to take the target state to be parameterised via a Hamiltonian split into two parts. The first is the so-called bare Hamiltonian H 0 (t), which can be time-dependent and describes the dynamics of the quantum system in question. The second part is then an additional driving term that includes a function f parameterised by the control parameters β, as well as operators O opt which provide additional degrees of freedom in the dynamics of the system. The full Hamiltonian of the control system is then: The parameters β can then be optimised for the optimal dynamics with respect to the metric defined by the cost function.
In this work, we generally use the Powell minimization [23] and dual annealing [24] approaches for the numerical optimisation as implemented in Python's SciPy [25]. When performing this optimisation without any CD terms in the Hamiltonian, we refer to them as bare Powell optimisation (BPO) and bare dual annealing (BDA) respectively, with bare referring to the lack of CD. Furthermore, we implement the chopped randomised basis (CRAB) approach [26,27] and combine its methodology with that of COLD, to obtain the method of COLD-CRAB. CRAB expands the size of the parameter landscape by employing randomisation, usually in the optimised pulse driving the system. The approach was first developed for quantum many-body problems whose simulation requires the time-dependent density matrix renormalization group, despite the fact that these were thought to not be tractable in the quantum control setting [27,28]. CRAB has benefits in that it can avoid traps in the control landscape [29], and has built-in flexibility for open-loop or closed-loop optimisation [27,30] although these advantages come at a higher computational cost due to requiring a far larger search-space for the optimisation.

B. Counterdiabatic Driving
An important class of optimisation problems deals with the case where the initial and final states are ground states of a Hamiltonian H 0 (t) at some initial t = t i and final t = t f time. In these cases, the adiabatic theorem guarantees that for an infinitesimally slow transformation of the system t f − t i → ∞, it should follow the instantaneous (non-degenerate) ground state of H 0 (t) and hence reach the target state with unit fidelity. This process is generally known as quantum annealing.
In large, complex systems with many degrees of freedom, quantum annealing tends to require prohibitively long protocol times due to vanishingly small gaps typically present in such systems. This often makes annealing protocols impractical [31,32]. It has been found that this problem can be formally overcome by using CD, where velocity-dependent terms are added to the Hamiltonian analytically enforcing the adiabatic wave function to be the solution of the time-dependent Schrödinger equation [7][8][9]. In this case, the dynamical state will follow the instantaneous eigenstate with no transitions regardless of the driving time. The form of the dynamical Hamiltonian enforcing this is [9]: with |n ≡ |n(t) the n-th eigenstate of the instantaneous Hamiltonian H 0 (t). The last term enforces the phases ( n|∂ t n ) on the instantaneous eigenstates, which are arbitrary and thus will be omitted. In general, knowledge of the CD Hamiltonian of Eq. (4) requires knowledge of the full spectrum of H 0 (t) at each instantaneous moment in time.

C. Counterdiabatic Optimised Local Driving
We will now introduce the main idea of COLD. The principle is to take the same approach as Sec. II B but with the original Hamiltonian given by H β (t, β), see Eq. (3). Quantum Annealing then applies to the whole family of control Hamiltonians H β (t, β) as long as the additional control function f (t, β) vanishes at the protocol boundaries: f (t i , β) = f (t f , β) = 0. This flexibility was explored in finding the optimal adiabatic path characterized by e.g. the shortest distance between the initial and the final states, i.e. a geodesic [33]. A similar geodesic approach was developed in the context of dissipative systems to minimize energy losses [34]. During the protocol, a dynamical Hamiltonian H β (t, β) generally induces transitions between the quantum states that it drives and the question about what is the optimal path remains open.
The Hamiltonian H β (t, β) and its eigenstates depend on time only through the driving parameters, which include β and any additional control terms in the particular protocol. This makes it convenient to introduce a path in the coupling space parametrized by a dimensionless parameter λ ∈ [0, 1] such that both H 0 and f are functions of λ satisfying H β (λ = 0) = H 0 (t i ) and H β (λ = 1) = H 0 (t f ), i.e. being equal to the initial and the final Hamiltonian at the protocol boundaries. By construction this implies that any additional fields introduced to the bare Hamiltonian H 0 must go to zero at the boundaries. In this way, any protocol can be uniquely characterized by first specifying the path f (λ, β) in the coupling space manifold and then the time dependence λ(t) along it. The path determines the sequence of couplings of the Hamiltonian during time evolution and hence the sequence of ground state wave functions followed by the driven state. Furthermore, the time dependence encodes the speed of traversing this path. We can then introduce a hermitian operator called the (path-dependent) adiabatic gauge potential [17]: A λ = i n |∂ λ n n|, which satisfies a closed form equation, where: with both H β and A having a dependence on λ and β(λ). We note that in the case of a nonlinear Schrödinger equation where the dynamics are described by classical Hamiltonian equations of motion, the commutators need only be replaced with Poisson brackets and the same idea applies [19]. Thus, the CD Hamiltonian reads and is equivalent to the CD Hamiltonian of Eq. (4) given knowledge of the exact adiabatic gauge potential. However, generally the adiabatic gauge potential is a very non-local object and solutions of Eq. (5) are unstable to small perturbations containing exponentially many terms in the number of degrees of freedom. LCD aims to find approximate gauge potentials that satisfy particular requirements like robustness and locality, thus circumventing many of the difficulties in determining the second component in Eq. (4) and (7) exactly. The goal, in essence, is to suppress the most relevant diabatic effects rather than completely eliminate them. This method has recently been experimentally implemented to speed up state transfer for synthetic lattices in ultracold atoms [35], for preparing states in nuclear-magneticresonance systems [36], and annealing protocols on an IBM quantum computer [32,37].
Following the methods of Ref. [17], the problem of finding the optimal adiabatic gauge potential can be cast as the minimisation of the Hilbert-Schmidt norm of G λ , which is equivalent to minimisation of the action with respect to A λ . In most cases, this is achieved by first choosing an operator ansatz -i.e. a set of linearly independent operators {O LCD } -and then using this set as an operator basis for the adiabatic gauge potential LCD . The action can then be minimized with respect to the the set of coefficients, α. The choice of operators {O LCD } can be made easier when noting that A λ acts as a generators of motion in the parameter space. This implies that for, say, real Hamiltonians like those we'll be exploring in the following sections, a good Ansatz for the adiabatic gauge potential is one which is non-interacting and imaginary. In the example of an Ising spin chain we may take A λ = N j α j σ y j , where j labels the N chain sites, and {O LCD } is a set the y-pauli matrices.
Without any additional control fields f (λ, β), LCD is essentially an informed choice of the operator set {O LCD } in a way that the resulting control protocol from the minimisation of Eq. (8) is optimal for a given H 0 (λ). In this case we explore the family of Hamiltonians The performance of such LCD protocols is determined by how accurately the variational manifold spanned by the set {O LCD } can approximate an exact A λ such that Eq. (5) holds.
In the case of the new protocol COLD, we allow for extra exploration of the family of Hamiltonians due to the additional control fields as in Eq. (3). This expands the family of Hamiltonians to Note, that the coefficients of the optimal control field change the form of the LCD driving coefficients, i.e. α → α(λ, β). The aim of COLD is then to optimise the coefficients β in such a way that the LCD term in the above equation allows for the greatest suppression of non-adiabatic effects for the dynamical Hamiltonian H 0 (λ) + f (λ, β)O opt . One can picture it as changing the path that the system takes between its initial and final states, with the express goal of picking a path that maximises the effects of the approximate counterdiabatic drive given by the second term in the equation. This path will depend on the form of the optimal pulse function, the operators O opt and on the values of β. We will focus on the optimisation of the control parameters β for a given choice of f (λ, β) and O opt , although the choice of operators O opt as well as the function of the control pulse f (λ, β) can also be optimised over as an extension.
With COLD, we have two methods to improve on the existing LCD protocol. As previously shown in Refs. [38,39], there is a possibility to add more terms to the LCD making it less local, e.g. through long-range interactions. In the spin chain case, we could take the aforementioned sum over σ y terms to be the first-order anzatz for the LCD, where higher-order ansätze might contain sets of operators {O LCD } with terms odd in σ y such as σ y j σ (z,x) j+1 . This procedure generally improves the performance of CD protocols at the expense of adding more complex operators which may be experimentally impractical depending on the scenario. Alternately, with COLD and the introduction of additional local control fields to the Hamiltonian, we can improve the performance of LCD at a fixed complexity of the CD term by significantly modifying the adiabatic landscape at intermediate couplings to enhance the performance of the given order of LCD.
In this work we pursue two directions of optimising the local control fields: numerical optimisation of the dynamics and minimisation of the higher order LCD terms. For the most part we will focus on numerical optimisation of the dynamics directly, as these will reach optimal values for specific protocols when implemented. COLD opens the possibility of minimising the higher order LCD terms instead, which benefits from not requiring calculation of the systems dynamics. This approach, as discussed in Sec. VII, allows for optimal control procedures using COLD to be implemented for large systems that would be cumbersome for procedures that require the numerical optimisation of the dynamics.
We also note that while it may seem prudent to treat the LCD coefficients α(t) as control pulses which may be parameterised and optimised in the same way that f is, we find that this method fails to perform well compared to using the variational form of the LCD as we have done. This is likely due to a difficulty in choosing an accurate form of the drive as well as parameterising it. On top of that, the loss function space in this case may become intractable when compared to that of COLD as we have presented it in this section.
Furthermore, we compare COLD to the use of CRAB, as discussed in Sec. II A. An advantage of COLD is that it can be combined with many advanced optimal control procedures, owing to the standard way additional control fields are introduced to the Hamiltonian. In this work we find the combination of COLD and CRAB particularly useful and we will refer to this as COLD-CRAB.

III. TWO SPIN QUANTUM ANNEALING
To showcase and explore the use of COLD in a relatively simple setting we will consider a two spin quantum annealing problem with bare Hamiltonian j , a ∈ {x, y, z} are the Pauli matrices applied to spins indexed by j. For the scaling function λ(t) we pick the form such that λ(0) = 0 and λ(τ ) = 1. We consider the case of J/h = 0.5, which means the spins start in the initial state of |↑↑ and finish in a superposition of all of the symmetric states.
As discussed in Ref. [17], since H 0 has a standard Ising spin chain form, the first-order LCD terms are given by the following ansatz for the adiabatic gauge potential: with the sum being over the full length of the N spin chain. Minimising Eq. (8) for this A λ with respect to the coefficient α gives To further improve on the first-order LCD we can implement COLD, as we will discuss shortly, or we can introduce higher-order terms to the ansatz for A λ . This second method serves as a good benchmark against COLD, since it offers an improvement to first-order LCD in the same way as COLD does, but requires more complicated interactions between the two spins increasing the implementation overhead. The second-order LCD can be found by taking an ansatz for the adiabatic gauge potential: where to solve for α, γ, and ζ we once again minimize the action given by Eq. (8) and obtain three coupled equations which can be solved numerically (see Appendix A for a detailed derivation). We now consider three distinct cases in this two spin quantum annealing example: no LCD, first-order LCD, and second-order LCD. The fidelity of the final state for each case over a wide range of driving times is shown in Fig. 1(a), with an easily distinguishable advantage in the case of LCD. The final fidelity where no LCD is implemented decreases rapidly as the ramp times are made short, with the system getting stuck in its initial state. On the contrary, first-order LCD retains good final state fidelities into short times, as the driving Hamiltonian becomes that of only the LCD term. The second-order LCD then gives unit fidelity, in agreement with previous observations [38], as for a two spin Hamiltonian the highest order corrections are that including two spin terms.
We now add an optimisable term, as described in Sec. II A, so that the new Hamiltonian reads: with N k the number of optimisation coefficients β, and β k ∈ β the coefficient of the kth frequency of the control function. Note that we consider with The form of the additional control field fulfils the requirement that the boundary conditions are H(0) = H 0 (0) and H(τ ) = H 0 (τ ). Note that Numerically optimising the β k for the best final state fidelity without adding LCD terms results in the BPO method introduced in Sec. II A. We show the results of BPO in Fig. 1(b), where it is observed that BPO gives better results than the case of no LCD in Fig. 1(a). However, for short times the BPO approach still results in the system getting stuck in the initial state.
Finally we present and compare the results of the new method, COLD. In this case the Hamiltonian before adding LCD terms is given by Eq. (16) and the coefficient of the first-order LCD is Note that the optimisation of the additional control field also feeds into the coefficient of the adiabatic gauge potential during the dynamics as discussed in Sec. II C. The results of the COLD approach for this two spin annealing protocol are shown in Fig. 1(b), where we observe an improvement of the final state fidelity beyond what is possible with first-order LCD alone in Fig. 1(a). In this example, LCD alone reaches a final state fidelity of 1 − F = 3% at short times, however COLD improves this error in the final state to 1 − F = 0.005%. This is due to the extended family of dynamical Hamiltonians in Eq. (10) owing to the addition of an optimisable control field. This result shows that COLD can provide an advantageous alternative to the addition of higher-order LCD which may be experimentally impractical.
We have found that COLD performs better than LCD of the same order or BPO when the system dynamics are calculated numerically. This does not, however, imply anything about the performance of COLD in more complex scenarios, like in the case of an unknown target ground state. In that case the fidelity is a poor optimisation metric. There is, however, a way to come to the same conclusions as those presented in Fig. 1 without the need to compute the dynamics exactly. We can do this by first using a guess for the COLD protocol to find the approximate adiabatic gauge potential and then minimising the integral of the norm of the second-order correction to the adiabatic gauge potential along the path. Note, that the ground state can be in turn obtained through first order COLD, so there is no need to diagonalize the Hamiltonian. This integral should be small if COLD has implemented a dynamical Hamiltonian that makes the first-order adiabatic gauge potential the leading term. It is effectively a measure of the error of COLD and can be given by with |ψ g (t) the instantaneous ground state along the path and one of the second-order correction terms. In order to confirm this is the case, we compare the different paths -COLD and LCD only -in the two-spin example in order to determine if I 1 is small for COLD. If I 1 is small when compared to the same measure for lower-order LCD as t → 0, then we know that COLD is enforcing a better dynamical Hamiltonian. In the case of the two spin annealing protocol we find that as t → 0, I 1 → 0.04 for COLD and I 1 → 0.2 for LCD, showing that COLD is minimising the second-order correction along the path. A simpler integral also reflects this correction in this two spin example, with I 2 → 0.03 for COLD and I 2 → 0.1 for LCD as t → 0. This is particularly useful in more complex scenarios as I 2 is relatively straight-forward to calculate, as we will demonstrate in the next section. We also observe the reduction of the corresponding integrals of the (σ y 1 σ z 2 + σ z 1 σ y 2 ) term of the second-order LCD. By minimising these integrals, it is possible to extend the COLD approach to more complex scenarios, including where the exact calculation of the dynamics is not possible.

IV. 1D ISING MODEL
In this section we apply COLD for state preparation on a 1D Ising spin chain in the presence of a transverse and longitudinal field. We consider an annealing protocol where the aim is to prepare the ground state across the Ising phase transition. The annealing Hamiltonian is given by where Z 0 is a small offset parameter to break ground state degeneracies and X f is the final x-field strength.
Note, the breaking of the ground state degeneracies is not a requirement but allows for easier consideration of the adiabatic path. As before, λ(t) is a scaling function that has the boundary conditions λ(0) = 0 and λ(τ ) = 1, with τ the driving time. This means we start from the ground state of all spins up and drive across the quantum phase transition to the ground state which is a superposition of all basis states. We again take the scaling function to be given by Eq. (12). In this example, we use X f = 10J and Z 0 = 0.02J. For the Hamiltonian of Eq. (23), the LCD to first and second order is well known, as the wave functions are entirely real. We take the first-order adiabatic gauge potential to be given by where the coefficients for the general periodic spin chain of Eq. (23) are Note that the quoted α above is technically for a periodic or infinite size system, with J 2 → J 2 (1 − 1/N ) for a finite system. However, we find that the inclusion of the factor for the finite system sizes we consider only changes the final achieved converged fidelities at short times by ∼ 10 −6 %. The second-order adiabatic gauge potential is of the form with the coefficients α, γ and ζ again obtained by minimising the action given by Eq. (8) and solving the coupled set of equations numerically (see Appendix A for a detailed derivation).
In this example, optimal control is implemented by introducing an additional driving field so that the dynamical Hamiltonian is given by Figure 3. Maximum amplitudes of CD terms in the Ising model annealing protocol for (a) first-and second-order LCD only with no additional optimal control fields and (b) the COLD approach optimised for the best final state fidelity implementing first-order LCD as shown in Fig. 2 (a). The plot shows the maximum amplitude at each driving time for the first-order α (red circles) and the two second-order terms γ (blue diamonds) and ζ (green triangles) as given in Eq. (26) (although the second-order LCD is not actually implemented in COLD). An inversion in the strength of the second-order and first-order LCD terms for (a) no additional optimal control fields and (b) the addition of optimal control fields shows that COLD implements a dynamical Hamiltonain which is favourable for the applied order of LCD (first-order in this case).
with β being the terms to optimise over. We take our additional terms to again respect the boundary conditions f (t = 0, β) = 0 and f (t = τ, β) = 0, meaning a natural choice is with ω k = 2πk the kth principal frequency and g(λ) given by Eq. (18). To implement the CRAB algorithm discussed in Sec. II A, we will use k → k(1 + r k ) instead with r k drawn from a uniform random distribution r k ∈ [−0.5, 0.5]. Note that there is a strong distinction between the CRAB, which is an established optimal control method in its own right and our own version COLD-CRAB, which includes an LCD term along with the optimal control field in the Hamiltonian. To be more precise, the COLD-CRAB Hamiltonian can be expressed as: where for each optimisable parameter β k associated with a k t h principal frequency we also assign a random value r k ∈ r as described earlier. Note that the dependence on r is inherited by the LCD drive term α, since it is a function of f (t, β, r).
As before, we choose the first order adiabatic gauge potential given by Eq. (24) and find that the coefficients are Note, with the introduction of the additional control fields f it is possible for α to be non-zero at the start or end of the protocol, asḟ is not fixed to be zero. However, this can be enforced by a suitable choice of the additional control field, we will consider replacing α → S(λ)α where S(λ) is a scaling function that tends to zero as λ → 0 and λ → 1. We find that the scaling function only has a minimal effect on the final fidelities observed. This issue could also be resolved by a suitable choice of f , with our example drive being an extreme case asḟ is maximal at the boundaries of the protocol. Note that this issue is present in LCD as much as in COLD and we have chosen to highlight it here as it may become a concern in an experimental setting where a discontinuous drive is simply impossible at the beginning and end of a protocol. The suitable choice of the form of f in a given example is a problem we will leave for future work, with our focus being on the introduction of the COLD protocol. We first compare the final state fidelity when using COLD versus BPO as shown in Fig. 2(a) for different driving times in a system of N = 5 spins and a single N k = 1 optimisation coefficient. As expected, at long timescales the two methods agree as we approach the adiabatic limit of the dynamics. However, at shorter time scales the difference in behaviour is dramatic. We observe that the BPO approach fails in the case of very fast driving as the state gets stuck in the initial state but the COLD approach converges to 1 − F ∼ 10 −3 . We note that the advantage achieved by COLD is not due to the introduction of first-order LCD terms alone, as in Fig. 2(a) we see that this will result in F = 0.0440 for τ = 10 −3 J −1 . COLD is instead achieving this by making the LCD term dominant for the dynamical Hamiltonian through the additional control fields.
To confirm this, we plot the maximum amplitudes of both the first-and second-order adiabatic gauge potentials in Fig. 3, where Fig. 3(a) shows the case of no optimisation and Fig. 3(b) the case of applying COLD. We can see that without COLD the second-order (σ x j σ y j+1 + σ y j σ x j+1 ) corrections to the LCD are far larger than the first-order, resulting in the small final state fidelities when only first-order LCD is implemented. In the case of COLD, this relationship reverses and the first-order LCD terms dominate the dynamics. This gives us some indication that one way to optimise the control pulse may be a minimisation of higher order LCD terms, which we explore further in Sec. VII.
We find that the results of BPO and COLD at short driving times are stable against increasing system size, as shown for τ = 10 −2 J −1 in Fig. 4(a), with only a small decrease in final state fidelity for larger systems with COLD. Similarly, increasing the number of optimisation coefficients N k results in little improvement in the values obtained at short times for this example, as shown in Fig. 4(b). It is possible that in more complex systems, more optimisation coefficients will be needed to gain a larger advantage. We also note that by increasing the number of coefficients, we are increasing the complexity of the cost function landscape to be explored by the minimisation procedure, hence leading to slightly worse final fidelities This can mean that alternative approaches than the Powell minimisation used so far, e.g. that of CRAB, could be better suited to probing the cost function for high N k . We also note that this lack of improvement in the results is likely the consequence of the form of the control field given by Eq. (28) rather than due to a failure of the optimiser in the face of a complex parameter space. We find that the parameter space is relatively smooth in the case of N k = 1, 2, 3 and a better solution for this form of control field does not exist.
We now consider the combined method of COLD-CRAB for this annealing example as shown in Fig. 2(b). We point out that with our application of CRAB in this scenario we are not enforcing β to be zero at the start and end of the dynamics, allowing for their to be a tuning of the z-field offset. This is consistent between CRAB and COLD-CRAB and therefore does not influence our comparison of the two. First, it is important to note that CRAB alone results in a overall speedup of the dynamics for a high final state fidelity 1 − F ∼ 10 −3 at long timescales. However, CRAB still suffers from getting stuck in the initial state at fast driving times and the final state fidelity again tends to zero. This is not the case for COLD-CRAB, which converges to large final state fidelities 1 − F ∼ 10 −3 at short driving times τ ≤ 10 −1 J −1 . Note that the difference between the convergence to final state fidelities are only marginally different between COLD and COLD-CRAB at longer times, but at short times COLD-CRAB performs a lot better. Further improvement could be gained by combining COLD with more advanced versions of CRAB or other optimal control methods.
As shown in Fig. 3, the amplitude of the driving required to achieve the fidelities discussed so far scales with the driving time. Practical scenarios will necessarily place limits both on achievable driving times and the maximum amplitude of any term that is being driven. However, the scaling of the drivings shown do not mean that everything diverges in the limit of τ → 0. To see this we can first write the Scrödinger equation for COLD as we then divide through byλ to get in the limit of τ → 0 thenλ → ∞ to result in the Hamiltonian term disappearing, or in other words, we turn off the Hamiltonian. We then only drive the system in the τ → 0 limit with the COLD or LCD driving term: In this limit then λ plays the role of time, and this could then be implemented in a practical scenario in finite time as it corresponds to some manipulation of the couplings in the system. This renormalised time cannot then be infinitesimally short if the couplings are bounded but we have shown that the protocol does not diverge as τ → 0.
In the case of a spin chain, evolving under Eq. (33) is effectively to first order in LCD implementing independent single spin rotations along the chain, and COLD can be easily applied [40,41].
If it is not possible to switch off the Hamiltonian as discussed above then as an alternative we can implement COLD with experimental constraints accounted for directly in the optimal control minimisation. We consider an extreme example of constraints to show that even in this scenario COLD can provide an advantage and corresponding speed-up. In the constrained case the annealing protocol remains that of Hamiltonian (23) but we choose to introduce a bound of X f on the maximum amplitude of all drivings. This makes it so that no optimal control or LCD term can go beyond the original amplitude of the x-field drive. We show the final state fidelities achieved for the constrained example in Fig. 5. As can be seen in Fig. 5(a), COLD provides a substantial improvement beyond what is achievable with BPO. BPO manages F < 0.5 for τ < 1J −1 , but COLD can reach final state fidelities F ∼ 0.9 for τ < 1J −1 . The real improvement, however, comes with the application of CRAB and COLD-CRAB. CRAB already improves the fidelities substantially, and would allow for a speed up in the annealing protocol but with COLD-CRAB the final state fidelities are even better, with F ∼ 0.99 achievable when approaching τ ∼ 0.1J −1 . Signs are seen of the onset of the convergence to small values for COLD-CRAB in Fig. 5(b) before the maximum amplitude required becomes too large and the short time results tend towards zero fidelity and the state being stuck again. With this example and the discussion on implementation via turning off the Hamiltonian, we have shown that COLD is capable of delivering improvements beyond other schemes even for practical problems with strict and rather extreme constraints.

V. TRANSPORT IN A SYNTHETIC LATTICE
The efficient transfer of states between opposite ends of a lattice could have future applications in the settings of quantum computation and simulation due to its promise of efficient transport of information [42]. This objective is often tackled in the setting of ultracold atoms in optical lattices. While the problem can be tuned to be a singleparticle system and the analytical solutions of the corresponding instantaneous Schrödinger equation are known [43,44] even for a finite system [45], efficient evolution for state transfer is not straight-forward. This is due to the fact that the majority of the states are delocalised across the lattice, meaning that the |ψ ψ| terms of the CD Hamiltonian of Eq. (4) are global in reach. It is normal to consider this system in the tight-binding limit where the implementation of global terms is not straightforward. Such terms can be generated via the interactions of the atoms with cavity modes [46,47] or from dipolar interactions [48][49][50]. However, it would be ambitious to expect this control to be general enough to implement the CD Hamiltonian of the exact solutions. This is one of the reasons that LCD has been pursued in this setting.
Recently, LCD has been successfully applied to improve an adiabatic rapid passage (ARP) protocol for population transfer across a synthetic lattice [35]. In this realisation, population transfer was achieved in a synthetic tight-binding lattice of laser coupled atomic momentum states. We will consider the same problem as in Ref. [35] but with the improvement that can be gained by COLD.

This system is described by the Hamiltonian
where J n (t) is the time-dependent tunnelling that describes the nearest-neighbour coupling, V n (t) is the onsite energy offset with respect to neighbouring sites and c † n (c n ) is the creation(annihilation) operator on a given synthetic lattice site. In the ARP protocol, the population gets moved from one end of the lattice to the other by linearly ramping the lattice from a positive tilt to a negative tilt via where V 0 = 4J 0 is the initial site energy slope and J 0 is the characteristic tunnelling scale of the lattice. The scaling function in this case is given by In order to implement LCD as shown in Ref. [35], the first order LCD can be accounted for by taking where J n,CD (t) = J n (t) 2 + (α n (t)/τ ) 2 , and α n (t) is the CD terms which can be found by solving a set of linear equations In order to implement COLD we include additional terms to the tunnelling of the lattice which can then be incorporated into the forms of both J n,CD (t) and φ n,CD (t). We again want the additional control terms to go to zero around the problem boundaries and a natural choice is the same as in the Ising spin chain example in Eq. (28). The parameters β are optimised as before by minimizing with respect to the fidelity of the final state, where the population has been fully transferred to the opposite lattice site. We first consider a system size of N = 7 sites which was successfully experimentally probed in Ref. [35], where final state fidelities of 0.75 were achieved for τ = 1ms with a final tunnelling strength of J/ = 1/2πkHz (equivalent to τ ∼ 1J −1 in our units). We initially confirm the breakdown of ARP in this setting for fast times, and the success of the LCD protocol at short times, as shown in Fig. 6 (a) and found in Ref. [35]. Implementing BPO on its own manages to enhance the achievable fidelities at intermediate times of τ > 0.03J −1 . However, eventually, as observed in all scenarios in this work, BPO becomes stuck in the initial state at fast times, and the fidelity goes to zero. Implementing the newly introduced COLD protocol achieves an order of magnitude improvement in the fidelity over LCD. This is also plotted in Fig. 6(a) alongside previous results of ARP and first-order LCD.
One concern could be that COLD is achieving this improvement by simply pumping power into the tunnelling term, but as we can see in Fig. 6(c) the maximum amplitude of the tunnelling term tracks that of LCD. A key issue for experiments is the maximum amplitude achievable by a driving term and with this result we can stipulate that COLD is likely to be feasible in the same regimes as LCD in this synthetic lattice system. There is single outlier at intermediate times as indicated by the single point peaking in maximum amplitude in Fig. 6(c), this is the exception to the rule, where the optimisation has found a marginally higher fidelity (see the offset point in Fig. 6(a)) by pumping in power.
A large concern for state transfer techniques is the robustness of a protocol with respect to an increasing system size. We show the best achievable fidelities with increasing system size for both BPO and COLD in Fig. 6(c). While both protocols show a decreasing fidelity with system size as is to be expected, once again COLD does not suffer from getting stuck in the initial state. This is shown by the BPO fidelities going to unity for large systems in Fig. 6(c), and is the same mechanism for this as for the short driving times in Fig. 6(a).
Another concern could be that BPO will beat COLD if enough parameters are allowed for the optimisation, i.e. if we increase N k enough. We observed no evidence of this for the Ising model example and we again do not observe this in this synthetic lattice example, as is shown in Fig. 6(d). Small improvements are made in the fidelities achieved with BPO and COLD for larger N k but this is not substantial.

VI. GHZ STATE PREPARATION
As a final example we focus on the preparation of multi-partite Greenberger-Horne-Zeilinger (GHZ) [51] states: in a system of frustrated spins (see Fig. 7(a)). We start out with a system of all spins pointing down and drive a bare Hamiltonian of the form: where J = 1 and h = 10J with the same λ(t) as used previously, given by Eq. (12). The form of the LCD to first and second order is the same as in the case of the Ising spin chain (see Eq. (24) and Eq. (26)) with the couplings in the case of the second order now including the additional terms between spins j and j + 2.
In order to explore the versatility of combining optimal control with LCD, for this example we design our optimal control drive according to the gradient ascent pulse engineering (GRAPE) method [53]. Our control coefficients β k are now discretised on a finite grid of N m time intervals t m with uniform step ∆t to obtain control sequences in which individual elements β k,m (t m ) are treated as continuous parameters where the total driving time τ = N m ∆t and k is used to denote a localised drive for a subset of spins where N k is the total number of control pulses. As in the Ising  Fig. 2(a). Plotted in (b) and (d) are final state fidelities for τ = 0.1J −1 for different system sizes N when optimising β either by minimising the maximum amplitude of a drive max β (·) or its integral I2(·). In (b) only first order COLD is applied post-optimisation while in (d) one of the second-order drives is also applied after minimising the other, e.g. if max β |ζ| is minimised to determine the optimal β, then both the first-order drive α and the other second-order drive γ are applied. For system sizes above N=10 we used ITensor [52] MPS calculations which were converged with a truncation level of 10 −14 per time step. at each site reaching a maximum bond dimension of D = 4. In all cases, a single optimisable parameter is used (N k = 1). model case, we take our optimal control Hamiltonian to be of the form in Eq. (27) with each k th drive acting on the specified subset of spins with local σ z operators. At each time interval t m the k th control drive strength is calculated as: with θ(t) = sin πt 2τ and κ = 30 an offset parameter used to control the shape of the drive. We use spline interpolation to calculate the derivatives of the control drive when they are required to obtain the LCD drives. The resulting function requires more parameters than the Fourier basis we chose to use in previous examples, however it also allows for more flexibility in the final shape of the drive. Furthermore, due the increased number of parameters and search space, instead of Powell optimisation as in previous examples we choose to instead implement dual annealing, which is a global optimiser and while computationally more costly, is far better in the case of a complex parameter space with multiple minima.
Since such a preparation of GHZ states involves the generation of entanglement in a system that initially contains none, we expect that the first-order COLD may not contain the leading order of the counterdiabatic drive and thus not be as effective. For this reason, we include second-order COLD terms as given in Eq. (26). We also explore the idea of using multiple control drives and localising them to parts of the system. Thus we implement both a global drive which is uniform across all spins as well as a 'corner' evolution, in which three different optimisable drives are used: one each for the first and last spin in the lattice as well as one for all of the remaining spins. This is depicted in Fig. 7(a), where different vertex (spin) colours represent different control pulses.
The results for a 5 spin system with control drives consisting of N m = 10 time intervals are plotted in Fig. 7(c), where we observe that first order COLD is indeed not particularly effective at short driving times and does not move the system out of its initial state (see density matrix plots in (b)), regardless of whether or not separate control is applied to the corner spins. This is very likely due to the fact that the local σ y terms are only a small contribution to the full counterdiabatic drive and thus we need to look to higher order LCD to see any improvements. This is exactly what the results indicate, as second order COLD shows a five-fold improvement over the first order when a global optimisable drive is applied and up to two orders of magnitude improvement when the corner spins are driven separately at short times (τ = 0.001J −1 ). We then run the optimisations for larger systems at time τ = 0.1J −1 and find that this advantage is retained even with increasing system size.
This is a big improvement over recent results in digitized adiabatic evolution with LCD [54], where optimisation was used to determine optimal coefficients for second order LCD in order to prepare a GHZ state on an Ising spin chain. At 10 spins the final state fidelity for τ = 1J −1 obtained in their paper was 0.18, while we reach a fidelity of 0.72 for 15 spins when using corner optimisation at τ = 0.1J −1 .
This example shows that COLD can be used to speed up protocols which generate entanglement and is further evidence for the benefits of experimenting with different optimal control methods such as GRAPE as well as optimisation algorithms like dual annealing.

VII. MINIMISATION OF HIGHER ORDER LCD TERMS
As alluded to in Sec. IV, the results plotted in Fig. 3 indicate that in optimising the control pulse through the parameters β we maximise the largest amplitude of the first order LCD drive and simultaneously reduce the second order drives. In the Ising spin chain case this corresponds to increasing the largest amplitude of α(λ, β) in Eq. (26) throughout the evolution while reducing the maximum amplitude of both γ(λ, β) and ζ(λ, β). These results are a further indication that the implementation of COLD through the minimisation of the second-order corrections discussed in Sec. III may be fruitful in more complex and/or larger systems, where the dynamics can not be calculated.
We thus investigate replacing the original cost function of Eq. (2) with one that depends a) explicitly on the maximum amplitude of the second-order drives γ(λ, β) and ζ(λ, β) and b) one that depends on the total power for either drive. Given that the LCD drives are functions of β, one can imagine that if there is indeed a relationship between minimising a higher order drive and how effective the lower order drive is in producing the target state as a result, then we can determine parameters of the control drive that lead to a better final state fidelity.
We take the Ising Hamiltonian from Eq. (23) and supplement it again with the parameterised control pulse from Eq. (28). We once again take our first-order LCD drive to be of the form α(λ, β) j σ y j and the secondorder drives to be γ(λ, β) j (σ x j σ y j+1 + σ y j σ x j+1 ) and ζ(λ, β) j (σ z j σ y j+1 + σ y j σ z j+1 ). In Fig. 8(a) we show the results when the cost function used to optimise the parameters β is the integral from Eq. (22) which captures the total power of the drive: while in (c) we instead choose to minimise the largest amplitude of the drive reached throughout the evolution: In both cases we plot the resulting final state fidelities for different evolution times τ and compare them to those obtained earlier in Fig. 2(a) for 5 spins. The results are surprising in that while optimising for fidelity, as was done previously, outperforms second-order minimisation in both the integral and amplitude cases at most times, there is a stretch of driving times aroundτ ∈ [0.05, 0.5] where second-order minimisation does better. This can be attributed to the fact that the parameter landscape for the new cost functions is completely different and allows for a more optimal value of β to be reached without being lost in some sub-optimal minimum during the optimisation. In Fig. 8(b) we plot the final state fidelities at evolution time τ = 0.1J −1 for up to 50 spins in order to check how this type of optimisation scales with system size and to compare the performance of both cost functions. We find that minimising one of the two second-order drives while driving with the other still leads to impressive fidelities, but not as good as those where only first-order COLD is used. Indeed we do not have any reason to expect an absolute optimum fidelity when using this method, however, the results in Fig. 8 are very encouraging.
While the new cost functions in Eqs. (47) and (48) may seem like a roundabout way to get to the same result -a better final state fidelity in shorter time -they have several particularly important advantages over the cost function given by Eq. (2). First and foremost, this approach does not require access to the wavefunction or experimental data at any point of the optimisation process. In optimising for final state fidelity directly we must compute the evolution of the system many times over in order to extract the fidelity at each iteration, but computing the drive integrals or their amplitudes is completely independent of the state of the system. This allows us to determine an optimal set of parameters β for an arbitrary system size extremely efficiently when compared to methods which require access to |ψ f . A single optimisation in their case, depending on the method used and the desired quality of the final outcome, may take hours or even days for larger system sizes. The new method allows us to perform an optimisation with good results within minutes regardless of the number of spins, only requiring the wavefunction in order to check the resulting fidelity after the optimisation is finished. This is a very useful tool given that most optimal control methods demand access to the wavefunction while sacrificing efficiency.
It is not obvious that such a relationship between lower-and higher-order COLD as well as the fidelity of the final state must exist. In fact, this may be a fruitful new research direction to explore, combining the results obtained in this work along with, e.g. the methods in [38], where an approximate gauge potential can be systematically built up as a series of nested commutators. This might be a way to determine which operator ansatz O LCD has a maximal amplitude for each driven Hamiltonian and lead to a systematic optimisation of control pulses without ever having to simulate the system evolution. There is clearly a lot of new territory to explore both in terms of optimal control and in understanding the adiabatic gauge potential a little better.

VIII. DISCUSSION AND OUTLOOK
We have introduced a new hybrid approach combining quantum optimal control and shortcuts to adiabaticity: COLD. Inspired by the successes of LCD, where diabatic transitions are suppressed and locality conditions can be met, COLD improves on its methodology by combining it with quantum optimal control. The natural way to enhance the performance of LCD is by introducing higher order CD terms, but these are often non-local and difficult to engineer in experiments. COLD circumvents this by allowing for additional control fields that extend the family of dynamical Hamiltonians which can be explored. In this way, our method may find the best possible path where the effect of lower-order LCD is most relevant and higher order corrections are suppressed.
COLD has a clear potential in efficiently speeding up adiabatic evolution in various settings. We demonstrate this numerically via several example protocols which indicate improvements beyond a classical optimisation approach BPO as well as LCD of different orders. Our work shows that COLD reduces the strength of higher order LCD corrections, and that it performs well for increasing system sizes. We have shown that COLD can be implemented in the limit of fast driving by a 'switching off' of the original dynamical Hamiltonian. For scenarios where removing the Hamiltonian is not possible, we have shown that an alternative way to implement COLD is to use a bounded optimisation where amplitudes are restricted. We find that both the COLD and COLD-CRAB protocols perform extremely well in this setting.
COLD will be most beneficial when the LCD is only realisable to a certain order but the higher order corrections are large. This means the diabatic transitions are not being sufficiently suppressed by the choice of LCD and COLD can be used to find the dynamical Hamiltonian for which the required order of LCD term dominates. Note, that this goes the other way too, with COLD not providing substantial improvements when the chosen lower order LCD is small across the path. This can be thought of as being the case in two limits. First is the adiabatic limit, for which any CD correction is small and COLD will tend towards the adiabatic result. Second, the low-order LCD terms can be small compared to the driving as the exact CD would be correcting transitions due to interactions at longer ranges. In this scenario, the order of LCD being implemented with COLD needs to be increased, so that the CD term is accounting for the longer range terms. We show this in Sec. VI where it is clear that the generation of correlations or entanglement requires the suppression of diabatic terms that are non-local and thus first-order COLD cannot achieve a notable speed-up. In this case, higher-order corrections would need to be implemented with COLD, and finding methods for executing these non-local terms will be beneficial in these scenarios.
A further option is to combine COLD with one of a large variety of numerical optimal control methods, as we have done for the example of CRAB and GRAPE. We have shown a substantial improvement for state preparation in the Ising model that can be obtained from the COLD-CRAB combination -particularly in the constrained case. Fusions of COLD with advanced optimal control methods for complex systems could prove even more fruitful with further study.
Another finding of our work is that COLD can be applied to more complex systems where exact dynamics are not possible, e.g. due to an excessively large Hilbert space. This may be achieved by variationally minimising the integrals and maximum amplitude of the driving coefficients for the higher order corrections to the LCD. This opens up a brand new research direction as it allows for the possibility to optimise the system's path without requiring access to the system's wavefunction or any sort of experimental resource. Note, this finding is more general than COLD itself, as it can even be used to optimise protocols that do not implement LCD terms, i.e. the menagerie of control procedures currently in use, providing a cost function that does not scale with the system size. This would be implemented by minimising the highest orders of the LCD in order to find a path which allows for the least diabatic losses.
where we will consider each coefficient to be homogeneous across the chain and dependent upon the scaling factor of λ which is itself time-dependent as noted in the main text. We take the second order ansatz of the LCD to be that given by Eq. (26). We then want to obtain G λ as given by Eq. (6), which requires utilisation of standard commutation rules and the commutation relations of the Pauli matrices. Following several pages of working, the following form of G λ can be obtained Note that the three spin terms would trivially go to zero for the two spin example considered in the main text. As Pauli operators are traceless, we can easily compute the Hilbert-Schmidt norm of G λ and we simply need to keep track of factors from the finite size of the lattice to get where the factor on the LHS comes size of the Hilbert space. To find the system of equations to be solved we need to minimise Tr G 2 λ with respect to α, γ, and ζ to obtain If only the first-order correction of α is needed, then this can be obtained by taking the first equation and setting γ and ζ to zero. From this, the 2-spin and Ising model first order corrections can be obtained. Note in the limit of periodic boundary conditions or an infinite system we can take N → ∞ to obtain the correct coefficients. We find that the coefficients that are proportional to system size only have a significant impact when the system is very small, e.g. in the two spin case,and, therefore, have little impact the results of the Ising model with N ≥ 5.

Appendix B: Optimisation distributions
The results presented in Figs. 2, 5 and 6 of the main text contain plots of the best (highest) fidelities from a number of optimisations in each instance. Multiple optimisation runs with different initial guesses for the optimisable parameters are included to avoid pitfalls such as local minima in the parameter landscape.
In the context of a physical implementation of one of these protocols, the optimal set of parameter values (ones which return the highest fidelity with respect to the target state) matter more than the average. However, in practice these optimisations can be very computationally costly, in particular for larger system sizes and higher numbers of parameters. This means that we need to understand the behaviour of the average and the worst case as they relate to the computational resources required.
If the parameter landscape is smooth and few local minima exist, then only a few optimisations are needed to determine the best values of the optimisable parameters. However, this is never a guarantee and particularly in the case of the CRAB protocol (along with COLD-CRAB), the behaviour of the optimisation is suboptimal when it comes to the number of optimisations needed to determine the parameter values which return the best fidelity of the target state. This is due to the fact that we modify the parameter landscape for every optimisation by randomly changing the frequency components in the control field. While this allows each optimisation to access a new solution space and thus increases the chances of converging to a more optimal form of the control field, it also increases variance in optimisation outcomes. Since we cannot know which frequency gives the best results a priori, the only way to really reap the benefits of CRAB and CRAB-enhanced COLD is to perform as many optimisations as possible.
This can be readily seen when we look at the standard deviation in the final fidelities over all optimisations. Fig. 9 depicts these for the Ising spin chain of Sec. IV, both in the unconstrained and constrained case as well as for varying number of spins N and parameters N k . We can see that in most cases for COLD the standard deviation of the fidelities stays below 10 −3 barring longer driving times in the unconstrained case in Fig. 9(a) as well as some in (b) for the constrained case. BPO generally displays slightly higher standard deviations, but neither shows very significant variations in the results post-optimisation. Note that the small variation in fidelity for increasing number of parameters in Fig. 9(d) gives further evidence for the fact that additional parameters do not improve the results of COLD or BPO in the case of the Ising chain.
When it comes to CRAB and COLD-CRAB, however, the picture is quite different. We find that the resulting fidelities are a lot more varied across optimisations, as would be expected given the additional component of randomness. Fig. 10 shows not only the mean fidelities across optimisations but also the interquartile range of the data and the maximum and minimum values for each driving time. We find that across optimisations we are just as likely -and in some cases far more likely -to get a much worse final fidelity as we are to get a better one. This is reflected in the large range between the maximum (worst) and minimum (best) fidelity for both methods as well as the interquartile range, which shows that the mean fidelity is a result of a large variation between large and small fidelities rather than a convergence to some inbetween value.
These results are useful in an assessment of computational resources for such optimisations as well as giving an insight into the range of possible outcomes, particularly when implementing more unpredicatble optimal control methods like CRAB.

Appendix C: Choice of LCD Ansatz
In determining the optimal choice of operator basis O LCD , we turn back to Eq. (5) and note that it gives us some clues about the form of the LCD. Firstly, we note that if we know nothing about the system other than, say, that it is a spin chain described by Pauli matrices, then we take the first order LCD to be all one-body terms while the second order can be two body terms and so on. This is a natural choice due to the locality of the terms but also with respect to their practical implementation in an experiment. Given these considerations, it makes sense that for a system of spins,the first order LCD is a set of local σ y terms.
To illustrate, in the case of the Ising spin chain case, we know that all wave functions have real coefficients, so we know that the exact CD is given by entirely imaginary terms. We can confirm this by attempting to use local σ x or σ z terms as our ansatz for the operator basis O LCD and find that their coeffcients are equal to 0 throughout the driving time. In the case of ansatz O LCD = α x j σ x j we find: according to Eq. (6), Then the action, as in Eq. (8), is: which, when minimised with respect to α x gives α x = 0.
The same procedure can be done for O LCD = α z j σ z j : where again we take the action: which minimised with respect to α z once again gives α z = 0.