Shortcuts to Adiabaticity Assisted by Counterdiabatic Born-Oppenheimer Dynamics

Shortcuts to adiabaticity (STA) provide control protocols to guide the dynamics of a quantum system through an adiabatic reference trajectory in an arbitrary prescheduled time. Designing STA proves challenging in complex quantum systems when the dynamics of the degrees of freedom span different time scales. We introduce Counterdiabatic Born-Oppenheimer Dynamics (CBOD) as a framework to design STA in systems with a large separation of energy scales. CBOD exploits the Born-Oppenheimer approximation to separate the Hamiltonian into effective fast and slow degrees of freedom and calculate the corresponding counterdiabatic drivings for each subsystem. We show the validity of the CBOD technique via an example of coupled harmonic oscillators, which can be solved exactly for comparison, and further apply it to a system of two-charged particles.


I. INTRODUCTION
Tailoring the nonadiabatic dynamics of quantum matter is an open problem at the frontiers of physics with important applications in emergent quantum technologies. Control protocols relying on adiabatic dynamics are natural to prescribe the evolution of a system along a reference adiabatic trajectory. While they are robust against uncontrolled errors in the experimental implementation, they are susceptible to decoherence. Driving protocols known as shortcuts to adiabaticity (STA) provide an alternative, by speeding up an adiabatic reference trajectory of a quantum system in a prescheduled amount of time [1].
STA have found broad applications in quantum mechanical systems of varying complexity. They can be used to guide the dynamics of systems with a discrete energy spectrum [2][3][4][5][6][7], as shown in the laboratory [8][9][10][11]. Similarly, STA can be used to control the degrees of freedom of continuous variables systems [1,[12][13][14][15][16][17][18][19] as demonstrated by the fast control of a trapped ion in phase space [20]. In the context of trapped ultracold atoms, early theoretical results indicated that STA could be applied to many-body systems [17,[21][22][23]. Ultrafast expansions and compressions of atomic clouds have by now been implemented in a wide variety of interaction regimes including thermal clouds [24], Bose-Einstein condensates well described by mean field theory [25,26], tightly confined quasi-one-dimensional atomic clouds with phase fluctuations [27], and a unitary Fermi gas as a paradigmatic instance of a strongly-coupled quantum fluid [28,29]. In addition, theoretical studies have shown that STA can be used to guide the evolution of many-body quantum systems that exhibit quantum critical behavior [30,31]. In this context, STA can be used to suppress excitation formation across a phase transition [32]. Implementing STA may require modifying the systems Hamiltonian with nonlocal interactions including high order terms [30,33]. However, the required controls may be simplified or absorbed into the form of the original system Hamiltonian [34][35][36][37]. Further efforts to control the dynamics of manybody systems have been put forward relying on integrability (e.g. the existence of Lax pairs) [38] or variational methods [39]. Despite this surge of progress, applications of STA remain mostly confined to systems with few degrees of freedom or the control of certain collective modes in many-body quantum systems.
Designing STA requires the ability to control and describe the time-evolution of a system, an ubiquitous challenge across a variety of fields when dealing with complex quantum systems. Among them, a prominent instance occurs in quantum chemistry, in the study of quantum systems with degrees of freedom spanning different time and energy scales. When the separation of scales is sufficiently large, it is possible to decouple the dynamics via the Born-Oppenheimer approximation (BOA) [40]. Born and Oppenheimer considered the description of a molecule and harnessed the separation of energy scales between the electronic and nuclear rotational and vibrational motions to simplify the description. As the electronic mass is much smaller than that of the atomic nuclei, the motion of the corresponding degrees of freedom occurs on vastly different time scales. When the mass ratio is large enough, the electrons move in an effectively static configuration of the nuclei. Further, this separation leads to the evolution of the nuclear component in the presence of a potential set by the energy of the electronic motion. When the assumptions of the approximation are relaxed, the nuclear motion is subject to a Berry vector potential due to the electronic motion [41][42][43][44][45]. Understanding the BOA and its limits has proved particularly fruitful in the field of spectroscopy [46,47], the study of molecular dynamics, and in computational models [48][49][50][51][52][53][54][55][56][57]. There has also been recent advances in the quantum simulation of molecular dynamics, with the recent proposal and experimental implementation of vibrational spectroscopy with trapped ions via boson sampling [58][59][60].
In this paper, we introduce Counterdiabatic Born-Oppenheimer Dynamics (CBOD) as an efficient technique for the fast control of complex systems that are well described by the BOA. To this end, we first provide a brief summary of the BOA and the engineering of STA by counterdiabatic driving in sections II and III, respectively. We then present CBOD in IV and demonstrate its validity with a paradigmatic example of coupled harmonic oscillators of unequal mass, in V.

II. BORN-OPPENHEIMER METHOD
In order to fix the notation, we will briefly discuss the Born-Oppenheimer approximation (BOA) that will be at the core of the CBOD technique. Throughout this work, we will consider two-body systems. However, the approach can be readily generalized to N -body systems with two sub-sets of slow (heavy) and fast (light) particles, that we indicate with the labels S and F , respectively. We will make the assumption that m S m F . In this section, we will introduce both the conventional BOA as well as the relaxed BOA, whereby, the fast degrees of freedom give rise to a Berry vector potential for the slow variables.

A. Conventional Born-Oppenheimer approximation
Consider a Hamiltonian of the form wherep i = −i ∇ i (i = S, F ) denotes the momentum operator and V (x S ,x F ) is a global potential term. The latter can generally be decomposed as with the first two terms acting exclusively on the slow and fast coordinates, respectively, and an interaction term, which is not separable in the coordinate representation in terms of {x S , x F }. The general form given by this Hamiltonian includes the usual molecular Hamiltonian for electronic and nuclear dynamics, frequently used in quantum chemistry [40,46,61]. We wish to obtain the solutions to the Schrödinger equation Given the difference in mass (m S m F ), the BOA suppresses the kinetic energy term of the slow variables [62] to obtain a reduced Hamiltonian for the fast ones. Aŝ x S commutes with this reduced Hamiltonian, we can simultaneously obtain the solutions of the reduced and full sub-systems. The reduced Hamiltonian governs the Schrödinger equation of the fast sub-system The slow coordinates x S can be regarded as a parameter on which the reduced system eigenvalues and eigenvectors depend . The solutions to the reduced Schrödinger  equation form a complete set, in terms of which the full  solution to the complete Schrödinger equation, Eq. (3), can be written as where n runs over the eigenstates of the reduced Hamiltonian. We will assume in this work that the fast subsystem is in a single eigenstate n, avoiding the need for the summation in Eq. (5).
Using the product expansion of the wave function, Eq. (5), and the reduced Schrödinger equation, Eq. (4), the full Schrödinger equation reads In the conventional BOA the derivatives of the fast subsystem wave function, φ (x F ; x S ), with respect to x S are neglected in the above equation, i.e., Integrating out the fast degrees of freedom is then straightforward and leads to a slow sub-system Schrödinger equation in the final form Note, E gives the full energy of the system while the full approximate wave function is given by Eq. (5). The conventional BOA involves truly two approximations: (1) the energy scales of the system are vastly different allowing for the suppression of one of the kinetic energies to obtain the reduced Hamiltonian and (2) corrections due to the elimination of derivatives in x S of the reduced wave function are small. These corrections are referred to as diagonal corrections and they are usually negligible in comparison to the energy scale of the fast sub-system [46]. The form of these corrections and their calculation is a vibrant area of research in its own right [63][64][65][66][67][68][69][70]. We have added another approximation to the conventional approach, that the fast sub-system evolves adiabatically, this was invoked when the sum over the fast sub-system states was neglected. This is a common approximation when using the BOA to describe timeevolution [46,71], as the fast sub-system is assumed to quickly relax to its ground state in the time-scale of the slow motion. Note, that when we combine the BOA with counterdiabatic driving, we will assume the fast subsystem either evolves adiabatically or, more importantly for our approach, that the fast sub-system is driven such that adiabaticity is enforced; in either case, the adiabatic approximation is met.

B. Relaxed Born-Oppenheimer Approximation
The approximation involving the elimination of the derivatives with respect to x S made in the conventional BOA can be relaxed [43,44]. Generally, the neglected diagonal term can couple arbitrary eigenfunctions of the fast degrees of freedom. A relaxed BOA consists of keeping the resulting cross-terms of the x S derivative while neglecting transitions between these different eigenstates. This leads to the appearance of a Berry connection between the two sub-systems in the full Schrödinger equation, which plays a role analogous to the vector potential in the quantum mechanics of a charged particle in an electromagnetic field [43,44,72].
Starting from the Schrödinger equation (6), our goal is to obtain a Schrödinger equation which is solely dependent on the slow degree of freedom, without invoking the previous approximation that disregards the derivatives of the momentum operator. To remove the fast degree of freedom we multiply Eq. (6) from the left by φ (x F ; x S ) † and integrate over x F to obtain There is only one non-trivial integral in the above Schrödinger equation, which is, It is an algebraic exercise to obtain the terms arising from this integral, which can be written in terms of a vector and a scalar potential in the slow sub-system Schrödinger equation with The vector potential A (x S ) is the familiar Berry connection. The scalar potential ε n (x S ) is local and is dictated by the fast variables. In addition, there is a contribution to the scalar potential experienced by the slow degrees of freedom that is given by g(x S ), which is the trace of the quantum geometric tensor [73] associated with the change of the eigenstates |φ with respect to the slow coordinates x S , treated as a parameter.

III. COUNTER-DIABATIC DRIVING
Among the variety of techniques available to engineer STA, counterdiabatic driving (CD) stands out as a universal approach. It relies on the use of auxiliary counterdiabatic fields to guide the evolution of the quantum system of interest through an adiabatic reference trajectory. CD was developed in the context of molecular dynamics by Demirplak and Rice [2,3,74], as an alternative to strictly adiabatic population transfers between molecular states; see also the independent and closely related work by Berry [4]. CD and related protocols have been recently implemented in a variety of platforms for quantum technologies including trapped ions [20], nitrogenvacancy centres in diamond [9,11], ultracold atoms in optical lattices [8] and as a method to speed-up stimulated Raman adiabatic passage in ultracold gases [10]. It has become a popular technique to control and engineer the nonadiabatic evolution of quantum systems while enforcing the following of adiabatic trajectories [2-4, 23, 39, 74].
Counterdiabatic driving relies on the spectral properties, eigenstates and energies, of the driven Hamiltonian of interestĤ According to the adiabatic approximation, the state of a system prepared in an eigenstate |n (0) at t = 0 evolves under a slowly-varyingĤ 0 (t) into which includes the dynamical phase as well as the geometric phase associated with the Berry connection i n (t ) |∂ t n (t ) . A STA protocol assisted by CD can be designed by identifying a modified driven HamiltonianĤ (t) such that the adiabatic evolution (15) becomes the exact solution of the corresponding time-dependent Schrödinger equation Hence, no matter how fast the system is driven, the evolution is described by the adiabatic trajectory (15), i.e., without the requirement of slow driving. The corresponding time-evolution operator also fulfills (16), which allows the identification of the modified driven Hamiltonian as the generator of evolution Making use of the following form of the time-evolution operatorÛ the modified driven Hamiltonian is found by explicit com- where we have defined The first term in (19) is recognized as the spectral decomposition of the original system HamiltonianĤ 0 (t). The second term,Ĥ 1 , is the auxiliary CD term required so that the adiabatic trajectory |ψ ad n (t) in Eq. (15) becomes an exact solution of (16), that is the Schrödinger equation for the full driving hamiltonianĤ =Ĥ 0 +Ĥ 1 .
When the energy spectrum ofĤ 0 is non-degenerate, the additional CD term can be recasted using the differential of the time-independent Schrödinger equation of the original system Hamiltonian [4],Ĥ 0 (t), which yields the following alternative expression for the auxiliary CD term in terms of the projectorP m (t) = |m (t) m (t)|.

IV. COUNTERDIABATIC BORN-OPPENHEIMER DYNAMICS (CBOD)
The Born-Oppenheimer method and the theory of counterdiabatic driving can be exploited jointly to engineer the superadiabatic control of complex systems, as we next discuss. The counterdiabatic drivings for the slow and fast sub-systems can be obtained via the BOA, either via the conventional or relaxed variants. These auxiliary control terms can then be used to drive the (exact) system Hamiltonian, a technique we shall term as Counterdiabatic Born-Oppenheimer Dynamics (CBOD).

A. CD with the conventional Born-Oppenheimer approximation
We first consider the conventional BOA, according to which the fast and slow sub-system Hamiltonians arê The required CD terms can be found via the general expression Eq. (20) and for the slow and fast sub-systems read respectivelŷ By contrast, the CD term for the full system iŝ The full form can be rewritten, exploiting the tensor product structure of the full wave function |Ψ = |φ ⊗|ψ (i.e. separable in the BOA). We note that while this separable form is natural in the BOA, it can be invoked generally [75,76]. Substituting in the factored form of the full wave function, the full CD can be written aŝ Therefore, the CBOD technique simplifies the global CD control by driving the two sub-systems separately, this is, by using the auxiliary control terms (25) and (26) as opposed to (27). We can gain further insight into the CBOD terms by assuming the spectra ofĤ S andĤ F to be nondegenerate, as this allows the recasting of the CD control terms into the form of Eq. (22), Note, that in the fast controlĤ F,1 , x S is treated as a parameter, which is the case of the reduced Hamiltonian in Eq. (4) in the BOA. As customary in STA assisted by CD, the required auxiliary terms are off-diagonal in state space. However, it is useful to focus on the dependencies on x S and x F . We wish to compare the CBOD terms with the exact CD term without resorting to the Born-Oppenheimer approximation. Assuming the system is exactly solvable the CD isĤ withĤ 0 |χ n = n |χ n , |χ the exact eigenstates and n the exact eigenvalues. Due to the interactions between slow and fast degrees of freedom, which are not separable in the {x S , x F } space, the operatorĤ exact,1 can be complicated to implement, as it involves a generally complex coupling between the two sub-spaces (or two particles).
Therefore, CBOD provides a computational advantage over the exact STA approach. Indeed, it circumvents the need to deal with the full spectra, proceeding, instead, in two-subsequent steps; treating first the fast degrees of freedom and then the slow ones. CBOD thus benefits from the dimensional reduction of the problem to engineer the CD term for the fast sub-system driving, Eq. (29).
In addition, CBOD may simplify the required drivings by potentially removing or reducing the coupling between the two sub-spaces. In such a case, CBOD controls will be simpler to implement than the exact CD terms. For the slow sub-system, Eq. (30), the off-diagonal terms are coupled via the time derivative of the energy for the fast sub-system. This is the potential surface which the slow sub-system experiences due to the fast sub-system, and it is of no surprise that to enforce adiabatic evolution it is required to drive off-diagonal terms with this as the scaling. The coordinate dependence of slow sub-system counterdiabatic driving, Eq. (30) is simplified in comparison to the exact driving, Eq. (31), as the potential energy surfaces can only depend on x S . Therefore, the slow sub-system driving only requires operators which act on the slow sub-system space and will have no cross terms between the two sub-spaces. Therefore, CBOD readily simplifies the required control term for the slow sub-system, without further approximations.

B. CD with the relaxed Born-Oppenheimer approximation
In what follows we derive the modified driving controls when the relaxed BOA is used, as discussed in Sec. II B. The required CD terms for the fast and exact systems are the same as in Eqs. (25) and (27); alternatively, by (29) and (31). The slow sub-system has a modified Hamiltonian resembling that of a particle in an electromagnetic field, see Eq. (11).
To obtain the CD under the assumption of no degeneracies in the spectra according to Eq. (22), we first need to obtain the time derivative of the Hamiltonian The latter admits the compact form with the potential term given by where {·, ·} denotes the anti-commutator andȦ ≡ ∂ t A.
The slow sub-system counterdiabatic driving is therefore given bŷ whereP m,S = |m S m S |. As a result, even with the relaxed BOA, the slow sub-system CD term only depends on operators related to the slow coordinate x S .

C. Applicability of Counterdiabatic Born-Oppenheimer Dynamics (CBOD)
CBOD, as an approximate technique, does not necessarily enforce the evolution of the system Hamiltonian to follow the adiabatic trajectory exactly. It resorts to the counterdiabatic driving terms constructed via the BOA to drive the (exact) system, which includes couplings between slow and fast subsystems beyond BOA. Said differently, CBOD is constructed to drive the fast and slow Hamiltonians of the BO Hamiltonian (as opposed to the exact system Hamiltonian) exactly through the adiabatic manifold. We will assess the validity of CBOD using the fidelity between the resulting state and that of the exact adiabatic evolution after a modulation of the system Hamiltonian in a prescheduled time.
The implementation of the CBOD technique, in general, involves the following steps: 1. Check the validity of BOA, i.e., that there are two separated energy scales in the region of interest.
3. Apply these (approximate) control terms to guide the dynamics of the (exact) system Hamiltonian.
We will consider next an example discussing each step in detail and certify CBOD by comparing its performance to the exact CD evolution.

V. COUPLED HARMONIC SYSTEM
To illustrate CBOD, we next consider the engineering of STA to drive two coupled harmonic oscillators, that can represent, e.g., two atoms in a harmonic trap interacting via a spring-like term. This model has been previously used to assess the BOA [77], admits an exact solution [78][79][80][81][82][83] and is realizable in controllable quantum systems of ion traps [84][85][86][87]. It, therefore, constitutes a natural test-bed for CBOD.
Specifically, we consider the Hamiltonian with continuous variablesx S andx F in one spatial dimension. Alternatively, it can be rewritten aŝ with κ S = k S + k I and κ F = k F + k I , which makes explicit the bilinear coupling. For the sake of generality, we first provide a derivation of the counterdiabatic drivings for this system when all spring constants are time-dependent. The spectral properties can be studied by diagonalizing the system in terms of two independent harmonic oscillators, the normal modes. We denote byŷ i ,p i , and κ i (i = 1, 2) the corresponding normal-mode coordinates, conjugate momentum and spring constants, for which explicit expressions are derived in Appendix A. In terms of them, the system Hamiltonian can be simply written aŝ The BOA leads to an approximation of the system Hamil-tonianĤ 0 also in terms of two independent harmonic oscillators, whose eigenstates under the conventional and relaxed BOA coincide, as the Berry connection identically vanishes and the quantum geometric tensor reduces to a time-dependent constant; see Appendix A for further details. Note, that in order for the slow sub-system of the BOA Hamiltonian to have a real harmonic frequency it is required that κ S (t) κ F (t) > k I (t) 2 . Knowledge of the exact and BOA eigenfunctions allows us to establish the validity of the BOA whenever m F /m S 1. To this end, we consider the fidelity between an exact eigenstate of the two coupled harmonic oscillators Ψ exact and the corresponding BOA Ψ BOA The Born-Oppenheimer method provides a good approximation for the coupled harmonic oscillator problem for a small mass ratio m F /m S , see Fig. 1. As the trapping frequency of any of the two sub-systems -slow or fast particles -is increased, the accuracy of the BOA increases, as shown in Fig. 1(a)-(b). This is consistent with the fact that the energy scale separation between the two subsystems is increased for a given mass ratio as the subsystem spring constants are made larger, making the two sub-systems more decoupled. By contrast, when decreasing the interaction strength, see Fig 1(c), the state Ψ BOA approaches Ψ Exact as quantified by the higher values of the fidelity. Naturally, the coupling between both subsystems increases with the interaction spring constant, leading to a breakdown of the BOA at large values of k I . We next derive and compare the auxiliary control terms required to enforce adiabaticity in an arbitrary prescheduled time using the exact CD and CBOD.

A. Exact Counterdiabatic Driving
The exact solution of the coupled system (36) can be written in terms of the two independent harmonic oscillators, the normal modes, described by Hamiltonian (38). Knowledge of the CD term for a single harmonic oscillator [23,88] readily yields the exact CD term for the coupled system as the sum of the generators of the squeezing operator for each normal mode. As such, they are spatially non-local due to the momentum dependence. Alternative controls can be obtained by means of the unitary transformation which acts on the position and momentum operators in the Hamiltonian aŝ Given that ∂ t U † = 0 the full driving HamiltonianĤ = H 0 +Ĥ 1 is transformed according tô  while the original wave function Ψ is mapped to Making use of (42), it is found that the transformed HamiltonianĤ T , unitarily equivalent toĤ, takes the form with the corresponding frequencies being . (46) Therefore, the exact CD of coupled harmonic oscillators can be implemented by a modification of the normal mode frequency of the original system Hamiltonian H 0 . However, this modification will require independent control of the slow, fast and interaction spring constants.

B. CBOD
Within the BOA the Hamiltonians of the slow and fast subsystems are that of two harmonic oscillators, and the corresponding CBOD terms are given bŷ wherex T = x F − k I (t) κ F (t) x S andp T =p F is the corresponding conjugate momentum operator, see appendix A. For the coupled harmonic oscillators the CBOD auxiliary controls under the conventional and relaxed BOA are equivalent, as the wave functions coincide. We will consider a case in which the fast sub-system is also driven to enforce adiabaticity by CBOD. However, within the BOA the fast sub-system is usually assumed to evolve adiabatically and, in that case, the controlĤ F,1 would not be required. The driving Hamiltonian with the CBOD control terms readŝ The evolution under this Hamiltonian is not necessarily adiabatic with respect to the exactĤ 0 eigenbasis, as the CBOD auxiliary terms are approximate. As a result, an STA designed by CBOD can not be arbitrarily fast. The direct exact solution of the Schrödinger equation with this Hamiltonian is hindered by the termx SpF resulting from the last anti-commutator. However, this term can be absorbed into the fast momentum. Dropping the time dependence for simplicity, one findŝ where in the final line we have made an approximation consistent with the BOA, that the fast sub-system momentum will dominate over the additional momentum term which is a function of the slow sub-system coordinate. The driving Hamiltonian takes then the form In a similar manner to the previous scenario, we can use the unitary transformation of to obtain the unitarily equivalent Hamiltonian which is that of two oscillators with bilinear coupling, i.e., the original system Hamiltonian (51), with modified spring constants Therefore, CBOD simplifies the engineering of STA in the system by driving the slow and fast sub-systems independently. In addition, it succeeds in doing so without the need to tailor the interaction term between the twosubsystems. This is contrary to the exact counterdiabatic driving, which involves a controlled modulation in time of all the potential terms in the original Hamiltonian, including the interaction. As CBOD relies on the BOA, the dynamics generated by the Hamiltonian (53) is not strictly adiabatic. The exact solution to the corresponding time-independent Schrödinger equation can be obtained in a similar manner to that shown for the exact solution in Appendix A, with a separation into two independent normal-mode harmonic oscillators. The non-adiabatic evolution of each normal-mode harmonic oscillator can then be described exactly by a self-similar transformation of the corresponding wave function Φ at the start of the evolution. For a harmonic oscillator of mass m and frequency ω (t) the scaling symmetry determines the evolution of the ground state according to [1,23,88,89] where ω 0 = ω (0) and b (t) > 0 is a scaling factor obtained by solving the Ermakov equation with boundary conditions b(0) = 1 andḃ(0) = 0. By solving the above Ermakov equation numerically for the parameters of Hamiltonian (53), the exact evolution of the system under CBOD can be obtained via Eq. (56).

C. CBOD Validity
To investigate the validity of the CBOD technique we consider the following modulation of the spring constant where k 0 and k 1 are offset and strength parameters of the ramping, respectively, and T f is the time-scale of the modulation, with a single modulation after t = T f . This ramp has first and second derivatives that vanish at t = 0 and T f , favoring adiabatic evolution [90], e.g., over a linear ramp. We will consider a ramping of the individual terms in Hamiltonian (36), i.e. κ S , κ F and k I . We show the fidelity of the ground state under a ramping of each spring constant in Fig. 2, i.e.
From the assessment of the validity of the BOA in this example, see Fig. 1, the CBOD technique is expected to be particularly sensitive to ramping the interaction spring constant. A significant drop off in fidelity is observed as the interaction strength is increased by the ramp. Overall, CBOD matches with high fidelity the exact adiabatic evolution. This is reflected by the values F ≥ 0.99 observed in Fig. 2.
For large ramps of the interaction spring constant, the dotted green lines of Fig. 2, a slight breakdown of the validity of CBOD for m F /m S ∼ 1 becomes manifest. This is most likely due to the momentum approximation of Eq. (50), which was required to decouple the fast momentum and slow coordinates in the CD control terms within the BOA. This approximation is only valid in the limit of m F /m S 1 and/or k I κ F , which are both broken in this scenario. Hence, the lower fidelity of the CBOD protocol for large mass ratio. The breakdown of the CBOD for m F /m S ∼ 1 is consistent with the breakdown of the BOA. Provided that the fast and slow degrees of freedom can so be defined, CBOD can generate high fidelity evolution of states under fast modulations.
However, by contrast to exact CD, CBOD is not valid for arbitrarily fast modulations. CBOD trades the possibility of engineering arbitrarily fast STA for the ability to treat interacting systems which may not be exactly solvable and to simplify the experimental implementation of the required control terms in these systems. The behavior of CBOD under faster modulations is investigated in Fig. 3, by the fidelity of the exact and CBOD techniques at the end of the process t = T f . As would be expected, the fidelity decreases with faster ramping times, T f but the fidelity remains high (F ≥ 0.9) and this trend continues for smaller T f , for which the fidelity exhibits a plateau. In Fig. 3, each ramp is of a moderate strength and as this strength is increased the fidelity decreases, as was shown in Fig. 2. Therefore, despite the dependence of the fidelity on the modulation time, CBOD can be used for fast driving of the system.

VI. CONCLUSIONS
Shortcuts to Adiabaticity (STA) provide control protocols to guide the dynamics of quantum and classical systems along an adiabatic reference trajectory, without relying on slow driving. A universal approach to designing STA is provided by the counterdiabatic driving (CD) technique that guides the evolution of an arbitrary quantum system by means of auxiliary control fields. However, determining the auxiliary controls requires knowledge of the spectral properties of the system, hindering its application to complex systems.
In this work, we have introduced Counterdiabatic Born-Oppenheimer Dynamics (CBOD) as a framework to design STA in complex systems. CBOD identifies the required controls to speed up the dynamics of the system by invoking the Born-Oppenheimer approximation (BOA) whenever a separation between fast and slow degrees of freedom is justified. In such a scenario, the required CD terms for the fast and slow variables can be obtained in two subsequent steps, which in the spirit of the BOA, avoids the need to diagonalize the high-dimensional Hamiltonian of the full system. Thus, CBOD facilitates the finding of the required Hamiltonian controls to speed up the dynamics, in scenarios where spectral properties are not readily available. In addition, CBOD also simplifies the implementation of the STA by reducing the need to control the coupling between fast and slow degrees of freedom. We have demonstrated the validity of CBOD by testing it in a paradigmatic test-bed of BOA, an exactly-solvable model of two driven coupled harmonic oscillators with unequal masses for which CBOD competes with the exact counterdiabatic driving in the preparation of a target state. We anticipate that the CBOD technique should facilitate the control dynamics of complex systems assisted by STA in the plethora of scenarios in which the Born-Oppenheimer approximation has proved useful. The exact solution to the Schrödinger equation of Hamiltonian (37) is well known [79]. First, we transform the position and momentum spaces canonically via the transformations and This is followed by a rotation of the coordinates under which the momentum is invariant. To diagonalize the system the rotation angle is found to be The Hamiltonian is then diagonalized into normal modes, i.e., two independent harmonic oscillatorŝ H 0 (t) =p 2 1 2µ +p with reduced mass µ = √ m S m F and spring constants The solution to the time-independent Schrödinger equation of Hamiltonian (A5) is that of two independent harmonic oscillators in the coordinates y 1 and y 2 with total energy ij (t) = ω 1 (t) i + with i, j = 0, 1, 2, . . . and frequencies ω 1,2 (t) = κ 1,2 (t) /µ. The full wave functions take the form of the usual Harmonic oscillator solutions, i.e. ψ ij (y 1 , y 2 ) = 1 With the exact driving utilised in the main text the time evolution will be the adiabatic solution ψ ij (y 1 , y 2 , t) = exp − i t 0 dt ij (t ) ψ ij (y 1 , y 2 ) . (A10)

Born-Oppenheimer approximation for two coupled harmonic oscillators
We now consider Hamiltonian (37) in the regime of m S m F , where the BOA is valid. Following the steps of Sec. II A, the time-independent Schrödinger equation of the fast subsystem reads with transformed coordinatex T = x F − k I (t) κ F (t) x S , momentump T =p F and slow sub-system spring constant κ S (t) = κ S (t) − k I (t) 2 κ F (t) . Eq. (A11) has solutions of a harmonic oscillator in the x T coordinate, with eigenvalues ε n (x S ) = ω F (t) n + 1 2 and frequency ω F (t) = κ F (t) /m F . The slow sub-system then follows the Schrödinger equation of with frequency ω S = κ S (t) /m S . Note, that κ S (t) can be negative, turning the frequency imaginary and the solutions considered in this work incorrect. We will take care to ensure that we remain in the real frequency limit, i.e. κ S (t) κ F (t) > k I (t) 2 . The total energy of the system will be and total wave function for a single modes is

Relaxed Born-Oppenheimer approximation
We next consider the application of the relaxed BOA to the coupled harmonic oscillators. The Berry connection in Eq. (12) and geometric tensor given by Eq. (13) are both integrals of the fast (or reduced) sub-system wave functions, φ u , which are given in Eq. (A12). In particular, we consider the fast sub-system to be in the ground state, The Berry connection identically vanishes and the geometric tensor simply reads (A20) Excited states of the fast sub-system also result in a vanishing Berry connection and a geometric tensor of similar form to above. Therefore within the relaxed BOA, the coupled oscillator problem has the same wave functions as the conventional BOA but with a different total energy. That difference is the value of the geometric tensor multiplied by a factor, as shown in Eq. (11). In this example, the relaxed and conventional counterdiabatic drivings are identical, as the wave functions coincide.