Failure of the geometric approach prediction of excess work scaling for open and isolated quantum systems

The task of finding optimal protocols that minimize the energetic cost of thermodynamic processes of long yet finite duration $\tau$ is a pressing one. We approach this problem here in a rigorous and systematic fashion by means of the adiabatic perturbation theory of closed Hamiltonian quantum systems. Our main finding is a $1/\tau^2$ scaling of the excess work for large $\tau$ in gapped systems. This result is at odds with the $1/\tau$ prediction of the geometric approach to optimization, which is predicated on the slow evolution of open systems close to canonical equilibrium. In contrast, our approach does not lead to an obvious geometric interpretation. Furthermore, as the thermodynamic work does not depend on how an isolated quantum system is split into a system of interest and its environment, our results imply the failure of the geometric approach prediction even for open systems. Additionally, we provide alternative optimization procedures, both for slowly-varying processes described by adiabatic perturbation theory and for weakly-varying processes described by linear response theory. Our findings are benchmarked and confirmed through the application to the driven transverse-field Ising chain.


I. INTRODUCTION
In the last decades, we have witnessed the development of several experimental techniques that enable the control and manipulation of few atoms, molecules or particles at the nanoscale. Using optical tweezers [1][2][3], ion traps [4], nuclear magnetic resonance [5], optical lattices [6] and other techniques [7,8], different experimental setups have been implemented to study several kinds of control and applications with potential development of new technologies.
The need to refine the manipulation of such small systems, increasing the level of control and decreasing the corresponding costs, has motivated intense theoretical activity on the topic. Among the prominent theoretical proposals are those of shortcuts to adiabaticity [9][10][11][12][13][14][15][16][17][18][19], which encompass many related methods. They all have the same goal, namely, the implementation of finite-time controls such that the final state of the system of interest is exactly identical to that produced by an equivalent quasistatic manipulation.
The remarkably good results of this theoretical framework do not hide, however, a few limitations. The implementations often require additional control mechanisms [9][10][11] or specific characteristics of the control fields [17,20], restricting the possible setups in the laboratory. Extensions to many-body systems may also face some difficulties, although it has been proved possible in some cases [20][21][22][23][24]. The different shortcuts to adiabaticity have been mostly restricted to isolated systems, i.e., systems whose dynamics is Hamiltonian. However, some extensions to open quantum systems have surfaced [25][26][27]. In addition, only recently the issue of quantifying * asorianialves@gmail.com † mbonanca@ifi.unicamp.br the cost to implement these control techniques has been addressed in the literature [28,29].
In parallel to this, the research on optimal finite-time thermodynamic processes has been seeking similar goals, although using a different perspective [30][31][32][33][34]. In this field, the minimization of the energetic cost during the finite-time control is the main target. In light of what was just presented, we might ask ourselves whether optimal finite-time control is also a kind of approximate shortcut to adiabaticity, and hence whether the two approaches can complement each other [35]. In the classical regime, the geometric approach to optimal slowly-varying processes [36,37] has been successfully applied to several interesting situations [38][39][40][41][42][43][44][45][46][47][48]. It has been also generalized to the quantum regime when the system of interest is weakly interacting with a heat bath [49][50][51][52][53]. In this approach, the energetic cost is expressed as an integral of a Lagrangian that is understood as a metric. The protocols minimizing such functional are then the corresponding geodesics.
In the present work, we show that the geometric approach cannot be applied to isolated quantum systems with a gap. We prove this rigorously through adiabatic perturbation theory (APT) [54]: for a slowly-varying process of duration τ , APT provides a systematic perturbative expansion in powers of τ −1 . The energetic cost is shown to decay as τ −2 , in clear contrast to the asymptotic τ −1 prediction of the geometric approach [36,37,49,50,52]. Additionally, we show that the latter cannot be applied to open quantum systems either, when the underlying microscopic dynamics of system plus environment is Hamiltonian, provided the full system is gapped. Our approach also allows for the optimization of slow processes. Taking the quantum Ising chain at zero temperature as a benchmark test, we compare our theoretical predictions with numerical simulations of the exact time-dependent dynamics, with excellent agreement. We also provide optimal finite-time protocols, following Ref. [55], in the complementary regime of fast but weak processes, where the variation of the external parameter is small compared to its initial value.
We will describe in detail the conditions in which our results apply. Foremost among them is the presence of a spectral gap. At this point, we just point out that gapped many-body systems are very common. Of course, in the important case of small systems, there are always finite-size gaps. However, even in the thermodynamic limit, gapped systems abound: (a) conventional s-wave BCS superconductors [56], (b) quantum Hall systems (both integer and fractional) [57], (c) integer-spin (Haldane) chains [58], and (d) quantum disordered paramagnets [59], to name a few. The quantum Ising model, whose one-dimensional version we study here, is just one example of a wide class of magnetic systems that are gapped except at a point in the phase diagram, namely, its critical point [see Ref. [59] for numerous examples]. In fact, gapless behavior often relies on additional conditions, as (a) the presence of exact symmetries, which give rise to Goldstone modes in phases with a spontaneously broken continuous symmetry or (b) some topological effect, as in half-integer spin chains [58] or gapless quantum spin liquids [60]. Departures from exact symmetries, as provided, e.g., by spin anisotropies induced by spin-orbit interaction, frequently end up generating a finite gap in these systems.
The performance of the two perturbative approaches applied here is tested far from equilibrium, and a clear qualitative difference is observed in the optimal protocols when the transition from one regime to the other is made. This highlights that different optimization strategies are required in different regions far from equilibrium. Only two parameters are necessary to delimit such regions, namely, the relative change of the control parameter λ (which is the transverse field in the case of the Ising chain) and the ratio between a characteristic time scale of the system and the duration τ of the process.
We show that the leading-order expressions for the energetic cost are quadratic forms in the speedλ in both perturbative approaches [55,61]. However, the notion of a metric, as it exists in the geometric approach, seems to be absent. In the regime of fast but weak processes, the reason for that is related to the memory kept along the process by the generalized force conjugate to λ. In other words, the average force at a given time t 0 depends not only on λ(t 0 ), but on the whole history of the protocol λ(t) up to t 0 . In the regime of slowly-varying processes, adiabatic perturbation theory provides a very simple quadratic form that depends only on the initial and final values ofλ. It is interesting to note that the metric of the geometric approach and the two quadratic forms we obtain here are derived from the same object, namely, an autocorrelation function of the generalized force.
The presentation is organized as follows. We briefly review the basics of the geometric approach to finite-time thermodynamics in Sec. II. The thermodynamic work in isolated systems is presented in Sec. III, and in Sec. IV we describe in detail the behavior that adiabatic perturbation theory predicts for the excess work in slow processes, highlighting the differences to the geometric approach and showing possible avenues of optimization. The complementary regime of weak processes is then discussed in Sec. V, together with the recently proposed method of optimization based on linear response theory. We briefly discuss the application of the presented methods to nonintegrable systems in Sec. VI and the failure of the τ −1 scale even for open systems in Sec. VII. A discussion of the results and the assumptions we made to obtain them is given in Sec. VIII and the paper is concluded in Sec. IX.

