Shortcuts to Adiabaticity in Krylov Space

Shortcuts to adiabaticity provide fast protocols for quantum state preparation in which the use of auxiliary counterdiabatic controls circumvents the requirement of slow driving in adiabatic strategies. While their development is well established in simple systems, their engineering and implementation are challenging in many-body quantum systems with many degrees of freedom. We show that the equation for the counterdiabatic term, equivalently the adiabatic gauge potential, is solved by introducing a Krylov basis. The Krylov basis spans the minimal operator subspace in which the dynamics unfolds and provides an efficient way to construct the counterdiabatic term. We apply our strategy to paradigmatic single- and many-particle models. The properties of the counterdiabatic term are reflected in the Lanczos coefficients obtained in the course of the construction of the Krylov basis by an algorithmic method. We examine how the expansion in the Krylov basis incorporates many-body interactions in the counterdiabatic term.


I. INTRODUCTION
In noisy quantum devices, dominant in the noisy intermediate-scale quantum (NISQ) era [1], the prospects of implementing exact adiabatic control protocols are dim.Noise generally lowers the fidelity of preparing a target quantum state, making the dynamics not unitary, and leading to a final mixed state.The presence of noise further limits the admissible operation time in adiabatic protocols, e.g., in adiabatic quantum computing and quantum annealing.In these devices, noise can act as a heating source leading to excitation formation [2,3], precluding the goal of finding the low-energy configuration of a given problem Hamiltonian.
The ubiquitous presence of noise in current NISQ devices forces us to rethink the use of adiabatic strategies.A natural approach is operating in timescales where environmental noise is negligible.A demonstration of this approach has recently been reported in quantum annealing devices, where noise-induced errors generated for moderate operation times [4] can be eliminated by shortening the duration of the process [5].However, this strategy generally limits the efficiency of the computation as a result of the adiabatic theorem, whether one considers the system closed [6] or open [7].An alternative approach relies on optimally tailoring the time dependence of the parameters that are varied in time in the system of interest (e.g., the harmonic frequency in a trapped system or a magnetic field in a spin system) [8].This is the principle behind the so-called boundary cancellation method that reduces excitations by devising smooth protocols in view of the adiabatic theorem in either isolated or open systems [7,9,10].Such an approach requires no additional control fields, easing the implementation of the driving protocols in the laboratory [11], but provides limited ad-vantages in the speedup, and technical assumptions in the adiabatic theorem may restrict its applicability.
While various techniques have been developed to engineer STAs, counterdiabatic (CD) driving stands out among them by providing a universal approach for any system in isolation.The early formulation due to Demirplak and Rice [28][29][30], independently developed by Berry [31], assumes the dynamics to be unitary and the system Hamiltonian to be diagonalizable at all times.However, progress over the past decade has shown that STAs can be applied to open quantum systems [32][33][34][35][36], as demonstrated in a pioneering experiment [37].
From the outset, the need for the system Hamiltonian to be diagonalizable at all times precludes the application of CD driving in important scenarios where such knowledge is unavailable, e.g., in quantum annealing.However, the development of approximate methods to engineer CD controls has challenged and de facto removed this requirement.Specifically, the early proposal of using digital methods for quantum simulation to realize CD controls [38,39], in combination with variational methods [40][41][42], has led to a framework for digitized-CD quantum driving for quantum algorithms [43,44], that include the use of STA in adiabatic quantum computation [45] and quantum optimization [46][47][48].
The nature of the CD controls remains currently an issue.In a system with many degrees of freedom, finding an efficient prescription to determine the CD fields is generally challenging.The first works exploring STA by CD driving in many-body quantum systems showed that the CD controls generally involved many-body interactions of arbitrary rank (one-body, two-body, etc.) [38,39,49,50].In addition, CD terms are generally spatially nonlocal [39].In systems of continuous variables, such as a harmonic oscillator or ultracold gases, CD terms cannot always be realized by applying an external potential [12,51] but may involve nonconservative momentum-dependent Hamiltonian terms [52][53][54].Likewise, in spin systems, CD terms may involve interactions among distant spins [38,39,49,50].As a result, one of the pressing problems in the development of STA is to find systematic approaches to tailor CD terms.One option is to find unitarily equivalent Hamiltonians for which STAs can be implemented exactly with experimentally available resources [53][54][55].Another relies on approximate protocols, determined through variational methods [39][40][41][56][57][58][59] or otherwise [49,60,61].
Current general approaches to engineering STA by CD driving are blind to any structure or symmetry in the actual dynamics.However, it is known that the presence of dynamical symmetries in a given process can significantly simplify the CD protocols required to control it and render the implementation of STA experimentally realizable.In cases where a dynamical symmetry is known, one can identify the CD controls in terms of the elements of a closed Lie algebra [62].However, the application of this approach has been limited to the restricted set of examples in which dynamical symmetries are known, i.e., few-level systems [63] and scale-invariant processes [64].
Further progress calls for novel approaches that systematically unravel and exploit any structure in the dynamics of the process to be controlled.This work introduces an approach that achieves this goal by formulating CD driving in Krylov space.Krylov subspace methods have a long tradition in numerical recipes and can be efficiently implemented using the Lanczos algorithm and its variants [65].In time-dependent quantum mechanics, Krylov space describes the minimal subspace in which the dynamics unfolds, greatly easing the computational resources to describe time evolution [66].They are further useful in foundations of quantum physics to characterize operator growth [67][68][69][70] and the fundamental speed limits governing it [71,72].Consider the case of the quantum dynamics in the Heisenberg representation.Given an observable of interest O 0 and a generator of evolution H, the evolution of the observable is set by where U (t) is the time-evolution operator.For time-independent Hamiltonians such evolution admits the expansion The dynamics generates the set of operators {L n O 0 } ∞ n=0 that are not orthonormal and generally live in an operator subspace known as the Krylov space.The Lanczos algorithm can be used to construct a basis in Krylov space and further provides the Lanczos coefficients that determine the entries of the matrix representation of the Liouvillian in the Krylov basis.

