Reservoir Engineering using Quantum Optimal Control for Qubit Reset

We determine how to optimally reset a superconducting qubit which interacts with a thermal environment in such a way that the coupling strength is tunable. Describing the system in terms of a time-local master equation with time-dependent decay rates and using quantum optimal control theory, we identify temporal shapes of tunable level splittings which maximize the efficiency of the reset protocol in terms of duration and error. Time-dependent level splittings imply a modification of the system-environment coupling, varying the decay rates as well as the Lindblad operators. Our approach thus demonstrates efficient reservoir engineering employing quantum optimal control. We find the optimized reset strategy to consist in maximizing the decay rate from one state and driving non-adiabatic population transfer into this strongly decaying state.


I. INTRODUCTION
Superconducting qubits, combining sufficient isolation from the external environment and very good scalability, are a leading platform for demonstrating quantum advantage of a quantum computer [1].The ability to quickly and accurately reset qubits is a key factor for reaching the thresholds on state preparation and gate errors required by conventional quantum error correction codes.Tunable environments [2][3][4][5] provide a convenient way to initialize qubits on-demand.A method utilizing such a tunable environment to efficiently prepare superconducting qubits in their ground state has recently been brought forward [6].It exploits the indirect coupling of the qubit to a low-temperature resistive bath via two intermediate resonators [6] and use a protocol that utilizes sequential resonances with resistive bath.Here, we use quantum optimal control theory (QOCT) to study the efficiency of this reset protocol to further optimize it.
For a given model of a quantum system and its dynamics, QOCT provides a set of tools for obtaining the shapes of pulses which maximize a desired objective such as a gate or state preparation fidelity [7].In contrast to dynamical-decoupling-like approaches [8], QOCT does not rely on any a priori assumptions on the timescales of correlations functions of system and environment, and it allows for continuous dynamical modulation with minimal restrictions on the shape, duration, and strength of the applied pulse [7].In general, QOCT methods can be distinguished into those that evaluate only the objective functional such as the Chopped Random Basis (CRAB) method [9] and those that make use of also the gradient of the objective functional [7].The latter require both forward and backward propagation of the system dynamics and update the pulse shape either sequentially in time, such as Krotov's method [10], or concurrently for all times at once, such as the Gradient Ascent Pulse Engineering (GRAPE) algorithm [11].QOCT is particularly useful to study the control of open quantum systems since it allows to determine fundamental performance bounds due to decoherence and decay processes [12].Remarkably, the latter are not necessarily detrimental but may also be desired, for example when export of entropy is required to reach the objective [12].This is true for cooling in general [13][14][15] and for reset of qubits to a pure state in particular [16,17].
Utilizing the coupling to environmental degrees of freedom is also at the heart of quantum reservoir engineering [18] which deliberately incorporates dissipation into the system dynamics.In its simplest form, it is realized by a switchable, constant-amplitude electromagnetic field that drives transitions into a fast decaying state [18].For open quantum systems without memory, the system is driven into the fixed point of the Liouvillian, with constant and positive decay rates, that governs the dynamics [19,20].This idea has found widespread application in quantum optical experiments, for example with trapped atoms [21], ions [22,23] and circuit-QED platforms [24].For trapped ions, combining reservoir engineering with QOCT has recently allowed to determine the field strengths required to reach the error correction threshold in entangled-state preparation [25].For superconducting qubits, major decoherence arises from twolevel fluctuators which also render the dynamics non-Markovian [26].This can be captured by a strongly coupled environmental mode [27] or negative and timedependent decay rates in a master equation [28].How-arXiv:1903.05059v2[quant-ph] 23 Apr 2019 ever, reservoir engineering protocols have thus far been limited to exploiting decay with constant or piecewise constant rates [19,[29][30][31].
Here, we lift the limitation of constant decay rates by combining reservoir engineering with QOCT and a master equation featuring time-dependent rates.The latter are both controllable and experimentally implementable with current technologies [6].Using Krotov's method for QOCT [32], we derive the optimal shape of the external control fields that determine the time-dependent decay rates in the master equation.
The paper is structured as follows: In Section II we introduce both the Hamiltonian of the model and the main features of the used quantum optimal control method.In Section III, we present the numerical results for the optimization of the original protocol and compare it with the previous solution [6].Moreover, we extend the original protocol by adding two additional sets of control fields and evaluate the influence of the initial fields with which the optimization is started.Finally, in Section IV we summarize our findings and present conclusions.

A. Model
We consider a three-partite system consisting of two harmonic oscillators, named left (subscript L) and right (subscript R) oscillator, and a qubit (subscript q) as sketched in Fig. 1 and previously discussed in Ref. [4,6].We assume the two oscillators to be linearly coupled to each other through quadrature operators and the qubit to be exclusively coupled to the right oscillator.This scenario is modeled by the Hamiltonian (using units in which = 1) where â † L , â † R and σ+ are the creation operators for the left oscillator, right oscillator and qubit, respectively.The first three terms in Eq. ( 1) describe the free evolution of the subsystems, with ω q/L/R (t) being the timedependent and controllable level splittings of the qubit, the left and the right oscillator, respectively.The fourth and fifth term describe how the right oscillator is bilinearly coupled to left oscillator and to the qubit with time-dependent interaction strengths g LR (t) and g Rq (t), respectively.
The Hamiltonian (1) can be simplified by applying a rotating-wave approximation, assuming g 0 LR ω R , where g 0 LR is the resonant coupling strength between the oscillators.This results in [6] Within this approximation, the number of excitations is a conserved quantity in case of unitary evolution.Therefore, the total Hilbert space H of the system can be conveniently divided into subspaces H N where the number of excitations N is constant.A state belonging to a subspace H N will then remain within the subspace during the evolution that is solely governed by Hamiltonian (2).However, we consider the three-partite system to be open, interacting with an environment through one of its subsystems.Specifically, we take the left oscillator to be linearly coupled to a thermal reservoir.The system-bath interaction Hamiltonian is of the form [6] where VR is an operator of the reservoir and α plays the role of an effective coupling strength.In order to derive a master equation for the open system, we employ its instantaneous eigenbasis {|Ψ n (t) }, defined by , with ω n (t) being the respective eigenvalue.In this representation, the system-bath interaction can be rewritten as where v mn = Ψ m (t)| (â † L + âL ) |Ψ n (t) .Using standard techniques based on a weak-coupling hypothesis and the Born, Markov and secular approximations [33], it is possible to derive a Markovian master equation for the open system.The decay rates, responsible for dissipation and decoherence, are given by where ω mn (t) = ω m (t) − ω n (t) and S R (ω) is the real part of the Fourier transform of the reservoir correlation function, where the average . . .R is taken over the thermal state of the reservoir.The corresponding master equation in Lindblad form reads L † mn (t) Lmn (t), ρ(t) , ( 7) 1. Schematic diagram of the considered physical scenario consisting of a qubit (q) linearly coupled to a harmonic oscillator (R), which in turn is linearly coupled to a second harmonic oscillator (L) that is in direct contact with a thermal bath.By temporally controlling the level splittings ω q/L/R (t) of the qubit, the right and the left oscillator, one can effectively tune the coupling strength to the bath and change the decay rates over several orders of magnitude.
where the Lindblad operators Lmn (t) = |Ψ m (t) Ψ n (t)| describe transitions among the eigenstates.The Hamiltonian Ĥ(t) can be directly controlled by tuning the level splittings ω q/L/R (t).Importantly, the Lindblad operators and decay rates inherit the temporal dependence from the instantaneous eigenstates and eigenvalues.As a consequence, Eq. ( 7) goes beyond the standard description based on static decay channels with constant rates.
Solving the full master equation ( 7) is a rather challenging task and we therefore limit our study to a finite number of subspaces H N .Specifically, we consider the open dynamics of the system in the two subspaces H 0 and H 1 , i.e., the subspace with no excitations, H 0 = span{|0, 0, g }, and that with a single excitation, In the restricted Hilbert space, the Hamiltonian reads in the basis {|0, 0, g , |0, 0, e , |0, 1, g , |1, 0, g }.This simplified model can be solved analytically in the basis of the instantaneous eigenstates |Ψ 1 (t) , |Ψ 2 (t) , |Ψ 3 (t) ∈ H 1 and the ground state |Ψ 0 = |0, 0, g .Accounting exclusively for population decay from the excited states in H 1 to the ground state |0, 0, g , but not for the reverse process of thermal excitation, we obtain the following Lindblad master equation, cf.Eq. ( 7), where and the three time-dependent Lindblad operators are given by Closed form expressions for the exact eigenvalues ω i (t) and eigenstates |Ψ i (t) , albeit rather lengthy, are straightforward to calculate with computer algebra.

B. Physical realization
The model introduced above is quite general.In the following, we focus on a possible experimental realization which implies certain constraints and specific functional dependencies between the bare frequencies of the three subsystems, {ω L (t), ω R (t), ω q (t)}, and their respective couplings.The model described by Hamiltonian (1) can be realized by means of a superconducting qubit coupled to two LC resonators [6].The resonators behave effectively as quantum harmonic oscillators, the tunable frequencies of which are determined by the capacitance C and a controllable inductance L, i.e., ω L/R (t) = 1/ L L/R (t)C L/R .In this implementation, the couplings between the components can be expressed as functions of the physical parameters of the system and the bare resonator frequencies [6].
The reservoir is realized by connecting a resistor to the left resonator with VR in the interaction Hamiltonian (3) describing voltage fluctuations over the resistor.The resistor can be modeled as a thermal bath of bosonic modes [34], with the bath correlation function (6) corresponding to the Johnson Nyquist spectrum, where R is the resistance of the resistor and T its electron temperature [34].At low temperature, the spectral function (12) strongly suppresses emission of thermal excitations from the resistor so that indeed population decay is the leading-order dissipative process for the studied three-partite quantum system.The decay rates can be expressed as [6] Γ where Γ 0 plays the role of a static decay rate.

C. Quantum Optimal Control Theory
In general, QOCT aims at finding the optimal external control fields to steer the dynamics of a quantum system in the desired way [7].The starting point is to express the optimization task, for example the preparation of a specific target state, as a functional of the yet unknown external control fields {E k (t)}, where α τ is the final-time target functional which may depend on one or several states ρl (τ ) and g describes constraints that are relevant also at intermediate times, e.g.constraints on intensities or spectra [35,36].A proper choice of the functional implies that the extremum is only attained in case of the optimal dynamics.The final-time functional α T quantifies how well the target is reached.For the reset task at hand, we seek to prepare the qubit in its ground state, irrespective of the initial state of the total system.This can be achieved by considering the dynamics for several initial states {ρ l (t = 0)}, making sure that all of them result in the desired target state [37].Moreover, no excitation should be left in any of the two oscillators, since otherwise these might get transferred to the qubit in an uncontrolled fashion later on.Hence, our set of initial states {ρ l (t = 0)}, is given by any complete basis of the excited subspace H 1 .The respective target state is the ground state of the total system ρtrg = |Ψ 0 Ψ 0 |.The final-time functional thus reads where Â, B = Tr{ Â † B} and D(τ, 0; {E k }) is the controldependent dynamical map.Here, α τ corresponds to the error of the reset protocol, since it measures the remaining population in H 1 at final time τ .An ideal protocol is given by α τ = 0.In order to minimize Eq. ( 14), we employ Krotov's method [10,38], an iterative optimization algorithm that comes with the advantage of monotonic convergence.Given Eq. ( 15), a choice of g in Eq. ( 14), and equation of motion for the system, Eq. ( 9), Krotov's method provides a recipe to derive an analytic update equation for E k [32].Here, we use the standard choice of minimal amplitude increase per iteration step [39], where E ref k (t) is a reference field for each E k (t), taken to be the field from the previous iteration, S k (t) ∈ (0, 1] a shape function to smoothly switch the field modulations on and off, and λ k a parameter that controls the update magnitude of E k (t) in each optimization step.With Eq. ( 16), the update equation for E k (t) reads [32] where {ρ The so-called co-states { χ(i) l (t)} in Eq. ( 17) are solutions of the adjoint equation of motion, with boundary condition The derivative of α τ with respect to ρl (τ ) can be turned into a usual gradient by representing the states in a complete, orthonormal basis [39].Note that the indices (i+1) and (i) indicate values for current and last iteration, respectively.
As with any optimization algorithm based on variational calculus, Krotov's method requires the calculation of gradients-one of them the gradient of the dynamical generator with respect to the controls, cf.Eq. (17).Peculiarly, not just the Hamiltonian Ĥ, cf.Eq. ( 8), but also the dissipator L D , cf.Eq. ( 10) depends on the controls  15), as a function of protocol length τ for different control fields.SR denotes the original protocol utilizing sequential resonances with the resistive bath [6], CP refers to a protocol with only constant fields, and OP1, OP2, and OP3 are results obtained with SR or CP to start the optimization.An optimization targeting equal dissipation rates, cf.Eq. ( 22), instead of minimizing ατ is labeled by ER.The inset highlights the speedup due to the optimization, by comparing the durations for which the optimized protocols and the SR reach an error of 10 −4 .All parameters are summarized in Table I. and thus contributes to the gradient, Whereas the gradient of the Hamiltonian with respect to ω L , ω R and ω q is straightforward to calculate, cf.Eq. ( 8), the gradient of the dissipator L D is rather lengthy to evaluate, cf.Eq. ( 10).This inconvenience is due to the dependence of the decay rates Γ i0 and Lindblad operators Li on the instantaneous eigenvalues ω i (t) and eigenstates |ψ i (t) .The required derivatives of ω i (t) and |ψ i (t) with respect to ω L , ω R , and ω q have been algebraically calculated using computer software.

A. Optimization of the original protocol
The original protocol [6] is based on a simple choice for the left-oscillator frequency ω L (t).Effectively, it consists of two stages, namely ω L (t) = ω + and ω L (t) = ω − , separated by an intermediate ramp.The entire protocol reads where τ is the total protocol duration, τ /2 is the hold time at each stage, and t R τ is the ramping duration.The operation points ω + − and ω − have been chosen such that Γ 20 (t) = Γ 30 (t) in the case of ω + and Γ 10 (t) = Γ 20 (t) for ω − .This choice guarantees that any excitation decays at some point of time during the protocol.
Figure 2 illustrates the performance of the protocol of sequential resonances (SR) with the resistive bath as a function of its duration τ .It shows a rapid approach towards errors α τ as small as 10 −6 for τ = 2000 ns for the parameters listed in Table I.Although this may be sufficient for most applications, the SR exhibits a plateau for longer durations, preventing it even theoretically to reach significantly smaller errors.The plateau is caused by population being locked in the excited state of the right oscillator -an unfavorable feature that is apparently not resolvable by simply extending the protocol duration.However, taking Eq. ( 21) as the initial guess for the above described optimization procedure, Fig. 2 shows that, depending on τ , an improvement of up to two orders of magnitude in the error α τ compared to the SR is possible.In addition, the optimized protocol (OP1) also resolves the issue of the plateau, reaching errors α tau < 10 −7 .The improvement with respect to the protocol duration is comparatively modest, as the inset of Fig. 2 illustrates.Taking, e.g., α τ = 10 −4 as a sufficiently small error, the speedup with respect to the SR is roughly ∆τ ≈ 280 ns.
Figure 3 compares the decay dynamics of SR and OP1, for τ = 1500 ns, showing the population of the excited eigenstates in Fig. 3(a) and the respective decay rates and control fields ω L (t) generating them in Fig. 3(b) and 3(c).We observe that the original two-stage protocol (SR) acts as intended, i.e., the population decays from all three eigenstates of H 1 .Since the intermediate ramp transfers a significant amount of population from |Ψ 1 to |Ψ 2 , Γ 20 (t) needs to be sufficiently large also during the second stage.Note that this population transfer between different eigenstates within H 1 occurs due to non-adiabatic transitions caused by changes of those particular eigenstates.These are caused by changes in the control function ω L (t), i.e., the ramps in the SR.
A similar reasoning readily explains also the OP1 field shown in Fig. 3(c).Compared to the SR, the optimization effectively shifts the base levels of both stages and adds oscillations on top.This results in an increase of Γ 10 (t) and decrease of Γ 20 (t), cf.Fig. 3(b), in particular during the second stage, directly causing the population of |Ψ 1 (t) (|Ψ 2 (t) ) to decay faster (slower).The additional oscillations, even though having small amplitude, drive non-adiabatic transitions between |Ψ 1 and |Ψ 2 which primarily transfer population to the fast decaying state |Ψ 1 , cf.Fig. 3(a).This becomes even more clear by inspecting Figs.3(d) and 3(e), which show the spectra corresponding to the insets of the OP1 field in Fig. 3(c).In both cases, the frequencies match differences between various eigenvalues ω i , evaluated at ω + and ω − for ω L in the left and right inset, respectively.Whereas the spectrum shown in Fig. 3(d) is dominated by a peak at ω 2↔3 , which does not seem to have a notable impact on the dynamics, Fig. 3(e) solely exhibits a peak at ω 1↔2 and is responsible for the above-mentioned population transfer between |Ψ 1 and |Ψ 2 .The combination of increasing decay rates and engineered population transfer results in the excitation to more efficiently decay from both states.
The optimization studied in Fig. 3 changes the coherent part of the evolution compared to the SR, creating non-adiabatic transitions by suitably modulating ω L (t) and adapting the decay rates Γ i0 (t) accordingly.Both effects are necessary to explain the observed improvement with respect to the SR.In contrast, Fig. 2 shows also optimization results where the system dynamics has been completely ignored in the optimization process.In this case, the minimization of Eq. ( 14) has been replaced by a functional targeting equal dissipation rates (ER).Namely, we have optimized ω L (t) to yield R 1 ≈ R 2 ≈ R 3 with each R i as large as possible, where are the time-integrated dissipation rates which are independent of the actual dynamics.The naive assumption behind this optimization is that, since all states ρ1 , ρ2 , ρ3 are equally weighted in Eq. ( 15), equal dissipation from all of them may be a good choice to decrease the error α τ .However, this is not the case, cf.Fig. 2, which emphasizes the interplay of coherent and dissipative dynamics in the problem at hand.

B. Optimization with an extended set of control fields
In the following, we extend the SR by assuming the frequencies of the right oscillator and the qubit, ω R (t) and ω q (t), to be temporally dependent and controllable.Since the eigenvalues ω i (t) and eigenstates |Ψ i (t) (i = 1, 2, 3) depend on all three frequencies, ω L , ω R and ω q , changing any of them may affect the dynamics.In other words, more control fields give the optimization more flexibility to steer the system dynamics in the desired way and engineer the dissipation rates more appropriately.
First, we inspect in Fig. 4 how the decay rates change as a function of the level splittings ω L , ω R , and ω q .Two important observations can be made from Fig. 4. On one hand, the decay rates are still mutually exclusive, in the sense that there exists no combination such that two of them are large at the same time.On the other hand, the attainable maximum does not change.Hence, adjusting ω R or ω q in addition to ω L does not yield essentially larger rates, and there will not be a significantly faster decay to the ground state.Although no improvement is to be expected due to increased decay rates, i.e., due to the dissipative part of the dynamics, one may still achieve an improvement by more appropriately steering the coherent part.
Figure 2 shows optimization results for the case that all three frequencies are time-dependent (OP2).The initial guess has been chosen according to the SR, i.e., Eq. ( 21) for ω L (t) and constant values for ω R , ω q .Despite the extended set of controls, the optimization does not yield errors significantly below the case where only ω L (t) is controlled.This finding is reproducible even when using different sets of controls, such as only ω q (t) and ω L or only ω q (t) and ω R (data not shown).We therefore expect that further controls beyond ω L (t) do not allow the coherent part of the dynamics to be steered more efficiently.
In order to study this expectation further and evaluate the impact of the guess fields, we have carried out optimizations with all three possible controls.Whereas ω R and ω q have been set constant as initial guess, ω L (t) has been chosen as ω L (t) = (ω + + ω − )/2 with additional ramps in the beginning and end.Due to this choice, Γ 20 is almost maximal during the entire protocol, whereas Γ 10 and Γ 30 are orders of magnitude smaller, cf.Fig. 4. Thus, only population in |Ψ 2 decays fast.Simply extending the protocol duration τ will not solve the problem of small Γ 10 and Γ 30 .Upon optimization, we are, however, able to find fields yielding similarly small errors α τ as before (OP3) in Fig. 2. We again analyze an exemplary dynamics for τ = 1500 ns in Fig. 5.The panels are as in Fig. 3 with the small insets in (b) and (c) providing a closer look at the shapes of the optimized fields, respectively decay rates, compared to their non-optimized, constant counterparts.Panel (d) shows the spectra of all optimized fields from panel (c).
guess fields, while |Ψ 3 exhibits only slow decay and the population in |Ψ 1 is almost conserved.The respective decay rates and control fields are shown in Fig. 5(b) and 5(c).Interestingly, the optimization leaves the base levels of each control field unchanged, again adding small oscillations on top.Consequently, the decay rates are unchanged in magnitude but exhibit small oscillations as well.Since Γ 20 is already maximal by choice of the guess fields, cf.Fig. 4, there is no possibility for the optimization to increase it.Instead, the optimization ensures that all excitations are coherently transferred to this strongly decaying state -in our example from |Ψ 1 and |Ψ 3 to |Ψ 2 , as evident from Fig. 5(a).Thus, we find a similar reset strategy as in Fig. 3: The control fields are tailored such that a single decay rate (not necessarily the same at different times) is maximal and population is transferred coherently into this strongly decaying state.We expect the reset strategies illustrated in Figs. 3  and 5 to be feasible for essentially any combination of control fields and choice of guess fields.This follows from the decay rates being mutually exclusive, cf.Fig. 4, i.e., there is always one state that decays rapidly.All that is hence required is to ensure coherent population transfer into this state which seems to be possible by tailoring the control fields.Remarkably, the addition of further control fields does not result in significantly smaller errors α τ , cf.Fig. 2. In fact, ω L (t) alone is already sufficient to fully control the decay rates and engineer the required population transfer.Nevertheless, adding more control options increases flexibility and is thus potentially beneficial in experiments, especially if certain control fields are convenient to implement experimentally.

IV. SUMMARY AND CONCLUSIONS
In summary, we have studied how optimization of external control fields speeds up the initialization of a superconducting qubit which is tunably coupled to a thermal bath via two resonators.The control knobs are the time-dependent level splittings of the qubit and the resonators.Starting from a protocol utilizing sequential resonances with the resistive bath and employing the level splitting of a single resonator as the only control field, while assuming the initial state to be confined to the single excitation subspace [6], we have replaced the analytically derived temporal dependence by a numerically optimized control field.This has allowed us to obtain an improvement in both reset speed and fidelity.
We have also tested whether adding multiple control fields, by explicitly accounting for the tunability of the level splitting of the qubit and the second resonator, results in additional improvements.This has turned out to not be the case.Moreover, we have found that in all control scenarios, the optimized reset strategy consists in maximizing the decay rate from a single state and driving non-adiabatic population transfer into the strongly decaying state by small oscillations in the control fields.Even for different combinations of control fields and various guess fields, the optimization has resulted in reset errors and times of the same order of magnitude.We thus suspect to have identified the quantum speed limit for qubit reset in this particular physical setup with tunable couplings, provided only single excitations are present initially.However, a more rigorous study exploring the full parameter space is required to prove that our solution represents indeed a global, and not only a local, optimum.
Whether the quantum speed limit identified in our study is related to invoking the rotating-wave or other used approximations remains an open question.In particular, it will be interesting to study whether the reset duration and error can be further decreased by utilizing couplings between the single-excitation subspace and higher-excitation subspaces.The rationale would be that highly excited states decay faster which might further decrease the protocol duration.The required transitions could again be driven by suitably shaped control fields determined by QOCT.
Our study is, to the best of our knowledge, the first demonstration of experimentally directly applicable reservoir engineering using quantum optimal control of time-dependent decay rates.It is related to earlier results obtained for controlling open quantum systems with non-Markovian dynamics which had shown, for example, improved cooling due to cooperative effects of control and dissipation [15] or better gate operations [27,40].Our approach differs from the more common scenario for the control of open quantum systems in which the external field modifies only the Hamiltonian and thus the coherent part of the dynamics, rather than the dissipator of the master equation [12].In contrast, in our example, both the coherent evolution and the decay rates change as a result of the field optimization.Specifically, the changes in the coherent dynamics are manifested in the occurrence of non-adiabatic transitions which go hand in hand with modifications in the time-dependent decay rates.Interestingly, coherent and dissipative dynamics are tightly intertwined and the optimization protocol affects both in a physically transparent way.Our study thus paves the way to explore quantum reservoir engineering in condensed phase settings.
} are the forward propagated initial states spanning H 1 , obtained by solving

FIG. 2 .
FIG.2.Excited state population ατ , Eq. (15), as a function of protocol length τ for different control fields.SR denotes the original protocol utilizing sequential resonances with the resistive bath[6], CP refers to a protocol with only constant fields, and OP1, OP2, and OP3 are results obtained with SR or CP to start the optimization.An optimization targeting equal dissipation rates, cf.Eq. (22), instead of minimizing ατ is labeled by ER.The inset highlights the speedup due to the optimization, by comparing the durations for which the optimized protocols and the SR reach an error of 10 −4 .All parameters are summarized in TableI.

FIG. 3 .
FIG.3.Dynamics of the SR (dashed lines) and its optimized version OP1 (solid lines) for a protocol duration of τ = 1500 ns, cf.Fig.2.(a) Population in the three eigenstates of the excited subspace H1.(b) Decay rates, cf.Eq. (13), from H1 into the total ground state |Ψ0 .(c) Left oscillator frequency ωL(t) following the original two stage protocol of Eq.(21).The two stages are still visible in the optimized version, with modulations on top, as highlighted by the two insets.The shaded area in the left inset corresponds to fast oscillations, which are not resolved due to the linewidth.(d) and (e) show the frequency spectra of the optimized splitting ωL(t) from the left and right insets of (c), respectively.The vertical lines indicate frequency differences, ωi↔j = |ωi − ωj| with ωi = ωi(ωL) being the instantaneous eigenvalues.

FIG. 4 .
FIG.4.Decay rates Γi0 from the excited subspace H1 into the total ground state |Ψ0 , cf.Eq. (13), as a function of level splittings ωL and ωR and for three different values of ωq.

FIG. 5 .
FIG.5.Dynamics obtained with the constant protocol CP (dashed lines) and its optimized version OP3 (solid lines).The panels are as in Fig.3with the small insets in (b) and (c) providing a closer look at the shapes of the optimized fields, respectively decay rates, compared to their non-optimized, constant counterparts.Panel (d) shows the spectra of all optimized fields from panel (c).