II. GEOMETRIC APPROACH
Consider a system with a collection of external parameters λ to be varied between times t i and t f , where τ = t f − t i is the process duration. A specific choice of the time dependence of λ(t) is called a protocol, and here we only consider protocols that can be written as a function of t/τ . In any process, the average work W delivered to the system can be separated in two parts, where W qs is the τ -independent (and protocol independent) quasistatic contribution, while W ex is the excess contribution, which embodies the extra energy we must give to the system in order to the carry out the process in finite time. The first part, W qs , only depends on the end points of the process and represents, therefore, an inescapable energetic cost; the second part, W ex , is protocol-dependent and can be minimized with a clever choice of how we drive the system. No matter how it is done, however, W ex → 0 as τ → ∞, the adiabatic (quasistatic) limit. The original formulation of the geometric approach to optimal driving [36,37] employs linear response theory (LRT) to describe the system's dissipation for long but finite process duration. It equates the excess work to where overdots denote time derivatives and the superscript T denotes the transpose of a matrix. The function ζ(λ) is the aptly named friction tensor, which is large for values of λ where the system dissipates more energy. Although Eq. (2) was initially obtained for classical systems, it should also be valid for quantum systems, as long as ζ is suitably defined from LRT's quantum formulation of the response function. Nevertheless, it is also possible to derive an equation with the same quadratic-in-λ form of Eq. (2) starting from a strictly quantum Lindblad master equation [49,50] and arriving at a different definition of ζ.
The friction has a clear geometric interpretation: it induces a Riemannian metric in parameter space, and the geodesics of such space (calculated with straight application of Euler-Lagrange equations) are paths of least resistance, i.e., protocols that minimize dissipation [36]. Additionally, Eq. (2) predicts that the excess work behaves asymptotically as τ −1 , as can be readily seen with a change of integration variable to s = (t − t i )/τ . Of course, this is consistent with the vanishing of the excess work in the adiabatic limit. Nevertheless, it predicts a regime in which the decay of W ex is universally given by τ −1 regardless of the shape of λ(t).
Both the LRT and the Lindblad derivation of Eq. (2) assume that the system of interest is in contact with a heat bath. They further assume that the slow evolution never takes the system too far away from canonical equilibrium. In what follows, we show that, for thermally insulated quantum systems with equilibrium initial states, the description of the excess work for slow processes is remarkably different. More precisely, we prove that, within the region of validity of our approach, the τ −1 behavior of the excess work for large τ does not exist in this context.

III. EXCESS WORK FOR ISOLATED SYSTEMS
Consider an isolated quantum system described by Hamiltonian H(λ) and evolving in time under unitary dynamics. Define the instantaneous eigen-equation H(λ)|n(λ) = E n (λ)|n(λ) . For simplicity, we focus on systems with discrete non-degenerate spectra (this assumption can be relaxed) and restrict ourselves to the case of a single external parameter λ, to be varied from λ i to λ f . Assuming that the system's initial density matrix is an equilibrium distribution, ρ(t i ) = n p n |n(λ i ) n(λ i )|, the density matrix at time t is where U is the unitary time evolution operator and |ψ n (t) = U (t, t i )|n(λ i ) is the solution to Schrödinger's equation with initial condition |ψ n (t i ) = |n(λ i ) . To have an initially canonical distribution, one would choose p n = Z −1 i e −βiEn(λi) , where β i is the initial inverse temperature and Z i is the associated partition function.
Defining the average work as the difference between average energies at the end and at the beginning of the process, W (τ ) = Tr where is a transition probability (or survival probability if m = n), i.e., the probability of finding the system in state |m(λ f ) at t f given that it was in state |n(λ i ) at t i . In turn, E m (λ f ) − E n (λ i ) is the energy difference of the transition. The average work, as written in Eq. (4), is equal to the average of the two-point measurement work, which is the quantum definition of work in a single realization of the process that obeys the quantum versions of Jarzynski's equality [62] and Crook's fluctuation theorem [63]. Note, however, that this is only the case because we assumed an initial equilibrium state. Alternatively, the same definition of average work used above can be put in an integrated power form, where we used the generalized force conjugate to λ, Equation (6) will be useful in the discussion of weak processes of Sec. V.
In the adiabatic limit, the solution to Schrödinger's equation is given by the adiabatic theorem [64] where we omit the time dependence of λ to evince the parametric evolution. In Eq. (8), the superscript zero signifies that we are in the adiabatic limit and φ n (t) contains the usual geometric and dynamic phases, which is notably process independent. The density matrix of Eq. (3), for a canonical initial distribution and in the adiabatic limit of Eq. (8), reads which is not a canonical distribution at time t, since the statistical weights are always evaluated at t i . Therefore, a generic isolated quantum system in a quasistatic process with an initial canonical distribution is taken away from canonical equilibrium at later times. This contrasts with the assumption of small deviations from canonical equilibrium taken in the geometric approach. In specific systems where E n (λ)/E n (λ i ) is independent of n (such as the harmonic oscillator with varying frequency or the ideal gas with varying volume), the density matrix keeps its canonical form with a time-dependent temperature, but this is not true in general. This serves as a hint that the geometric approach is unsuited for isolated systems. We know from Eq. (1) that the excess work is obtained by subtracting Eq. (10) from (4). Using the identity p n|n = 1 − m =n p m|n , we arrive at where E mn (λ) = E m (λ) − E n (λ) and the prime indicates that the diagonal term m = n is not included. It is noteworthy that the energy difference appearing in Eq. (12) is evaluated solely on λ f . This means that the parametric energy variation of the eigen-states is entirely accounted for in the quasistatic work of Eq. (10), which leaves only the energy variations from transitions (and not survivals) to be considered in the excess work. Not only does the expression in Eq. (12) vanish in the adiabatic limit, but the minimal work principle [65] ensures that it is always non-negative for non-degenerate systems with initial density matrices obeying p n ≥ p m if E n (λ i ) ≤ E m (λ i ), such as the canonical distribution or the density matrix of a pure ground state. Equation (12) contains the expression for the excess work we will consider in slow processes. Of course, exact evaluation of this equation requires solving Schrödinger's equation to determine the transition probabilities. In the next section, we introduce adiabatic perturbation theory and its approximate expressions for the excess work.