II. OUTLINE
In this work, we introduce a formulation of CD driving in Krylov space using the celebrated Lanczos algorithm.In Sec.III, we briefly review the key concept of STA.In the CD driving, for a given time-dependent Hamiltonian, the dynamics is assisted by an auxiliary control field known as the CD term.The CD Hamiltonian acts as the generator of adiabatic continuation, discussed in proofs of the adiabatic theorem, e.g., by Kato [6] and Avron and Elgart [73].Similarly, it has been discussed in the context of quasiadiabatic continuation by Hastings [74][75][76][77][78][79][80][81].Recent literature refers to the CD term as the adiabatic gauge potential (AGP).It is further related to the Berry connection, and its norm gives the real part of the quantum geometric tensor, i.e., the quantum metric tensor or fidelity susceptibility, as discussed in Refs.[30,38,82,83].Thus, the AGP has broad applications beyond quantum control, extending to quantum state distinguishability, quantum state geometry, adiabatic theorems, critical phenomena, quantum thermodynamics, etc.
Finding the explicit form of the AGP is a fundamental problem for practical applications and has been discussed from various viewpoints.The spectral representation obtained in the original studies [28][29][30][31] has a disadvantage as it is generally difficult to obtain the corresponding operator form in systems with many degrees of freedom.In that case, we can start the analysis from the operator equation for the AGP.The equation is solved approximately using the variational method [40].
In the variational method, the validity of the approximation strongly depends on the chosen ansatz.The operator equation is recast into an integral representation.It gives a nested commutator expansion and indicates possible operator forms of the AGP.Combining the nested commutator expansion with the variational method offers a systematic method for finding the AGP for complex systems [41].In most applications, the expansion is truncated to obtain an approximate result.It is implied that taking into the infinite-order expansion gives the exact result, although rigorous proof has not been shown.In the present work, we circumvent the need for the variational method by providing an exact closed-form expression for the AGP in the Krylov method.
In Sec.IV, we introduce the basic concept of the Krylov subspace method and develop a general framework for finding the exact AGP from the Krylov expansion.The Krylov expansion is formulated by defining a proper inner product and a Liouvillian superoperator for a target system.For a given initial seed operator, the Krylov subspace where the dynamics unfolds is determined from the Krylov algorithm.We show that a specific choice of the seed operator is useful to solve the equation for the AGP.The AGP is expressed by the Krylov basis and the Lanczos coefficients obtained from the Krylov algorithm.We find that the AGP is classified into two categories.They are characterized by the parity of the number of the Krylov basis.
Comparing the exact form of the AGP with the integral representation with the nested commutator expansion gives a close relation of the AGP to the complexity of the Krylov space.We discuss that the properties of the AGP can be understood directly from the series of the Lanczos coefficients and the operator wave functions defined from the general framework of the Krylov method.We also discuss how the variational method with the nested commutator expansion is justified.
In Sec.V, we apply the general framework to various canonical examples, including two-and three-level systems, and the harmonic oscillator.To be instructive, we demonstrate those well-known examples by using several different ways to determine the AGP.
The full potential of the present framework is displayed when it is applied to systems with many degrees of freedom.In Sec.VI, we treat integrable, nonintegrable, and disordered quantum spin chains.We first apply the method to a one-dimensional transverse Ising model without a longitudinal magnetic field.The AGP of the system is well known in that case [38,49,50] and we rewrite the result with respect to the Lanczos coefficients.We find that the quantum phase transition can be identified from the Lanczos coefficient series.When we apply the longitudinal magnetic field, the system becomes nonintegrable and the exact solution is not available.We consider a truncation of the Krylov expansion, and the result is shown to be equivalent to that of the variational method.In reporting explicit expressions for the AGP in many-body systems, our work advances the study of STA beyond the large body of literature focused on leadingorder truncations of the CD term [40,41,57,59].
We also discuss in the same section the onedimensional isotropic XY model.We treat several cases where the interaction couplings are uniform or random.Although the model can be mapped onto a free fermion model, the explicit construction of the AGP for a given set of coupling constants is a difficult task.We can formulate the expansion systematically and demonstrate the expansion up to a considerably large system size.We discuss closely how each order of the expansion affects the result.We also consider the case where the system is equivalent to the integrable system described by the Toda equations.We discuss the implications of the integrability condition on the Krylov expansion.
The present study is concluded with final remarks in Sec.VII.

III. ADIABATIC GAUGE POTENTIAL AND COUNTERDIABATIC DRIVING
Consider a closed quantum system described by the Hamiltonian operator H(λ) depending on the set of parameters λ = (λ 1 , λ 2 , . . .).Throughout this paper, a capital letter denotes an operator or a matrix.Let |n(λ)⟩ represent an eigenstate of the Hamiltonian with the eigenvalue ϵ n (λ).The time-independent Schrödinger equation, and the equation for adiabatic continuation read, respectively, (2) The phase of the eigenstate is fixed by requiring the relation ⟨n(λ)|∂ λ n(λ)⟩ = 0.The AGP operator A = (A 1 , A 2 , . . . ) is introduced by differentiating the eigenstate with respect to λ and enforces adiabatic continuation for all eigenstates; i.e., it is independent of n.
One of the prominent applications of the AGP is the CD driving [28][29][30][31]84].For time-varying parameters λ(t), we consider the time evolution Here, the CD term is introduced as , where the overdot denotes the time derivative.It prevents the nonadiabatic transitions among eigenstates |n(λ)⟩, which means that the solution of the Schrödinger equation ( 3) is exactly given by the adiabatic state of H: Even when the implementation of the CD term is challenging, we can use the CD term to assess the nonadiabatic effects [85].
While it is a nontrivial problem to obtain the explicit form of the AGP for a given Hamiltonian, its matrix elements can be formally written in terms of the spectral properties of H as The main aim of the present work is to find a systematic way to obtain the operator form of the AGP.The AGP satisfies where . This relation was used to find an approximate AGP by variational methods [40,41,86].We exploit this relation to obtain the exact form of the AGP.
As an alternative useful relation, one can invoke the integral representation introduced by Hastings in the context of quasiadiabatic continuation [74,75,77]: The integrand is proportional to the operator ∂ λ H(λ) conjugated by a unitary.Using the unitary operation is represented by L λ , we can perform the integration over s to write This formal expression motivates us to use the expansion [41,87,88] The construction of the AGP in Krylov space that we present in the following follows solely from using the expansion (9) in combination with Eq. ( 6).Its importance relies on the fact that it removes the need for the spectral properties of H(λ) in determining the CD term and shows that the operators in the AGP are generated from the nested commutators L 2k−1 ∂H at odd orders.We note that the variable s in the integral representation represents a fictitious time.The unitary e −iH(λ)s is interpreted as the time evolution operator in the fictitious time with no need for the time-ordered product, as H(λ) is independent of s.When we keep all possible operators generated from the nested commutators, the exact AGP can be obtained by solving the equation for α nc k from Eq. ( 6).Practically, a truncation of the operator series yields an approximate AGP.The infinite series by nested commutators produces the same type of operators many times, and it is not clear how many terms should be kept to obtain a required accuracy.To treat the AGP systematically, we rearrange the expansion in Eq. ( 9) and represent the AGP in a finite series by using a set of orthonormal Krylov basis elements.