IV. EXCESS WORK IN SLOW PROCESSES
Adiabatic perturbation theory [54] provides perturbative corrections to the adiabatic theorem, valid when the process duration τ is large but finite. In practice, it gives where |ψ (r) is the rth order correction written in the instantaneous basis of H(λ). The expansions of Eqs. (13) and (14) are purposefully constructed to include the adiabatic limit in its r = 0 term, with C  mn (t) can be systematically calculated [54]. For example, the expression for r = 1 and m = n is where φ mn (t) = φ m (t) − φ n (t), and F mn (λ) = m(λ)|F (λ)|n(λ) are the instantaneous matrix elements of the force, defined in Eq (7). Since λ is a function of t/τ , it givesλ ∝ τ −1 , which sets the order of C (1) mn (t) in Eq. (15). Similarly, C (2) mn (t) containsλ andλ 2 , both of which are proportional to τ −2 , and this continues on for higher orders [54]. Thus, a superscript (r) means that the given quantity is proportional to τ −r . However, τ −1 , by itself, is not the proper quantity to determine the validity of the theory. To assess how accurate first-order APT is, it is more appropriate to consider the inequality which, looking at Eq. (15), is always true as long as an inequality known as quantitative adiabatic condition [66,67], a validity condition for the adiabatic theorem itself.
At this point, we describe in more detail the conditions under which we expect our approach to be valid. From Eq. (18), it is clear that APT breaks down if some relevant E mn vanishes. However, there are some important cases in which this will not happen. First, in finite systems, the energy differences E mn are bounded from below by finite-size gaps. Second, systems that are gapped in the thermodynamic limit (see several examples in Sec. I) will not violate Eq. (18) if they are driven at T = 0 or at temperatures much smaller than the gap. This can be seen from Eq. (12), where it is clear that transitions other than from the ground state are absent or strongly suppressed by the Boltzmann weight. This will be the case of all our numerical simulations later on.
In possession of Eq. (13), one can easily write down the transition probabilities appearing in Eq. (12). We have and, in practice, one calculates Eq. (19) to the desired order. We are now able to determine the corrections to the excess work coming from the corrections to the timedependent state. It must be pointed out that, since the sum in Eq. (19) is squared, a specific correction C (r) mn may emerge in many orders of correction to W ex , not only in order r.
We already know that the r = 0 term of the theory reproduces the quasistatic work, and thus does not contribute to the excess work. Following our superscript convention, truncating Eq. (19) at r = 1 leads to This is the central result of this paper. According to it, in a slow but finite time process of an isolated gapped system, there is no first order correction to the average work in the inverse process duration τ −1 . This result, which relies on the condition of initially diagonal density matrices, contrasts with the geometric approach prediction of W ex ∼ τ −1 for slow processes. We stress that an initial diagonal state describes an initial thermodynamic equilibrium between the internal parts of the isolated system. Furthermore, Eq. (20) relies neither on a specific number of degrees of freedom nor on a certain partition of the isolated system. Thus, it also applies to open quantum systems whose underlying microscopic dynamics with its environment is Hamiltonian and gapped, as further explained in Sec. VII. Continuing on for higher orders, we see that the r = 1 term of Eq. (19) is enough to describe the second order correction to the work, which is the first non-zero correction. Such correction has already been considered in the context of tradeoffs between power and efficiency of quantum heat engines [68,69]. Equation (21) is the starting point for the optimization of the energetic cost of slow processes in isolated quantum systems.

A. Optimizing the excess work
In contrast to Eq. (2), APT results for the excess work have no functional dependence: Eq. (21) only depends on λ andλ at t i and t f , not on the entire history of the protocol λ(t). In fact, ifλ(t i ) = 0 =λ(t f ), Eq. (21) also vanishes (see Eqs. (15) and (16)) and the first non-zero correction would be where C Emn(λ) . This continues on for higher orders: zeroing out the first r derivatives of λ at the beginning and at the end of the process makes the first non-zero correction to the work be W (2r+2) . This approach to reducing the order of the excess work is known as boundary cancellation method (BCM) [70][71][72][73][74]. There is a simple formula for generating a BCM protocol of order r, namely [74] λ which is a polynomial for any integer r > 0 with all time derivatives up to order r vanishing at t i and t f . However, there is a caveat: BCM presupposes the validity of APT. First, adiabaticity must be secured, and forcing time derivatives of the protocol to be zero at the end points might even be detrimental to reaching that goal. This is because the protocol will have to "speed up" during the process, which might spoil inequalities (17) and (18).
Conversely, we can abuse inequalities (17) and (18) to find protocols that adhere to APT as soon as possible. The fast quasi-adiabatic (FQA) strategy [75][76][77][78][79] amounts to setting the LHS of inequality (18) equal to an initially undetermined constant, and solving this first-order differential equation for λ with boundary conditions λ(t i ) = λ i and λ(t f ) = λ f . Since there is only one constant of integration, c must also be used to impose the boundary conditions, which ultimately leads to c ∝ τ −1 . Comparison of Eq. (24) to Eq. (15) leads to the conclusion that the FQA protocol distributes first-order transition probabilities uniformly throughout the entire duration of the process.
As it is, Eq. (24) only works when considering a single pair of energy levels m and n, but we can use the expressions for the work in Eqs. (10) and (21) to take every level pair into account. To this end, and inspired by Eq. (24), we set the ratio |W (2) /W (0) | equal to a tobe-determined constant c and disregard the terms that are not fully t-dependent, which leads tō Equation (25) should then be solved with λ(t i ) = λ i and λ(t f ) = λ f , a strategy which we also refer to as FQA for consistency. In contrast to BCM, FQA does not cancel Eq. (21), but it may guarantee the validity of such equation for the smallest τ possible. In essence, BCM is an answer to the problem "given APT validity, find protocols that minimize the cost order by order". Conversely, FQA is an answer to the problem "find protocols that ensure APT validity and are as fast possible", "fast" meaning small process duration. Both strategies have their advantages, as we shall see next.
In the following analysis, it will be useful to write the external parameter as where ∆ ≡ λ f −λ i is the total variation of λ in the process and g(t) is a function obeying g(t i ) = 0 and g(t f ) = 1. The discussion has been general so far, and any discrete and non-degenerate quantum system can be subjected to any of these strategies. However, to compare their efficiency, we apply them to a specific model: the transverse field Ising model (TI) [80], a one-dimensional chain of N spins, initially prepared in its ground state. The specifics of the system are detailed in appendix A. The external parameter is the magnetic field λ = B, to be compared with the fixed coupling J between neighboring spins. Here, for the sake of illustrating our results using this model as an exactly diagonalized testbed with a well defined thermodynamic limit, we consider only processes that do not cross the system's quantum critical point (at B = J), because it brings unwanted difficulties (that we have considered elsewhere [24,61,81]). The results for processes entirely contained in the paramagnetic phase (B > J) or in the ferromagnetic phase (B < J) can be seen in Fig. 1. Figure 1a shows four driving protocols for a paramagnetic process with B i = 10J. LIN is the naive linear protocol, which we will always use as a basic protocol we want to outperform. BC1 and BC2 are the two lowest order BCM polynomials, given in Eq. (23) with r = 1 and r = 2, respectively (r = 0 reproduces LIN). Lastly, FQA is the solution to Eq. (25) and, unlike the other three protocols, it depends on ∆ - Fig. 1a was generated with ∆ = 10J. Figures 1b and 1c show the comparison between numerically obtained exact results (symbols) and APT (lines) predictions for the excess work per spin. Figure 1b shows the excess work per spin vs τ and was constructed using the protocols of the previous plot with ∆ = 10J. From this figure, it is apparent that the secondorder APT result of Eq. (21) for the FQA protocol approximates the dynamics "early": already at Jτ = 10 −2 we see good agreement with the exact result, which is a consequence of FQA's design. It, however, has the standard τ −2 decay for large τ , since it does not necessarily slow down near t = t i and t f . The LIN protocol has the same decay and, in this case, its excess work is always greater than FQA's excess work. Conversely, among the considered protocols, BC1 and BC2 take the longest to be well approximated by APT, but, once they enter APT regime, they behave as τ −4 and τ −6 , respectively.  Fig. 1c is exactly the one shown in Fig. 1a, the same for every value of ∆ plotted -this was done for ease of representation. This plot shows how the FQA protocol is sensitive to the energy spectrum: it gives higher excess work than LIN when the field variation is small. On the other hand, BC2 already outperforms every other protocol shown for this not so large value of τ . Figures 1d-1f mirror the analysis of the previous plots, but for a ferromagnetic process with B i = 0.1J. Many of the features of the paramagnetic process are repeated in the ferromagnetic process. Worthy of note is the larger time decades considered in Fig. 1e (when compared to Fig. 1b) and an apparent coalescence in the excess work of Fig. 1f for every protocol considered, both of which are consequences of the proximity to the critical point. Indeed, on the far right of Fig. 1f, the last points in the plot represent processes in which the critical point is crossed, but this is of no consequence to our analysis since this region of the plot is far from the region of validity of the perturbation theory used.
Overall, we can see that the τ −1 behavior of the excess work is never observed as the leading order behavior for slow processes in Fig. 1.

V. EXCESS WORK IN WEAK PROCESSES
To complement the analysis of slowly-varying processes, in this section we study weak processes. To be specific, "weak" means that the total variation ∆ of Eq. (26) is small when compared to initial value of the external parameter, λ i . Such processes are well described by LRT [82], in which the Hamiltonian of the system is expanded to first order in ∆, where H i = H(λ i ) is the initial Hamiltonian and F i = F (λ i ) is the initial force of Eq. (7). Equation (27) is plainly written in standard timedependent perturbation theory form, whose first order result for the density matrix of Eq. (3) is [64] where |n i and E i n are initial eigen-states and eigenenergies. Then, the average force, itself expanded up to first order as where we defined the response function and Φ i (t) = Φ(λ i ; t).
Finally, the average work is obtained by placing Eq. (29) into the integrated power expression, Eq. (6). After an integration by parts of the last term, the boundary term and the previous two terms combined reproduce the weak limit of the quasistatic work in Eq. (10). Thus, the leftover term must be the lowest order contribution to the excess work, which we write as [55,83] The relaxation function Ψ relates to the response function as Φ(λ; t) = −∂ t Ψ(λ; t) and, using Eq. (30), we can write and Ψ i (t) = Ψ(λ i ; t). It is noteworthy that Eqs. (30) and (32) display the LRT functions Φ and Ψ in alternate forms, compared to the usual correlation functions (see Sec. VI). For isolated quantum systems, both forms are equivalent. Equation (31) is the staring point for the optimization of the energetic cost of weak processes in isolated quantum systems.