IV. KRYLOV EXPANSION
A. Inner product, basis operators, and vector representation of operators In the Krylov method [66], we use a set of operators satisfying an orthonormal relation.To define the orthonormality of operators, we first introduce the inner product for an arbitrary pair of operators X and Y as The operators are not necessarily Hermitian.In addition, the measure ρ(H) is a positive-definite Hermitian operator but not necessarily normalized.We note that the present method is applicable even when the Hilbert space dimension is infinite and the energy spectrum is continuous, provided that ρ(H) is chosen appropriately.We see in the following that the result of the AGP is independent of the choice of ρ(H).What is important is that ρ(H) commutes with H.We have for Hermitian operators X and Y .
To find an explicit representation of the superoperator L, we introduce a set of basis operators X = (X 1 , X 2 , . . .), that are Hermitian and orthonormal with each other: The number of operators is not specified here and is discussed in the following after clarifying the aim of the analysis.Generally, for a given quantum system, it is equal to or smaller than the square of the dimension of the Hilbert space.One of the aims of introducing the basis operators is that the superoperator L can be represented by an antisymmetric Hermitian matrix L. It has elements and satisfies L † = L and L T = −L.The diagonal components are equal to zero and each of the off-diagonal components is purely imaginary.Corresponding to the matrix representation of superoperator, a vector represents an operator.We write the Hamiltonian and the AGP Then, we obtain a vector representation of Eq. ( 6) as It is not a difficult problem to obtain the formal solution of this equation by using the spectral representation of L. However, L is generally a matrix of large size and the diagonalization is much more difficult than that of the Hamiltonian H.We resolve this problem in the following by introducing an algorithmic method.

B. Lanczos algorithm and Krylov basis
Equation ( 17) implies that the AGP |A⟩ is constructed from a linear combination of L n |∂H⟩ with n = 1, 2, . . . .We prepare the normalized vector |θ 0 ⟩ from the relation The coefficient b 0 represents the normalization factor and is written as The zeroth-order normalized vector |θ 0 ⟩ and the coefficient b 0 are defined for each component of λ.The same applies to the quantities introduced in the following.We abbreviate the component index to simplify the notation.Then, the new normalized basis |θ 1 ⟩ is defined from By construction, |θ 1 ⟩ is orthogonal to |θ 0 ⟩.We repeat the same procedure by using the relation with n = 2, 3, . . . .The positive coefficient b n is chosen so that |θ n ⟩ is normalized.Thus, the introduced vectors satisfy the orthonormal relation ⟨θ m |θ n ⟩ = δ m,n .When the dimension of the Hilbert space is finite, the number of basis elements must be finite, which means that there exists an integer The number of the basis vectors is given by d, which we refer to as the Krylov dimension.This way of constructing a basis set is nothing but the Lanczos algorithm since the matrix L is brought to a tridiagonal form T , satisfying We can also write Generally, for a given matrix L and an initial basis element |θ 0 ⟩, we can render the matrix in tridiagonal form algorithmically.We find in the following that the present choice of the initial basis in Eq. ( 18) is convenient to solve Eq. ( 17).
The introduction of the orthonormal basis vectors corresponds to that of the orthonormal basis operators |O n ⟩ = |θ n ⟩.In the original representation, with n = 0, 1, 2, . . ., d − 1.They are generated by the procedure and satisfy . This set of operators represents the Krylov basis.In the present choice of O 0 , the operators of even order O 2k (k = 0, 1, 2, . . . ) are Hermitian, and those of odd order O 2k−1 (k = 1, 2, . . . ) are anti-Hermitian.We note that the introduction of the basis operators X is not necessary, since we can construct the Krylov basis directly from Eq. ( 25).The introduction of the basis operators makes it clear that the introduction of the Krylov basis is equivalent to the Lanczos algorithm.
The following examples illustrate that the two options can prove convenient.
The advantage of the basis operator representation of L by L is that we do not need to calculate the nested commutators L n ∂H once we can construct a single matrix L. We also see that the number of the basis operators X is not necessarily equal to the square of the dimension of the Hilbert space d H .For the present purpose, we need operators in L n ∂H and the dimension of L, denoted by d L , satisfies Thus, the Krylov dimension d is defined by the minimum number of the basis elements.When the matrix L is block diagonalized, we may treat only the block in which the operators in ∂H are included.A good choice of the basis reduces the computational cost.It is known that the general upper limit of the Krylov dimension is given by the relation d ≤ d 2 H − d H + 1 [69].Generally, the Krylov method is useful when we treat the Heisenberg representation of a normalized operator O 0 , O(s) = e iHs O 0 e −iHs [66][67][68][69][70].We can represent the operator by a finite series as where φ n (s) is known as the operator wave function.The time dependence of O(s) can be conveniently studied by using the operator wave function, a feature we next apply to the computation of the AGP from the integral representation in Eq. (7).

C. Adiabatic gauge potential
We are now in a position to solve Eq. ( 6), or the equivalent Eq. ( 17), by using the Krylov basis.We use and solve the equation for {α k } d A k=1 .It is important to notice that A includes the Krylov basis at odd order, O 2k−1 .This property is a direct consequence of the representation in Eq. ( 9).The number of the operators is denoted by d A and is related to the Krylov dimension d as We first consider the case of even d.In this case, one finds Setting each side of this equation to zero yields where That is, we can find the AGP satisfying the relation |∂H⟩ − iL|A⟩ = 0, which is a suffi-cient condition of Eq. ( 17).We also see that the relation |∂H⟩ − iL|A⟩ = 0 represents the equation for a dynamical invariant, when the eigenvalues of the Hamiltonian, ϵ n [λ(t)], are time independent [89,90].In this case, diagonal components of ∂ λ H(λ) in the eigenstate basis are equal to zero, i.e., ⟨n(λ)|∂ λ H(λ)|n(λ)⟩ = 0. Next, we consider the case of odd d.In this case, an additional term appears in Eq. ( 30) and no solution exists for |∂H⟩ − iL|A⟩ = 0. We examine L 2 |A⟩ = −iL|∂H⟩ to find the expression Inverting the matrix in this expression, we can obtain the explicit form of the AGP.In the following, we solve this equation by using a different approach which proves illuminating.
We conclude this part by stating that the AGP can be constructed systematically by using the Krylov basis.The AGP is represented by an expansion of the Krylov basis, and the coefficient of each term is obtained as a function of b 0 , the scale of ∂H, and the set of Lanczos coefficients {b n } d−1 n=1 .When d is even, the instantaneous eigenvalues of the Hamiltonian must be time independent.Conversely, the Krylov dimension d is even (odd) when the eigenvalues of the Hamiltonian are time independent (dependent).The flowchart of the algorithm is presented in Fig. 1.Equation ( 28) is compared with Eq. ( 9).The former is expanded by orthonormal operators and the total number of series elements is finite, if the resulting AGP is given by a finite number of operators.The ex-pansion is also applicable to systems with a continuous spectrum.Thus, the Krylov method offers a general systematic method for constructing the AGP.