A. Optimizing the excess work
An effective strategy to minimize (31) consists of expandingġ(t) in some basis of functions [55]. To do that, we first rewrite Eq. (31) in a more symmetric form, where the property Ψ i (−t) = Ψ i (t) was used. Note the similarities and differences between Eq. (33) and Eq.
(2) of the geometric approach -both of them present quadratic functional dependence on the first derivative of the external parameter, but the kernel of Eq. (33), Ψ i (t − t ), shows explicit time dependence. The derivativesġ(t) andġ(t ) can now be properly expanded in the interval [t i , t f ]. Due to their convenient mathematical properties, the Chebyshev polynomials T l (x) are a good choice of efficient basis of functions [84]. Following Refs. [55,84], the truncated and regularized expansion ofġ(t) in a finite number n of polynomials where a l are the coefficients to be determined in the optimization and the factors g n,l regularize the truncated expansion [84] to avoid the Gibbs phenomenon at the extremities of the expansion interval. Their expression is [84] g n,l = n − l + 1 n + 1 cos πl n + 1 Substituting expression (34) into Eq. (33), we obtain where the A lj are given by The excess work (36) becomes then a finite multidimensional quadratic form in the a n whose minimum we want to find. The boundary conditions g(t i ) = 0 and g(t f ) = 1 are additional constraints in our optimization problem. Using the method of Lagrange multipliers, we can find then the coefficients a l that provide the optimal protocol by solving a set of linear algebraic equations. Naturally, we call this optimization strategy Chebyshev expansion up to n modes (CEn). We remark that the relaxation function Ψ i (t) is the main physical input in this procedure. The factors A lj crucially depend on the switching time τ and Ψ i (t). Due to its relation with the response function, Ψ i (t) can be obtained from experiments when it is not accessible theoretically.
To test this optimization strategy with the relaxation function of an isolated system, we once again employ the TI chain. Its relaxation function Ψ i (t) is given in Eq. (A7) of appendix A when the chain is initially prepared in its ground state. Deep in the paramagnetic phase (B i J) and in the thermodynamic limit N → ∞, we can swap sum by integral and approximate the quantities of Eq. (A7) as This leads to where J n is the Bessel function of the first kind (the caligraphic J is used to avoid confusion with the coupling constant J).
The function of Eq. (38) has a clear decay factor in J 1 (2Jt)/Jt, indicating a relaxation time scale of J −1 ; and a pure oscillating factor in cos(2B i t), indicating an oscillation time scale of B −1 i (see Fig. 2). From Eq. (A2), it can be seen that the relevant energy gap of the system in the paramagnetic phase is B i . Thus, in the TI chain, it is the oscillation time, not the relaxation time, that determines its closeness to adiabaticity. The results using the relaxation function of Eq. (38) can be seen in Fig. 3. The symbols in the plots for the excess work per spin depict the numerical results obtained from the integration of the exact time-dependent dynamics, while the lines represent the LRT prediction. Figure 3a shows the linear protocol and three protocols obtained from the Chebyshev expansion for Jτ = 2. Somewhat surprisingly, these protocols look identical to the BCM protocols we presented in the context of slow processes. Indeed, a glance at Fig.1b reveals that for this process duration, all of the protocols considered in that plot already agree with the APT prediction. This indicates that Jτ = 2 can be considered slow, and the minimization of the LRT functional exploits this fact to generate protocols that follow the BCM guideline. However, note that for realistic values of ∆, i.e., 10 −1 ≤ ∆/B i ≤ 10 1 , increasing the number of modes worsens the performance, as can be seen on Fig. 3b. The work only decreases with increasing number of modes for much smaller values of ∆ (see Fig. 3c), closer to the point where LRT starts agreeing with the numerical data. In any case, all three CE protocols outperform the linear one for the entire ∆ range of Fig. 3c.
The situation is considerably different for Jτ = 1. Figure 3d compares the linear protocol once again with three Chebyshev expansions, with the same number of modes considered in the previous case. In this case, we start seeing deviations from BCM protocols, but only when using 15 modes or more. Essentially, the LRT functional cannot output better protocols than BCM ones when given a small number of degrees of freedom to optimize. When it has enough degrees to work with, it generates oscillating protocols that outperform BCM-like protocols while escaping APT description. Performance-wise, Fig. 3e shows that increasing the number of modes does seem to decrease the excess work around ∆/B i = 10 −1 , but not beyond this point. In fact, the oscillating protocol CE17 is outclassed even by the linear protocol for ∆/B i values as low as 10 0 , demonstrating that such oscillating protocols are only reliable in the LRT regime ∆/B i 1. For Jτ = 1, increasing the number of modes past 17 starts giving non-realistic oscillatory protocols, with amplitudes several orders of magnitude higher than the endpoint. The maximum number of modes that still gives realistic protocols in fact decreases with decreasing τ . Figure 3g shows LIN and three CE protocols, the highest number of modes used being 11, for Jτ = 0.5. CE11 already shows considerable oscillations -the protocols for higher number of modes are impractical. Once again, the performance of these oscillatory protocols is only good for small values of ∆: Fig. 3h reveals that, when compared to CE9, CE11 has a marginal advantage at best and a large disadvantage at worst.
The trend of the Jτ = 0.5 (Fig. 3g) and Jτ = 1 (Fig. 3d) cases studied seems to indicate that, using a high enough number of modes, oscillating protocols can be generated from the LRT functional for any value of τ . However, we have not found oscillating protocols for Jτ = 2 (Fig. 3a) for a number of modes of up to 23, where numerical errors started becoming a factor. Figure 4 shows the excess work per spin (obtained numerically) vs τ for the protocols shown in Fig. 3d and CE11 shown in Fig. 5a. These plots demonstrate how specific are the protocols obtained from the LRT functional for the values of τ used to generate them. consequence of their BCM-like appearance. As ∆ is increased, the valleys become dispersed and less profound, as can be seen in Fig. 4b for ∆/B i = 1. When the system is far from LRT description, as for ∆/B i = 10 in Fig. 4c, the valleys still exist, but the protocols CE15 and CE17 lose any advantage they had -another sign that the oscillatory CE protocols only guarantee a decrease in the excess work when one stays in the weak regime, as mentioned before.
Note that, in accordance to what we discussed in the context of slow processes, for every protocol considered in Fig. 4, the excess work decays as τ −2 for large τ . This includes CE11 and CE13, which one would naively predict to have steeper decays, as expected from BCM protocols. The reason for this is the fact that CE11, for instance, is not exactly a BCM protocol -after all, unlike Eq. (23), it has to accommodate for the specific properties of the system, imported from the relaxation function. Figure 5 compares CE11 with BC3, obtained from Eq. (23) with r = 3. They are almost indistinguishable from the point of view of Fig. 5a, but Fig. 5b reveals that CE11 has finite first derivatives at the endpoints, which means APT gives W ex ∼ τ −2 . Conversely, BC3 has a considerably steeper decay of τ −8 in the APT regime, as can be seen in Fig. 5c. Nevertheless, such small derivatives are the reason why CE11 consistently beats LIN in Fig. 4.
The absence of the universal τ −1 decay for large τ is hence corroborated using another perturbation theory. We emphasize that the analytical results of this section considered the thermodynamic limit of the Ising chain through expressions such as (38). Hence, the τ −1 scaling does not seem to be related to this asymptotic limit. The agreement between both perturbation theories used in this paper in the appropriate weak and slow regime is elucidated in appendix B.
In passing, we remark that the expression for the relaxation function of Eq. (A7), deep in the ferromagnetic phase (B J) and in the thermodynamic limit, is ap-proximated by Disregarding the constant prefactor, Eq. (39) is very similar to Eq. (38), only differing by the substitution B i ↔ J. This is related to the Kramers-Wannier duality between the two phases of the system [85]. Consequently, the LRT results for the ferromagnetic phase are the same as those for the paramagnetic phase, given suitable values of B i and τ . For instance, the protocols obtained from the LRT functional for B i = 0.1J and Jτ = 10 are identical to the protocols obtained for B i = 10J and Jτ = 1, given in Fig. 3d.

VI. NONINTEGRABLE SYSTEMS
The optimization strategies presented and discussed in this paper are naturally written in terms of the eigenquantities of the Hamiltonian. Yet, unlike the transversefield Ising Hamiltonian of Eq. (A1), the vast majority of Hamiltonians cannot be diagonalized exactly, i.e., most systems are nonintegrable. When the system is numerically integrable, the strategies contained herein are readily available. For the application of the Chebyshev expansion method, the eigen-equation need only be numerically solved at λ i . However, the eigen-equation must be solved for every λ value traversed in the defined process for the application of the FQA method, since it requires knowledge of how the spectrum changes with λ.
The optimization strategies contained herein can also be applied in experimental situations where the linear response of the system can be determined through measurements. To elaborate, the methods can be used if the generalized force F -which is λ-independent if the external parameter is linearly coupled to the system -is known and if the response function can be determined from the correlation function where is the (equilibrium) density matrix in the adiabatic limit, [·, ·] denotes the commutator and is the generalized force in the instantaneous interaction picture of H(λ).
Note that the trace in Eq. (40) is an equilibrium average, since the state ρ (0) is diagonal. The success of linear response theory is largely based on the measurability of this average. For instance, with an initially canonical distribution, one needs to determine Φ i (t) = Φ(λ i ; t) and calculate Ψ i (t) to apply the Chebyshev expansion method. On the other hand, to determine fast quasiadiabatic protocols, one first rewrites Eq. (25) as where Υ(λ; t) is a function described in appendix B, while the right-hand side contains the average energy in the adiabatic limit. Then one needs to determine the average in Eq. (40) for every λ in the process, which might not be an easy task since the state of the system may not maintain canonical form throughout the evolution, as discussed in Sec. III. If this can be done, one then calculates Υ from Φ(λ; t) = −∂ 3 t Υ(λ; t) and solves the differential equation (43).

VII. OPEN SYSTEMS
In this section, we briefly consider the implications of our analysis to open quantum systems. We stress that at no point during our derivation of the excess work in slow processes (Sec. IV) did we assume that the external parameter λ interacts with the entire isolated system. Thus, in a gapped bipartite system (call it U ), if λ acts only on one of the parts (call it S), the excess work exerted on U still decays asymptotically as (at least) τ −2 for slowly-varying processes (we also assume the total energy of U is well below its gap). However, since the force applied on S must also be the total force applied on U (there are no other external influences), it follows that the work exerted on the open system S is equal to the work exerted on the closed system U . Consequently, through a Hamiltonian description of a closed and gapped quantum system U that includes the open quantum system S, we can see that the excess work done on open systems also scales as τ −2 in slow processes.
The above result can be made precise with the operational definitions of average work used in this paper. The Hamiltonian of the isolated system U , whose control parameter couples only with the part S, can always be written as where H S (λ) is the Hamiltonian of the open system we manipulate (thus dependent on λ), H R is the Hamiltonian of the rest R (thus independent of λ) and H I is the Hamiltonian of the interaction between S and R (also independent of λ). Now, using Eq. (6) to calculate the work exerted on U , we see that where we used Eq. (44) to write ∂ λ H U (λ) = ∂ λ H S (λ), the fact that the trace over the entire system can be  6. Excess work vs process duration for a 3-spin Ising chain calculated numerically from the exact dynamics. The magnetic field applied on one of the spins was varied linearly in time from 2J to 3J, while the magnetic fields in the other two spins were kept constant at 3J. S+R represents the excess work done on the entire chain (calculated as the difference between average energies at the end and at the beginning), S represents the excess work done only on the manipulated spin (calculated as the integral of the force applied on that spin) and APT represents the first-order adiabatic perturbation theory prediction (stripped of its oscillations). The main plot was constructed using the ground state of the total Hamiltonian (46) at ti as the initial state, thus leading to Wex ∼ τ −2 . Meanwhile, the inset was constructed using a tensor product state: the ground state of the manipulated spin times the ground state of the other two spins combined, thus leading to Wex ∼ τ −1 .
decomposed into partial traces (Tr U = Tr S Tr R ) and ρ S (t) = Tr R ρ U (t) is the reduced density matrix of S.
In the last line of Eq. (45), we used Eq. (6) once again to identify the average work exerted on S. Therefore, the general conclusions we obtained for the work in closed gapped quantum systems also applies to open quantum systems (assuming system plus environment is also gapped), including the absence of the asymptotic τ −1 behavior of the excess work in slowly varying processes.
We emphasize that the above arguments are rather general. In particular, we make no assumption about the size of the system, nor do we assume specific features of the interaction between the open system and its environment. However, we do assume that the initial state of the isolated system U commutes with the initial Hamiltonian (44), a choice that exactly describes an initial thermodynamic equilibrium between S and R.