D. Operator wave function and adiabatic gauge potential
The AGP is closely related to the operator wave function φ n (λ, s) defined from the Heisenberg representation (34) where the initial condition is chosen as b 0 (λ)O 0 (λ) = ∂ λ H(λ).Substituting this representation into Eq.( 7), we obtain for k = 0, 1, . . ., d A , and for k = 1, 2, . . ., d A .This relation between φ n (λ, s) and α k (λ) shows that the latter is obtained from the Laplace transform of the former.The behavior of the operator wave function has been studied in the context of the Krylov complexity, and we can exploit the properties obtained in that context [66][67][68][69][70][71].For example, the operator wave function |φ(λ, and the initial condition |φ(λ, 0)⟩ = (1, 0, 0, . . . ) T .Here, iB is related to the matrix T in Eq. ( 22) under a unitary transformation.Since the equation for |φ(λ, s)⟩ is interpreted as a Schrödinger equation with a Hamiltonian iB(λ) independent of the fictitious time s, the solution is obtained by solving the eigenvalue problem The form of the Hermitian matrix iB indicates that the eigenvalues come in pairs ±ω n where ω n ̸ = 0, and the zero-eigenvalue state exists only when the size of the matrix d is odd.We refer to the details on the pairing of eigenstates in Appendix A. Here, we look at only the zero-eigenvalue state |ϕ(λ)⟩ satisfying B(λ)|ϕ(λ)⟩ = 0 for odd d.We can solve the eigenvalue equation to obtain the normalized solution Since the matrices L and B are constructed from the commutator L(•) = [H, (•)], the eigenvalues are related to the energy eigenvalue difference ϵ m −ϵ n .The zero-eigenvalue state of M implies the existence of the diagonal components The contribution on the right-hand side is absent for even d with ∂ϵ n = 0. We can also use the equation for |φ(λ, s)⟩ to obtain the explicit form of α k .The 2kth component of the equation is given by Using the integral representation in Eq. ( 7), we obtain The left-hand side is calculated by using the integration by parts to give where the second term exists only for odd d and we write |ϕ⟩ = (ϕ 0 , ϕ 1 , . . . ) T .In the odd-d case, we obtain with k = 1, 2, . . ., d A − 1.It is not a difficult task to confirm that this relation is consistent with Eq. ( 33).The use of the operator wave function also allows us to obtain where |0⟩ = (1, 0, . . ., 0) T and Q = 1 − |ϕ⟩⟨ϕ| represents the projection operator onto the nonzero-eigenvalue states.We show the derivation in Appendix A. This representation is useful when we evaluate the norm of the AGP.
It is instructive to compare the present result for the odd-dimension case to that for the even-dimension case.Equations ( 31) and (32) show that, when the Krylov dimension is even, each order is calculated without using the higher-order contributions.This property is practically useful for systems with many degrees of freedom.
As we discuss in the next sections, we frequently consider the truncation of the series expansion as an approximation.By contrast, for an odd Krylov dimension, all the Lanczos coefficients are required to construct each term of the AGP, as we see in Eqs. ( 43) and ( 44) and the zero mode |ϕ⟩ in Eq. ( 39).However, each component of |ϕ⟩ takes a small value and could be negligible for large systems.
We can estimate a contribution from each term of the expansion in Eq. ( 28) by the Lanczos coefficient.When b n is an increasing function with respect to n, the corresponding α k is a decreasing function.The typical global behavior of the Lanczos coefficients has been discussed in many-body systems.It was found that b n ∝ n for chaotic systems and leads to a maximal pace of operator growth.Likewise, b n ∝ √ n for integrable systems, and b n ≈ const for noninteracting systems [67].In Fig. 2, we show the behavior of α k in the case of a linear and square-root growth of b n .The constant case is found in the examples in Sec.VI.In the figure, we also show a special case b n ∝ n(d − n) where the operators defined from the Krylov complexity theory form a SU(2) algebra [67][68][69][70][71][72].We also note that the series of Lanczos coefficients typically shows an oscillating behavior, as shown in the examples below.Given that the coefficients α k in the AGP expansion involve the ratio b 2k /b 2k+1 as in Eq. ( 44), a regular oscillation series of b k leads to a decreasing series on α k .These observations indicate that the property of the CD term is closely related to that of the operator growth in the Krylov subspace.

E. Classification of basis operators
It is instructive to notice that the AGP consists of the nested commutators at odd orders.When the original Hamiltonian is real symmetric, the nested commutators at even orders L 2k ∂H are real symmetric and those at odd orders L 2k−1 ∂H involve the imaginary unit.This means that the basis operators are classified into two parts: X represents basis at even orders and Y at odd orders.These Hermitian operators satisfy LX ∈ Y and LY ∈ X. Accordingly, the matrix L, the basis operator representation of L, takes the form where M µν = (X µ , LY ν ).We note that M * νµ = (Y µ , LX ν ).The size of the matrix M is determined by the numbers of the basis operators d X and d M is generally a rectangular matrix, and the Lanczos algorithm is applied for even d = 2d A as where each of {|θ 2k ⟩} d A −1 k=0 has d X components, each of {|θ 2k−1 ⟩} d A k=1 has d Y components, and the size of the lower triangular matrix on the right-hand side is d A × d A .Since the numbers of the basis operators must be large enough to span the operator space, we find d X ≥ d A and d Y ≥ d A .In the case of odd d = 2d A + 1, M is decomposed as where each of k=1 has d Y components, and the size of the matrix on the right-hand side is (d A + 1) × d A .We also find d X ≥ d A + 1 and d Y ≥ d A .In this case, the minimum number of d X is larger than that of d Y .

F. Relation to the variational method
The orthonormal relation of operators is useful to understand the relation between the present method and the variational method [40].In the variational method, the AGP is obtained by minimizing the cost function for a given operator ansatz of A with undetermined coefficients.In our notation, this can be written as with ρ(H) = 1.One of the systematic methods for obtaining the AGP is to use the nested commutator series in Eq. ( 9) and carrying out the minimization procedure with respect to the coefficients α nc k (λ) [41].Practically, the number of series elements in Eq. ( 9) is restricted to a finite value, and the approximate AGP is obtained from the minimization.
It is not obvious that the variational method can give the exact AGP even when all of the possible operators are incorporated in the trial form of the AGP by nested commutators.When the AGP satisfies |∂H⟩−iL|A⟩ = 0, which is a sufficient condition of Eq. ( 17), we have discussed that the Krylov dimension is even and that the matrix L as well as B are invertible.Then, the minimization procedure gives the exact AGP |A⟩ = −iL −1 |∂H⟩.
On the other hand, when the Krylov dimension is odd, L is not invertible and special care is required for the zero-eigenvalue state.The zero-eigenvalue state of L denoted by |ϕ L ⟩ is obtained in the same way as that of B, Eq. ( 39), and is written by a linear combination of even basis {|θ 2k ⟩} . It is orthogonal to the AGP as ⟨ϕ L |A⟩ = 0, and the cost function is decomposed as where P = |ϕ L ⟩⟨ϕ L | and Q = 1 − P are projection operators.The first term does not affect the variational procedure, and the minimization of the second term gives This is the solution of L (|∂H⟩ − iL|A⟩) = 0, which means that the variational method gives the exact AGP.We note that the trial form of the AGP must include all possible operators from the nested commutators at odd order to find the exact AGP from the variational method.When we consider a restricted number of operators, the minimization gives an approximate AGP.This procedure is essentially equivalent to considering a restricted number of basis operators for the Krylov expansion.However, as we explicitly show in the following examples, the variational procedure does not necessarily lead to the result from the Krylov expansion.This is because the coefficients of the AGP in the variational method are optimized in the truncated space.
A significant fact is that we can find the exact AGP associated with a generalized cost function of the form This form is useful when the dimension of the Hilbert space is infinite and when the spectrum is continuous.
Although the exact AGP must be independent of ρ(H), the approximate AGP is generally dependent on it.In the variational method, we usually set a constant ρ(H).
It may be possible to use a different ρ(H) for the variational calculation.However, even when all possible operators are incorporated, it is not evident that the variational method gives the exact result, which is independent on the choice of ρ(H).The Krylov method states the requirements for ρ(H) explicitly and clarifies that the result is independent on that choice.

V. APPLICATIONS TO SMALL SYSTEMS
In the construction of STA, one is interested in the time dependence of the Hamiltonian.We set λ(t) = t and identify λ as time t.Then, the AGP is equivalent to the CD term.In the applications discussed below, we write the Hamiltonian as H(t) and use the CD term H CD (t) instead of the AGP A(λ).For the small systems discussed in the present section, it is not a difficult task to calculate the CD term explicitly.We study how the CD term is obtained by the Krylov method in typical small systems.

A. Two-level system
The study of STA by CD driving in the canonical twolevel system [28][29][30][31] was soon followed by its experimental demonstration [19,24,91], often in a rotating frame, i.e., making use of a unitarily equivalent CD Hamiltonian.To illustrate the engineering of STA in Krylov space, consider the two-level Hamiltonian in terms of the positive scalar h, the unit vector n = (n 1 , n 2 , n 3 ), and the vector Σ = (X, Y, Z) with Pauli operators as entries.In this case, we have essentially no other choices than to set the basis operators as (X, Y, Z).We choose ρ(H) = 1/2 for the inner product.We assume that n(t) depends on t; otherwise, the CD term trivially gives zero.However, the explicit parameter dependence of h(t) is not necessary, as h determines only the overall scale of the Hamiltonian and the resulting CD term is independent of h.
It is a simple task to calculate the L matrix explicitly.We have Then, we set the initial basis vector to generate the Krylov basis and the Lanczos coefficients We find LO 2 − b 2 O 1 = 0, which shows d = 3 and d A = 1.
For ḣ ̸ = 0, the Krylov dimension d = 3 equals the number of basis operators.It is reduced to d = 2 and d A = 1 for ḣ = 0 where b 2 = 0. We note that the eigenvalues of the Hamiltonian, ±h/2, are time independent when ḣ = 0.In any case, the dimension of the AGP is given by d A = 1.We find the CD term This result is consistent with the known result [28][29][30][31].

B. Driven harmonic quantum oscillator
The driven harmonic oscillator is a workhorse in nonequilibrium quantum dynamics and, not surprisingly, has played a key role in the development of STAs [12,92] and their experimental demonstration [16,26].Although the dimension of Hilbert space is infinite, it is not difficult to treat the system analytically, since the system is a single-particle one, and the spectrum is discrete.In addition, its dynamics is described by a closed Lie algebra [93].

Consider the Hamiltonian
where Q and P are the position and momentum operators, respectively.Modulations of the time-dependent frequency ω induce expansions and compressions, while transport processes are associated with variations of the trap center q 0 [26,94,95].We use the creationannihilation operator representation where The CD term, in this case, is given by [54,55,92] Thus, the CD term involves a term proportional to the momentum operator, the generator of spatial translations, and a second term proportional to the squeezing operator, which is the generator of dilatations.
In the present case, we can explicitly calculate all the nested commutators.As mentioned, using the basis operators X is unnecessary.We find for k = 0, 1, . . ., and For ω = 0, the eigenvalues of the Hamiltonian are time independent, and the Krylov dimension is given by an even number.The explicit form of the Krylov basis is given in Appendix B. It is instructive to see how the exact AGP in Eq. ( 67) is obtained in the expansion.In the case at d A = 2, the CD term is expanded as CD , and the first term where . This result shows that each term of the expansion is dependent on the definition of the inner product in Eq. (10).

C. STIRAP
As a practical application of three-level systems, we next discuss the stimulated Raman adiabatic passage (STIRAP) [96,97].It is a method of population transfer between two states.We introduce an additional state and apply two external pump fields to the system.The states are given by |1⟩, |2⟩, and |3⟩, and we consider population transfer between |1⟩ and |3⟩.The simplest STIRAP Hamiltonian is given by A typical protocol is given in Fig. 3. Since we treat three-level systems, the number of independent operators is given by eight, except the identity operator.We also see that the Hamiltonian is real symmetric and the L matrix is written as Eq. ( 47).Possible basis operators for ρ(H) = 1/2 are given by and Using this basis, we find Eq.( 47) with We note that the number of operators in X can be reduced to four, as we can understand from the general discussions in the previous section.Since they are not much different, we use the 5 × 3 matrix here.We apply the Lanczos algorithm for the given M matrix with the protocol in Fig. 3(b) to calculate α k shown in Fig. 4(a).The Krylov dimension is given by d = 7.
The expansion is compared with the exact result [98] H where θ(t) = arctan(ω p (t)/ω s (t)) and ϕ(t) = [arctan( ω 2 p (t) + ω 2 s (t)/δ)]/2.In the Krylov method, the CD term is given by the form µ , where and are plotted in Figs.4(b)-4(d).For short and large times, the adiabatic condition is approximately satisfied, and the CD term is well approximated by the first term of the Krylov expansion.

VI. APPLICATIONS TO MANY-BODY SYSTEMS
The exact AGP or CD term is known for a limited number of many-body systems, and we expect that the   Krylov method gives advantageous results that cannot be obtained from other methods.For many-body systems, the required number of operators is large, and it is still a formidable task to find the exact CD term even in the present method.In this section, we treat one-dimensional spin systems where the exact CD term is known for some examples.
A. One-dimensional transverse Ising model

Exact result without a longitudinal magnetic field
We first treat the Ising spin chain in a transverse field: Many spins are aligned in a chain, and the number of spins is denoted by n s .We consider the periodic boundary condition, and the subscript is interpreted as mod n s .We are interested in the large-n s limit where the system at g = 1 shows a quantum phase transition [99,100].
It is also known that the system is equivalent to the free fermion system [101,102].Then, the Hamiltonian is represented as an ensemble of two-level systems, and the CD term for each two-level system can be found from the result in the previous section.Here, to study properties for many-body systems, we do not use the mapping and treat the spin operators.
Under the setting ρ(H) = 1/(2 ns n s ), we define or-thonormalized operators with k = 1, 2, . . ., n s − 1.These operators take bilinear forms in fermion operators and can be used as basis operators.Since the Hamiltonian commutes with (−1) P = ns n=1 Z n , the number of independent operators is reduced.We have the relations , and W ns−k = (−1) P W k .Then, the subscript takes k = 1, 2, . . ., n s /2 for even n s .
Acting L to these operators, we obtain Since the Hamiltonian is real symmetric, the nested commutators at odd order involve only W k .In Appendix C, we show The Krylov dimension is odd in this case.Using the result b 0 = √ n s v ġ(t)/2 and b 1 = √ 2v, where we assume v > 0 and ġ(t) > 0, we can calculate all of the Lanczos coefficients.We plot the Lanczos coefficients in the left in Fig. (5) for several values of g.The asymptotic forms satisfy This result implies that we can generally find a phase transition point from the Lanczos coefficients.
The coefficients of the CD term are obtained from Eq. (33).Since the matrix in the equation has the diagonal components 4v 2 (1 + g 2 ) and the nonzero off-diagonal components 4v 2 g, we can invert the matrix by using the discrete Fourier transformation.The details are given in Appendix C. We obtain The CD term is given by This result entirely agrees with that in Refs.[38,49,50].Since the Krylov basis at odd order in the present case is given by a simple form W k , it is natural to find the same result without using the Krylov method.
We plot α k for several values of g in the right in Fig. 5.For g ̸ = 1, the Lanczos coefficients at even order are larger than those at odd order; we see from Eq. ( 44) that α k is a decreasing function in k.Then, we find that α k rapidly decreases, which means that few-body interactions become the dominant contributions to the CD term.It was discussed in Refs.[38,50] that α k ∼ g k−1 for |g| < 1 and α k ∼ g −k−1 for |g| > 1 at the thermodynamic limit n s → ∞.For g = 1, the Lanczos coefficients take a constant value, and α k decreases slowly as a function of k.At n s → ∞, an infinite number of operators contributes to the CD term as discussed in Ref. [38].

Approximate result with a longitudinal magnetic field
The free fermion representation is possible only when the direction of the magnetic field is perpendicular to the direction of the interaction operator.Here, we treat the nonintegrable Hamiltonian with g(t) = t/t f .This is the standard form of the quantum annealing Hamiltonian and a test bed for universal dynamics of quantum phase transitions with a bias [103].
We consider the time evolution from the ground state of H(0) toward that of H(t f ).The nonadiabatic effect makes the system deviate from the instantaneous ground state, and we apply the CD term to prevent the transition.
It is a complex problem to find and implement the exact CD term for this Hamiltonian, and we consider a restricted number of operators.As we discussed in the above example, each term of H CD (t) involves an odd number of Pauli Y operators.We keep the Y basis up to two-body operators with ρ(H) = 1/(2 ns n s ).We act with L on these operators to construct the M matrix.Using the X basis we can construct the M matrix with the size 9×3 and apply the Lanczos algorithm to calculate the approximate CD term.We note that the size of the M matrix can be reduced to 3 × 3, keeping the result unchanged.The Krylov dimension is given by d = 7.We calculate the approximate CD term and numerically solve the time evolution with and without the approximate CD term by setting the initial state as the ground state of H(0) for the system size N = 6.In Fig. 6, we plot the approximate CD term as a function of t and the fidelity The time evolution for the quantum annealing Hamiltonian in Eq. ( 92).We set N = 6, γ/v = 1.0, and h/v = 1.0 for the left and h/v = 0.1 for the right.In the upper, we plot (a1(t), a2(t), a3(t)) for the approximate CD term In the lower, we plot the fidelity f for the time evolutions with H(t) (red open circle) and with H(t) + HCD(t) (blue filled circle).
as a function of the annealing time t f .Here, |ψ(t)⟩ represents the time-evolved state with or without the approximate CD term, and |ψ gs (t)⟩ represents the instantaneous ground state of H(t).f is close to unity when the adiabaticity condition is satisfied.
We numerically confirm that the present method gives the same result as the variational method, which is expected from the general discussion.When h/v takes a large value, the ground-state energy is well separated from the other ones, and we find a large fidelity.As h/v decreases to zero, the fidelity becomes smaller and finding the exact ground state becomes difficult even with the approximate CD term in use.As we see from the figure, many-body terms in the CD term become important for small h, and we require higher-order terms to improve the result by the approximate CD driving.

Krylov method
In the example of the Hamiltonian in Eq. ( 79), O 2k−1 incorporates k + 1-body interactions only, which is a specific property for the transverse Ising model.Here, to study more complicated situations, we treat the onedimensional isotropic XY model (XX model) with the open boundary condition This model is also known to be equivalent to a free fermion model.However, the coefficients are dependent on the site index n and it is generally a difficult task to find the exact CD term.The mapping to a bilinear form in fermion operators denotes that we can find a closed algebra within a limited number of operators.For ρ(H) = 1/2 ns , we define with n = 1, 2, . . ., n s and k = 1, 2, . . ., n s − n.Since the Hamiltonian is real symmetric, we define X = ({Z n }, {V k n }) and Y = ({W k n }) to construct the M matrix with the size n s (n s + 1)/2 × n s (n s − 1)/2.We have For the Hamiltonian coefficients {v n (t)} ns−1 n=1 and {h n (t)} ns n=1 , we consider an annealing protocol.We take v n as a constant value and consider the two cases: (i) uniform distribution v n = v 0 > 0 and (ii) random distribution v n = v 0 r n , where r n is a uniform random number with r n ∈ [−1, 1].For a given set of {v n }, we change the magnetic field h n (t) from h n (0) = h 0 ≫ |v n | to h n (t f ) = 0. We take h n (t) = h 0 (1 + tanh f n (t))/2 with f n (t) = n − 1 + x 0 − (n s − 1 + 2x 0 )t/t f .The system becomes adiabatic for large t f .The parameter x 0 takes a large value so that the conditions at t = 0 and t = t f are satisfied.We plot the protocols used in the following calculations in the upper in Fig. 7.The instantaneous eigenvalues for (i) and (ii) are, respectively, plotted in the lower in Fig. 7.The Hamiltonian commutes  with M = ns n=1 Z n , and we consider the block with M = n s − 2. The Hamiltonian takes a tridiagonal form, and the size of the matrix is given by n s .
Although it is not difficult to implement the Lanczos algorithm numerically for a considerably large value of n s , we here take n s = 6 to keep a good visibility of the plotted points.In Fig. 8, we display the Lanczos coefficients and the coefficients of the CD term for a fixed t.The Lanczos coefficients show an oscillating behavior.For the uniform case (i), the even order tends to be larger than the odd order, and correspondingly |α k | shows a decreasing behavior.For the random case (ii), we observe a complicated behavior denoting that the higher-order contributions of the Lanczos expansion are important.
The typical behavior of b n for a large n s is shown in Fig. 9.For any choice of parameters, we observe a flat band, which is considered to be a property for "noninteracting" systems.
To study the properties of the obtained CD term, we next decompose it.The Lanczos basis at odd order is written as where |θ which weights the contribution of the p-body term.This is a function of t and is plotted for each value of p in Fig. 10.We can understand from the comparison between the energy levels in Fig. 7 and q (p) in Fig. 10 that the many-body interaction terms cannot be neglected when the energy levels significantly change as a function of t.
It is generally understood from Eq. ( 5) that the CD term gives a large contribution when some of the energy levels are close to each other.We calculate the amplitude of the CD term ⟨H CD |H CD ⟩, which is plotted in Fig. 11.As we have discussed, the CD term is decomposed as CD and In the same figure, we plot the result where the first term of the expansion is kept for each decomposition.We also plot the first contribution of the approximate CD term H nc CD obtained from the expansion in Eq. ( 9).The coefficient α nc 1 is obtained from the minimization condition of Tr[( Ḣ − iLH nc CD ) 2 ] [ 40,41].The result implies that the approximate CD term underestimates an abrupt growth of the CD term.The norm of the CD term for the XX model at ns = 6.Bold solid curves "exact" (black) represent σ = ⟨HCD|HCD⟩.We use the replacement HCD → H (2) CD for dashed curves "2-body" (blue), HCD → ib0α1O1 for dotted curves "k = 1" (green), and HCD → iα nc 1 L Ḣ for thin solid curves "variational" (red) where α nc 1 is obtained from the variational method.
We note that this result does not necessarily lead to the failure of the approximation method.The CD term is independent of the choice of the initial condition of the time evolution.In the present examples, as we see from Fig. 7, the energy level crossings occur at higher energy levels.The ground-state level is isolated from the other levels, and we can expect that a large amplitude of the CD term is not required for the control of the ground state.The situation is, in this sense, opposite to that across a quantum phase transition, discussed in Sec.VI A 1, in which the gap between the ground state and the first excited state closes, and the norm of the CD exhibits a singularity [38,50].
In the present analysis, we set the measure ρ(H) in the inner product as ρ(H) = 1/2 ns .When we control a state with an energy level well separated from the other levels, it is reasonable to consider a weighted measure such as the Gibbs-Boltzmann distribution ρ(H) ∝ e −βH .Although the exact CD term is independent of the measure, each term of the expansion Therefore, when we consider a truncation approximation in the Krylov expansion, the choice of the measure strongly influences the result.Since a nontrivial choice of the measure makes the calculation of the inner product difficult, it is practically an interesting problem to find a convenient form of the measure.For the ground state, it is tempting to take the limit β → ∞ for ρ(H) ∝ e −βH .However, the measure is not positive definite in this limit, and we find unexpected behavior, such as the vanishing of the Lanczos coefficient b n for n < d − 1.

Toda equation
It is known that the exact CD term is analytically obtained when the coefficients of the Hamiltonian satisfy the Toda equations [104] ḣn The CD term is then given by with W 1 n as in Eq. ( 98).This term satisfies Ḣ −iLH CD = 0, which means that the instantaneous eigenvalues of the Hamiltonian are time independent.
Despite this simplicity, the corresponding Krylov expansion is generally involved and is required at each time.The time evolution of the Hamiltonian is reflected only in the choice of the initial basis b 0 O 0 = Ḣ, and the expansion is essentially insensitive to the choice.Except for the special cases discussed below, we find that the Krylov dimension is as large as the number of odd basis elements n s (n s − 1)/2.Highly nontrivial cancellations should be observed when we calculate α k to give the result with q (p) = δ p,2 .
As a very special case, we can find the result with d = 2 and d A = 1 when the coefficients are written as Using the Toda equations, we obtain that a single equation describes the time evolution For a given h 1 and a initial condition θ(0), the coefficients evolve, keeping the equidistant of h n (t) and a quadratic form of v n (t).In this case, we find that the first-order term in the Krylov expansion is proportional to the exact CD term, i.e., Since LO 1 − b 1 O 0 = 0, the expansion terminates at this order, and we obtain a simple result with d = 2.It was discussed that the present choice of parameters saturates the operator speed limit, i.e., the quantum speed limit in unitary operator flows [72].The nested commutators span the Krylov space within a limited number of operators.

VII. DISCUSSION AND SUMMARY
The use of the integral representation of the CD term introduced in Ref. [41] has eased the study of STA in many-body systems by removing the requirement for the exact diagonalization of the instantaneous system Hamiltonian.In its place, the CD term can be expressed as a series of nested commutators that follows from the Baker-Campbell-Hausdorff formula.The coefficients in such an expansion can be determined through a variational principle [40].
In this work, we have introduced Krylov subspace methods to provide an exact expression of the CD term.The Krylov algorithm identifies an operator basis for the terms generated in the series by nested commutators, along with the set of Lanczos coefficients.Using these two ingredients, we have provided an exact closed-form expression of the CD term, circumventing the need for a variational approach.When the dimension of the Hilbert space is finite, the series by the Krylov basis is finite, which is in contrast to the expansion using nested commutators.
We have shown the applicability of our method in the paradigmatic models in which the CD term admits an exact closed-form solution.This includes single-particle systems such as two-and three-level systems and the driven quantum oscillator.Although the applications of the present method can be laborious, the implementation of our approach in these systems is a straightforward task, applicable to any Hamiltonian with no special symmetry.We have further applied our formalism to a variety of quantum spin chain models encompassing the cases in which the system is integrable, nonintegrable, and disordered.Specifically, we applied the construction of the CD term in Krylov space to the one-dimensional transverse-field quantum Ising model as a paradigmatic instance of an integrable and solvable quasi-free fermion Hamiltonian.We have further demonstrated our approach in the presence of a longitudinal symmetry-breaking bias field that breaks integrability.A similar study is possible for the XX model where the free fermion representation is available.However, in that case, the site-dependent couplings make the explicit construction of the CD term by the standard method difficult, except for the special case when the coupling constants are varied in time according to a Toda flow, and the resulting CD term takes a simple local form.We have explicitly constructed the CD term for general and disordered couplings, and the result was compared to the approximation method [40,41].
The main task of the Krylov algorithm is to construct the basis operators for the minimal subspace in which the dynamics unfolds and to determine the associated Lanczos coefficients by iterations.The CD term is constructed as a series involving only the odd-order operators of the Krylov basis, with the corresponding coefficients in this compact expansion being determined in terms of the Lanczos coefficients.These properties imply that we can find some implications by comparing the Lanczos coefficients of odd order and those of even order.As we see from Figs. 5 and 8, the Lanczos coefficients typically show an oscillating behavior.When the coefficients at even order are larger than those at odd order, the corresponding coefficients of the CD term show a decaying behavior, and the CD term is well approximated by the first several terms of the expansion.We also find that the Krylov dimension is even when the instantaneous eigenstates of the original Hamiltonian are time independent.Thus, we can directly find the dynamical properties of the system from the Krylov algorithm.
The Krylov expansion is dependent on the choice of the inner product.Although the CD term is independent of the choice, each term in Eq. ( 28) is sensitive to it, a feature that can be relevant when considering truncating approximations.It is an interesting problem to find a proper choice depending on the situation to treat.We can also consider the truncation of the Krylov subspace.To this end, it suffices to restrict the basis operators and to construct the L matrix in the truncated space.Then, we can apply the Lanczos algorithm to find an approximate CD term.We note that even in that case the coefficients of the CD term are obtained without using any variational procedures, which gives a different result from the variational method.
Our primary emphasis has been on exact and analytical results formulating the CD term in Krylov space.In addition, there exist powerful numerical algorithms that largely simplify the computation of the Lanczos coefficients, used in our methodology.These are well established in the literature on Krylov subspace methods and numerical analysis [65].They are further available in popular software and numerical routines.They hold for any matrix of Hessenberg form and replace the Gram-Schmidt diagonalization by the use of Householder reflections, making the implementation numerically stable and computationally efficient.
Our work offers an interesting prospect to improve state-of-the-art quantum algorithms by combining the formulation of the CD term in Krylov space with the digital approach for quantum simulation.This approach may prove advantageous over the current formulation relying on variational methods [43][44][45][46][47][48], suggesting the need to generalize the error scaling in digitizing the CD term [105] to Krylov space.
In summary, we have proposed a technique for constructing the CD term exploiting the Krylov operator space.The method is flexibly applied to systems with many degrees of freedom and can be a powerful general method for understanding the dynamical properties of the system in control.Suppression of nonadiabatic transitions in discrete systems with many degrees of freedom is one of the dominant problems in quantum technologies, such as quantum simulation and quantum computing, and we hope that our method will be an efficient technique inspiring further studies.
Note added.-Recently, related results were reported in Ref. [106].and The Krylov basis is constructed by choosing the initial normalized vector |θ 0 ⟩ = (x, y, z) T .We obtain (B19) The expansion terminates at the fifth order, which means d = 5 and d A = 2 assuming x, y, and z are nonzero.As we discuss in the main body of the paper, the condition q0 = 0 gives y = 0, and we obtain d = 3 and d A = 1.For ω = 0, one finds x = z = 0 and obtains that d = 2 and d A = 1.
Appendix C: Krylov basis for the one-dimensional transverse-field Ising model For the Hamiltonian in Eq. ( 79), we apply the Krylov algorithm to find ) For a given set of Lanczos coefficients, α k is obtained by solving Eq. (33).Here, we denote the matrix in the equation by K and represent its spectral decomposition as K =