A. Example: small Ising chain
To illustrate our finding, we have included Fig. 6, which shows W ex vs τ when we manipulate only one spin of an Ising chain with three spins, with Hamiltonian (46) The magnetic fields B 1 and B 2 are constant, while the magnetic field B 3 (t) of the manipulated spin is varied linearly in time. For the main plot, we chose the ground state of Eq. (46) as the initial state. Hence the excess work on the entire chain decays as τ −2 for large Jτ . The excess work on the manipulated spin behaves exactly the same (apart from numerical errors). On the other hand, in the inset, the initial state was chosen to be the tensor product of the ground state of spin 3 and the ground state of spins 1 and 2 together, which is still a pure state but noncomutable with the initial Hamiltonian. In this case, the excess work decays as τ −1 for large Jτ , which agrees with the APT prediction for nondiagonal initial states. The notable oscillations are natural in such a small chain and they are predicted by APT, but we removed them from the APT result in this figure for ease of representation (see Ref. [61]).

B. Example: driven quantum Brownian motion
As a second illustrative example, we consider quantum Brownian motion in a harmonic trap whose minimum location varies in time. The Hamiltonian we use to model this phenomenon is given by where P , Q and M denote, respectively, the momentum and position operators of the Brownian particle and its mass whereas κ 0 denotes the stiffness of the trap. By λ(t), we denote the time-varying position of the minimum and H CL stands for the Caldeira-Leggett model [86], which describes a collection of N harmonic oscillators with mass m k , frequency ω k and position and momentum operators q k and p k interacting with the Brownian particle. According to Eq. (45), the work performed reads where the non-equilibrium average Q(t) , (see App. C) was obtained using the equilibrium initial state where Z i = Tr[exp (−βH BM (λ(t i )))] and β is the inverse temperature.
Inserting expression (50) in Eq. (49) and performing an integration by parts, it is straightforward to show that where Φ BM (t) = −dΨ BM dt. Hence, the excess work is exactly given by the LRT functional (31) derived in Sec. V.
To investigate the large τ behavior of Eq.(52), one can first take the limit N → ∞. In this case, an explicit expression for the so-called spectral density J(ω) has to be chosen in the continuum limit. Using the standard gapless Ohmic spectral density, with a high-frequency cutoff, one obtains an exponentially decaying Ψ BM (t) (see App. C for the details), where ω 2 0 = κ 0 /M , Ω 2 = ω 2 0 − γ 2 and γ = η/M . Plugging Eq. (54) in Eq. (52), the excess work W BM ex for the linear protocol, λ(t) = λ 0 + δλ (t/τ ), reads, where the leading order is τ −1 .
On the other hand, if N is large but kept finite, Ψ BM (t) will be a quasi-periodic function as expressed in Eq. (32). In this case, the excess work (52) for the same linear protocol will be given by a sum of terms proportional to where ω mn = E mn /h are the frequencies related to the spectral gaps of the isolated system, as defined before in Sec. III. The previous expression shows that the τ −1 term is absent and suggests that the limits of large τ and large N might not commute. This is particularly relevant for numerical simulations where N is always finite.