Figure 1 .
Figure 1.The flowchart of the Krylov algorithm to obtain the AGP A(λ) for a given Hamiltonian H(λ).

Figure 5 .
Figure 5.The Lanczos coefficients (left) and the coefficients of the CD term (right) for the one-dimensional transverse Ising model in Eq. (79) at ns = 200.For the Lanczos coefficients, we denote b1, b3, . . .by the filled circle and b2, b4, . . .by the open circle.

Figure 8 .
Figure 8.The Lanczos coefficients bn and the coefficients of the CD term, α k , for the XX model with ns = 6 at t/t f = 0.5.The top are results for case (i), uniform distribution, and the bottom are for (ii), random distribution.For the Lanczos coefficients, we denote b1, b3, . . .by the filled circle and b2, b4, . . .by the open circle.

Figure 9 .
Figure 9.The Lanczos coefficients of the XX model with ns = 100.We consider the random case (ii) and take t/t f = 0.5.

Figure 10 .
Figure 10.Time dependence of the norm fraction {q (p) (t)} ns p=2 in Eq. (101), weighting the contribution of the p-body interactions to the CD term, for the XX model with ns = 6.Many-body CD terms are shown to be necessary when the energy levels change significantly along the driving protocol.

1 ⟩
has n s + 1 − p components and the corresponding Krylov basis involves p-body interactions.We decompose the norm of the CD term as ⟨H CD |H CD ⟩ =

Figure 11 .
Figure 11.The norm of the CD term for the XX model at ns = 6.Bold solid curves "exact" (black) represent σ = ⟨HCD|HCD⟩.We use the replacement HCD → H )where ⟨•⟩ denotes the average Tr[ρ(H)(•)].Since the dimension of the Hilbert space is infinite in the present system, we need to choose the density operator ρ(H) in the inner product in a proper way.For example, we can use the canonical Gibss-Boltzmann distribution ρ(H) = e −βH /Tr e −βH .Now, we obtain