VIII. DISCUSSION
In this section, we summarize the assumptions we made and the assumptions we did not make to arrive at the asymptotic behavior of the excess work for large process duration, namely, W ex ∼ τ −2 for a generic protocol.
It is clear from Sec. IV that the important ingredients are Hamiltonian dynamics and isolated gapped systems at temperatures well below their gap or open systems treated as part of a larger closed gapped system (where APT can be applied). Furthermore, we also assume initially diagonal density matrices. This second condition reflects initial thermodynamic equilibriumρ(t i ) = 0, as per von Neumann's equation. Note that tensor products of diagonal states of open system and environment do not lead to diagonal states of the total composite system, given the presence of interactions between the parts. Although such product states are often cited as equilibrium states when the open system is weakly coupled to the environment (as in the context of the geometric approach [36,37,49,50]), the mild effects of the interaction can accumulate over long process durations and thus should not be ignored. Whenever the initial state is not diagonal, the excess work does scale asymptotically as τ −1 , as seen in the inset of Fig. 6.
As mentioned previously, our arguments are independent of the size of the system. Moreover, we assumed no relaxation mechanism to be present, such as observed in systems obeying the eigenstate thermalization hypothesis or in the chaotic regime [87,88]. This means our approach describes the asymptotic behavior of the excess work in any gapped system initially in thermodynamic equilibrium. Of course, that does not rule out the possibility of other excess work scalings happening before the asymptotic behavior takes place -indeed, two of us verified this in Ref. [61]. In any case, with or without relaxation, the excess work eventually behaves as τ −2 for τ large enough in gapped systems.
A connection between our analysis and relaxation mechanisms can be made through Eq. (33) written in the following form where the integration variables were changed to s = t/τ and s = t /τ , and the derivatives are nowġ(s) = dg/ds andġ(s ) = dg/ds . It is important to stress that Ψ i (t) can either describe few-body or many-body correlations depending on the observable that couples to the control parameter. Equation (57) suggests that the excess work scaling for weak and slowly-varying processes (the large τ limit) is ruled by the long time tail of Ψ i (t) (or, equivalently, by the small frequency behavior of its spectrum), which in turn is related to the long time tail of the response function (see Sec. V). It was verified in Sec. VII B that, using Eq. (57) for linear protocol, a simple exponential decay leads to an asymptotic τ −1 scaling. This type of response decay can be obtained in general from Lindblad master equations. However, we should keep in mind that large times predictions of this effective description eventually fail (see Ref. [89] for a recent discussion).
Although the TI chain, the system we used to corroborate our findings, is integrable and consequently nonrelaxing, the metric of the geometric approach can still be suitably defined for it. For example, following the original LRT derivation of the geometric approach, the metric in the paramagnetic phase of the chain can be written as [36,37] where Ψ(B; t) is the relaxation function of Eq. (38) with B i → B, Ψ(B, 0) = N J 2 /4B 3 and the so called relaxation time is This integral vanishes when the full paramagnetic relaxation function (38) is used, owning to the presence of the oscillatory cosine factor. Heuristically, we might argue that these oscillations should not contribute to a measure of relaxation of the system, and hence we can drop the cosine factor in this calculation. Considering just the Bessel envelope of Fig. 2, we arrive at τ R = J −1 , which is the natural decay time scale of the Bessel factor. Nevertheless, simply putting Eq. (58) into Eq. (2) does not reproduce the correct excess work in any of the numerical simulations shown in this paper.

IX. CONCLUSION
Despite increased recent interest in the development of a unifying theoretical framework to investigate the minimization of energetic costs in finite-time thermodynamic processes of quantum systems [52], it seems that such unification still evades us. Here, we have provided extensive evidence to support the claim that the celebrated geometric approach to optimal driving, in its current formulation, does not apply to quantum gapped systems evolving under unitary dynamics. As far as adiabatic perturbation theory and linear response theory reach, there is no sign of a thermodynamic metric from which geodesics can be obtained, neither is there the predicted τ −1 asymptotic behavior of the excess work. However, our analysis does not exclude the possibility of studying the so-called thermodynamic length along the lines of Refs. [90,91]. Attempts at a geometric approach for isolated quantum systems can be found in Refs. [92,93].
Although some elements of our numerical analysis of the transverse-field Ising chain may be system-specific, we expect most prominent features highlighted in this paper to be universal. In the context of slowly-varying processes, we have discussed and exemplified the minimum exponent decay τ −2 of the excess work, as guaranteed by adiabatic perturbation theory. We then tested two well-known optimization strategies, enabling us to indefinitely increase the exponent or ensure the minimum exponent decay for reasonably small process duration τ .
In the complementary regime of fast but weak processes, we applied a Chebyshev expansion approach to optimization, which is based only on the knowledge of equilibrium correlation functions of the system, in typical linear response theory fashion. The protocols obtained are surprisingly sensitive to the time scales present in the relaxation function, and they reproduce the boundary cancellation method in the slow limit.
In both regimes studied, the perturbation theories do not give lower bounds for the excess work of the protocols generated: higher order BCM protocols can always be used to guarantee steeper decay in the APT regime, just as higher order Chebyshev expansions seem to guarantee ever-increasing performance in the LRT regime. Overall, the present results seem to suggest that, in isolated quantum systems, there is no single way of obtaining optimal protocols for every scenario and, as a matter of fact, it is not clear if the notion of a unique optimization procedure exists.
Finally, we have briefly analyzed the consequences of our findings to open quantum systems with two simple but illustrative examples, showing that the asymptotic τ −1 scaling is indeed absent for initially diagonal states and that the limits of large N and large τ might not commute. More detailed work is however necessary to understand how the geometric approach predictions can be reconciled with the results presented here. This further analysis might focus on the role of the initial state of the full isolated system and on the reliability of the long time predictions of effective descriptions such as Lindblad master equations. With little effort to go beyond what we have shown here for initially diagonal states, it can be shown that an initial state that does not commute with the full initial Hamiltonian leads to a non-zero first-order APT correction to the excess work (as demonstrated in the inset of Fig. 6). Concerning Lindblad master equations, they often imply exponential decays of correlation functions which, according to our LRT approach, would lead to an asymptotic τ −1 scaling in the large τ limit.
Our analyses certainly leave several important open questions concerning the conditions required for the τ −1 scaling of the excess work. In this sense, this paper consists of a critical assessment of the claims that have been made so far within the geometric approach in quantum systems. In particular, predictions coming from effective descriptions of open quantum systems often employ phenomenological assumptions that should be very critically analyzed, especially those that can impact the scaling discussed here. It remains to be proven then under which conditions truly Hamiltonian quantum dynamics might lead to the τ −1 scaling. so, we first place Eqs. (15) and (16) into Eq. (21), Now, since ∆ is small, λ f ≈ λ i and we can replace every instance of the former with the latter (including in the energy gap appearing in the adiabatic phases, inside the integral of Eq. (9)). Doing this carefully leads to (witḣ λ = ∆ġ) Lastly, expanding the squared absolute value in Eq. (B3) reveals that this expression is indeed equal to Eq. (B1). Both LRT and APT predict an excess work that is quadratic in their small perturbative parameter (∆ and τ −1 , respectively). Their agreement in the simultaneously weak and slow regime -which can be shown to hold for the microscopic state of the system, instead of the excess work -is a consequence of the fact that both are well defined perturbation theories for the same equation, namely, Schrödinger's equation.