Inertial Theorem: Overcoming the quantum adiabatic limit

We present a theorem describing stable solutions for a driven quantum system. The theorem, coined inertial theorem, is applicable for fast driving, provided the acceleration rate is small. The theorem states that in the inertial limit eigenoperators of the propagator remain invariant throughout the dynamics, accumulating dynamical and geometric phases. The proof of the theorem utilizes the structure of Liouville space and a closed Lie algebra of operators. We demonstrate applications of the theorem by studying three explicit solutions of a harmonic oscillator, two-level and three-level system models. These examples demonstrate that the inertial solution is superior to that obtained with the adiabatic approximation. Inertial protocols can be combined to generate a family of solutions. The inertial theorem is then employed to extend the validity of the Markovian master equation to strongly driven open quantum systems. In addition, we explore the consequence of geometric phases associated with the driving parameters.


I. INTRODUCTION
Closed-form solutions for the propagator are of utmost importance for quantum control [1][2][3][4][5][6]. Generally, any Hamiltonian that allows control does not commute with itself at different times, leading to a time-ordering procedure in the evolution operator [7,8]. Hence, the development of closed solutions, for control Hamiltonians, faces the formidable task of time ordering [9,10]. Moreover, one desires that such solutions are stable under variations in the driving protocol and external noise.
The present paper is devoted to the construction of closedform, stable solutions of driven quantum systems. Currently, the popular theoretical as well as experimental approach utilizes the adiabatic theorem, constrained to slow driving [11][12][13][14][15]. Here, we propose approximate solutions for rapid processes based on the inertial theorem. In the appropriate limit, these solutions incorporate the adiabatic approximation.
The inertial theorem utilizes a timescale separation between variables to generate the system's evolution for slow acceleration of the external driving. The derivation subsides in the Liouville space and requires the existence of a timedependent operator basis. Formally, the theorem is similar to the adiabatic theorem, where adiabatic states are replaced by the time-dependent eigenoperators of the generator of the propagator [16]. The inertial solutions remain precise for rapid driving of the system under the condition of slow acceleration relative to the system's dynamics.
Inertial protocols can generate a diverse family of solutions, as the fast degrees of freedom are arbitrary smooth functions while only the rate of change of the slow variables is restricted. These solutions identify invariant operators, which are time-dependent constants of motion [17,18]. Moreover, inertial solutions can be combined to generate an inertial Hamiltonian, extending the family of possible solutions.
We demonstrate applications of the inertial theorem by studying three physical models: a time-dependent harmonic oscillator, a driven two-level system, and a three-level system. These models are the building blocks of both experimental and theoretical studies performed in the quantum regime [19,29].
We utilize the inertial solution to derive the equations of a motion for a driven open quantum system. We then analyze the geometric phase associated with the inertial solution. This phase differs in its physical role from the Berry phase of the adiabatic solution. In contrast to the Berry phase, the former phase directly influences observables and can be witnessed for nonclosed circuits in the driving parameter space. For systems interacting weakly with the environment, the geometric phase induces a shift in the decay rates.

II. THE INERTIAL THEOREM
Obtaining a closed-form solution for dynamics generated by a time-dependent Hamiltonian is a difficult task.
The inertial theorem constructs a solution by incorporating the time dependence within a time-dependent operator basis and scaled time.
The derivation of the inertial theorem is conducted in Liouville space, a state space of system operators {X }, endowed with an inner product (X i ,X j ) ≡ tr(X i †X j ) [35][36][37]. In Liouville space, the system's dynamics are represented in terms of a basis of orthogonal operators {V }, spanning the space. For example, for a two-level-system, the operator basis can be the Pauli operators. An arbitrary operatorÂ, spanned by operator the basis {V },Â(t ) = N 2 i=1 a i (t )V i , is represented by the vector A(t ) = {a 1 (t ), a 2 (t ), ..., a N 2 (t )} T in Liouville space, where a i (t ) are time-dependent coefficients and N is the size of the Hilbert space.
Employing the Heisenberg equation of motion, the dynamics of systems in Liouville space are determined by where v the expansion elements in the basis {V }. The conventionh = 1 is used throughout this paper and the superscript H designates that the operatorsV i are in the Heisenberg picture, i.e.,V H i (t ) =Û † (t, 0)V iÛ (t, 0), whereÛ (t, 0) satisfies the Schrödinger equation with respect to the HamiltonianĤ (t ) [at time t = 0, the Schrödinger and Heisenberg pictures coincide, v We consider a finite time-dependent basis, forming a closed Lie algebra; this guarantees that Eq. (1) can be solved within the basis [38]. This property applies trivially for any finite Hilbert space or else when a closed subalgebra can be found, for example, the Heisenberg-Weyl group that defines the Gaussian states of the quantum harmonic oscillator [39].
It is useful to limit the description to the minimal subalgebra required to solve the system's dynamics. In the case of compact algebras, this greatly simplifies the analysis, while for noncompact algebras, finding a subalgebra is a prerequisite for constructing the inertial solution.
For a closed Lie algebra, Eq. (1) has the simple form where M is a N 2 by N 2 matrix with time-dependent elements and v H is a vector [40]. For compact algebra, M has real eigenvalues, for noncompact algebras complex eigenvalues are possible (see Sec. IV A) [41]. A formal solution for Eq. (2) requires a time-ordering procedure v H (t ) = T exp (−i t 0 M(τ )dτ ) v(0), where T is the chronological time-ordering operator. This formal expression is impractical, as it includes an infinite sum of integrals [9].
The current derivation bypasses the time-ordering procedure by the following strategy: We search for a driving protocol that allows solving Eq. (2) explicitly and then extend the solution to a broad range of protocols employing the inertial approximation. By choosing a suitable time-dependent operator basis, the generator of the dynamics in Liouville space can be expressed as Here, P ( χ ) is an invertible matrix (unitary for an Hermitian M), dependent on constant parameters {χ }. D is diagonal real matrix with time-dependent eigenvalues, which are a function of both {χ } and time-dependent parameters { (t )}. The parameters are expressed in short notation as χ = {χ 1 , χ 2 , ..., χ m } T and (t ) = { 1 (t ), 2 (t ), ..., N 2 (t )} T . In the first two examples presented, there is a single parameter χ = χ , which is equal to the adiabatic parameter μ.
More general examples include multiple inertial variables, cf.
Sec. IV C. We will prove in Sec. V that decomposition Eq. (3) can always be achieved for any time-dependent analytical Hamil-tonianĤ (t ). Nevertheless, the existence of a solution is not constructive, therefore the suitable time-dependent basis associated with a general protocolĤ (t ) is not straightforward. Once a suitable protocol and time-dependent operator basis is found, for which M obtains the required form, Eq. (3), the solution becomes Next, we write the eigenvalues of D as a product of a χ -dependent function and a time-dependent function, leading to D = diag(λ 1 ( χ ) 1 (t ), ..., λ N ( χ ) N (t )). Since λ j , j , and χ are not specified, such decomposition is general and does not enforce further restrictions on our result. The solution of where the scaled-time parameters are θ k (t ) = t 0 dt k (t ) and c k = i P ik are constant coefficients. The Liouville vector F k corresponds to the eigenoperatorF k = i P −1 kiV i , where P −1 ik are elements of P −1 . For a Hermitian M, the eignvalues λ k are either zero or are pairs with equal magnitude and opposite signs. The structure of Eq. (4) allows an explicit solution for the dynamics, including cases where the operator basis is time dependent. As a result, the solution circumvents the timeordering operation. However, the approach is limited by the condition that P is a constant unitary operator, i.e., χ = const. For a set basis of operators {V }, this condition restricts the relevant driving protocols.
For general protocols, when χ (t ) varies with time, the solution can be extended to processes of slowly varying χ (t ) (inertial). The entire proof is reported in Sec. III, and follows a mathematical construction similar to that of the adiabatic theorem [11,14,42].
For a state v, driven by a inertial protocol, the system's evolution is given by where the first exponent is determined by the dynamical phase and the second includes a new geometric phase: Here, G k are the biorthogonal partners of F k . The system's state follows the instantaneous solution determined by the instantaneous χ (t ) and phases associated with the eigenvalues λ k k and eigenoperatorsF k . We restrict the analysis to the case where λ k k do not cross, hence, the spectrum of D remains nondegenerate throughout the evolution. Substituting the inertial solution, Eq. (6), into Eq. (4) enables assessing the validity of the approximation. The quality of approximation is expressed in terms of the magnitude of the inertial parameter (Sec. III) which implies that the inertial solution, Eq. (6), remains valid when χ follows a path in the parameter space of {χ }, where the eigenvalues λ k and λ n are distinct [15]. The factorization of the eigenvalues of D recognizes two timescales. A short timescale related to the (rapidly changing) frequencies { (t )} and a long timescale associated with the change in {λ( χ )} variables. This separation of timescales allows containing the rapid change in the protocol within the inertial solution. The rapid time dependence of the solution is effectively absorbed into the frequencies and time-dependent operator basis. In addition, we identify the dynamical invariant operators [43][44][45], which are associated to the eigenvectors F k with vanishing eigenvalues, λ k = 0.
The inertial theorem incorporates the adiabatic theorem, since for slow driving { } are constant and {χ } are slowly varying, therefore the matrix M(t ) can be diagonalized at each instant. We then obtain the eigenoperators associated with M(t ) and the decomposition in Eq. (3).

III. INERTIAL THEOREM PROOF
The following derivation is in the spirit of the adiabatic theorem, as presented by Schiff [11], and the generalization for a non-Hermitian Hamiltonian has been presented by Ibanez [42]. We formulate the derivation in Liouville space, a Hilbert space of operators. These operators are defined in terms of the underlying Hilbert space of wave functionŝ X |ψ = |φ with a scalar product ψ|φ .
For a system described by a finite algebra of operators, the Liouville generator M( χ, ) is a diagonalizable rank N 2 , parameter-dependent matrix, where the elements of χ and are real parameters which can be viewed as coordinates of a parameter space. We assume the N 2 instantaneous right eigenvectors of the Liouville generator M are nondegenerate (at all times, i.e., there is no level crossing). These are denoted by F k ( χ ), k = 1, 2, ..., N 2 , and are associated with the eigenoperators of M which satisfy an eigenvalue equation [16], where β k are time-dependent complex functions [46]. For example, when decomposition Eq. (3) holds, the eigenvalues are β k (t ) = exp (−iλ k t 0 dt k (t ) ). In addition, for compact algebras, the matrix M is Hermitian and the left and right eigenvectors are conjugates.
We introduce a second set of biorthogonal partners { G( χ )}, satisfying The two sets are biorthogonal, meaning that ( G k , F n ) = δ kn , where (,) is the scalar product of the two vectors in Liouville space.
The set of instantaneous eigenvectors constitutes a complete basis of the Liouville space, allowing us to expand the quantum state in terms of the basis elements. We therefore propose a solution for Eq. (2), which is a superposition of the eigenvectors F k with additional dynamical and geometric phases, Eq. (6).
The orthonormal property of the eigenvectors, ( G k , F n ) = δ kn implies that for all n and k, therefore ( G k , ∇ χ F k ) is pure imaginary and, as a result, φ k are real. Similarly, by differentiating the identity ( G n , M F k ) = 0, for n = k, we obtain We proceed by inserting Eq. (6) into Eq. (2). We then project G k from the left and utilize the orthogonality property and the derived identities Eqs. (11) and (12) to obtain a set of differential equations, with ξ nk ≡ t 0 dt [λ n n − λ k k ] − (φ n − φ k ). Typically, the dynamical phase is the dominant contribution to ξ nk , allowing us to neglect the geometric phases. Integrating Eq. (13) and solving iteratively leads to The term ( d χ dt ) −1 diverges in the inertial limit, inducing rapid oscillations in the last term. Assuming the integrand of the last exponent is integrable in the interval [ χ (0), χ (t )], the Riemann-Lebesgue lemma implies that when the phase of the last exponent, or ξ nk , change rapidly relative to the integrand, the sum in Eq. (14) vanishes [47]. This implies that the inertial solution is valid when for all n = k.

IV. INERTIAL EXAMPLES
The first two examples, a driven harmonic oscillator and two-level system, demonstrate the theory in the simple framework of Liouville space, with a single inertial variable. These examples allow a closed-form analysis. The third example, a driven three-level system, is more involved, leading to a generalization of the STIRAP process.

A. Parametric driven harmonic oscillator
To demonstrate the inertial theorem, we begin by analyzing the dynamics of a driven harmonic oscillator. Physically, the system can be realized by a particle in a varying harmonic potential [29] and is represented by the Hamiltonian whereq andp are the position and momentum operators, m is the particle mass, and ω(t ) is the oscillator frequency. We consider an initial Gaussian state, which is fully defined by the set of time-dependent operators: with where χ = μ =ω ω 2 . When μ → 0, the adiabatic solution is exact.  The eigenoperators and eigenvalues,F k and λ k , are obtained by diagonalization (see Appendix D). The matrix B has real eigenvalues that possess a non-Hermitian degeneracy for χ = μ = 2 [41]. This limits the solution to avoid proximity to the degeneracy point. The inertial parameter, Eq. (8), takes the form ϒ ∼ (μ 2κω ) 2 , where κ = 4 − χ 2 , which explicitly becomes When ϒ 1, the inertial solution, Eq. (6), is a good approximation of the true dynamics. The quality measure A ≡ − log 10 (1 − F ), of the inertial solution as a function of time. As the fidelity converges to unity, A increases. The increase in the fidelity at small times can be explained by Eq. (20). As t f decreases χ (t ) becomes constant, (χ (t f ) → χ (0)). Calculation parameters for the HO are ω(0) = 20, ω(t f ) = 10, and a = −5 × 10 −3 . The parameters for the TLS are¯ (0) = 20,¯ (t f ) = 10, ε = 8, andā = −5 × 10 −3 , with initial state v(0) = {4, 1, 1} T .
For the demonstration, we consider a particle of mass m = 1 in a varying harmonic potential. The particle is initialized in the ground state ρ(0) = |0 0|, associated with the initial frequency ω(0) = 20.
The inertial approximation is evaluated by comparison to a converged numerical solution, denoted byρ N . The fidelities F of the inertial and adiabatic solutions are calculated in terms of the Bures distance with respect toρ [49]. The fidelities are compared in Fig. 1. For the analysis, we use the protocol ω(t ) = , which satisfies The protocol is designed by imposing ω(0), ω(t f ) and parameter a, while t f and χ (0) are adjusted accordingly. Modifying the protocol duration t f interpolates between the sudden and adiabatic limits [50].  The comparison indicates that the inertial approximation outperforms the adiabatic approximation. This is in accordance with Fig. 3, demonstrating that the inertial parameter is always smaller than the adiabatic parameter ϒ < μ. The evolution of the state is presented in Fig. 2 in terms of the expectation values Ĥ , L , and Ĉ . For small a, the inertial solution remains accurate relative to the converged numerical result, see Fig. 2(a). Increasing a leads to the breakdown of the inertial solution, as shown in Fig. 2(b). The stability of the inertial approximation can be checked by adding random noise at each time step to χ (0) and a in Eq. (20). As expected, the solutions were stable to Gaussian noise with a standard deviation in the scale of ∼0.1 · a. A general protocol ω(t ) with a duration τ can be expressed in terms of a dimensionless parameter s = t/τ . The adiabatic parameter then becomes μ(s) = (ω(s) − ω(0))/(ω(s)ω(0)sτ ). For a sufficiently large protocol duration the inertial condition is satisfied and the inertial solution constitutes an accurate solution for the system dynamics.
The relative accuracy of the associated inertial and adiabatic solutions can be compared using ϒ, Eq. (19) and μ. When the oscillator frequency is accelerated slowly (ω < ω 3 μ 2 , cf. Eq. (19)) the inertial parameter behaves as ϒ ∝ μ 2 . Hence, for slow acceleration, even in the adiabatic regime the inertial solution possesses superior accuracy.

B. Two-level-system model
The driven two-level system, characterized by SU(2) algebra, is utilized as an additional demonstration of the inertial theorem. The system Hamiltonian reads whereŜ i is the i = x, y, z spin operator and the time-dependent Rabi frequency reads¯ (t ) = ω 2 (t ) + ε 2 (t ).
The dynamics of the system is analyzed employing a timedependent operator basis {H ,L,C,Î}, Here, χ =χ =μ =ω ε−ωε 3 , whereμ is the adiabatic parameter of the two-level system. To transform Eq. (22) to the factorized form, Eq. (3), we introduce a scaled vector for which the dynamics obtains the desired form, d dt u H = −i¯ B u H . This procedure is not limited to the two-level system and relies on the fact that the identity I commutes with any operator. The inertial solution is obtained by diagonaliz-ingB, Eq. (23), leading to the form of Eq. (6) (Appendix D).
We consider a protocol with a constant ε and a linear change inχ ,χ (t ) =μ(t ) =χ (0) +ā · t. This leads to the following protocol: Using this protocol, the exact adiabatic and inertial solutions were calculated. The results are shown in Fig. 1(b), illustrating the superiority of the inertial solution over the adiabatic result.
An experimental verification of the inertial solution of a two-level system has been demonstrated experimentally, employing a ytterbium ion 171 Yb + in a Paul trap [51]. The inertial protocol was realized in the experiment, demonstrating high accuracy of the inertial solution with respect to the measurements. In addition, deviations from the inertial solution were explored. Showing that when the studied protocol slightly deviates from the inertial condition, the phase still follows the inertial solution. This phenomena is witnessed as well in Fig. 2 for the harmonic oscillator. The experiment confirms the stability of the inertial solution under external noise.

C. Three-level system model
The three-level atom is one of the most extensively studied driven systems [52][53][54]. This is an elementary example of SU(3) algebra, which is abundant in many branches of physics [55][56][57][58][59]. In addition, the model serves as a template for a basic experimental technique in atomic and molecular physics: STIRAP. Based on an adiabatic approach and a dynamical symmetry, the STIRAP is a technique to efficiently transfer population between two quantum states via an intermediate state [32][33][34]. Amendments to the adiabatic scheme have been suggested, forcing the system to follow the adiabatic evolution by adding counteradiabatic terms to the Hamiltonian [60,61]. Insight on this well-established system can be gained from the inertial approach based on new invariants.
The basic model considers an atom with allowed transitions between the first and second levels, as well as the second to third levels. The atom interacts with an incident electromagnetic field, coupling the levels (|1 ↔ |2 , |2 ↔ |3 ): E (z, t ) = E 12 e i(ν 12 t−k 12 z) + E 23 e i(ν 23 t−k 23 z) + c.c. Under the two-photon resonance condition (Gell-Mann symmetry for N = 3) [52,53,55,56], the Hamiltonian within the rotating frame approximation has the form with = 21 = 32 , where i j = ν i j − ω i j is the detuning between the laser frequency ν i j and the Bohr frequency ω i j . In the adiabatic STIRAP, remains constant throughout the procedure. The Rabi frequencies are defined as α(t ) = d 12 · E 12 (t ) and β(t ) = d 23 · E 23 (t ), where d i j is dipole moment between levels i and j. The two-photon resonance model describes either a or ladder linkage pattern, where the photon energies correspond to the energy gap between state |3 and |1 , that is, ω 31 = ν 21 − ν 23 for the system and ω 31 = ν 21 + ν 23 for the ladder linkage, Figs. 4(a)-4(c).
With the aim of constructing an inertial solution for Hamiltonian Eq. (25), we begin by analyzing this model in Liouville space. The SU(3) algebra is characterized by eight orthogonal operators along with the identity. As a consequence, the solution requires a methodical approach to obtain a suitable protocol and the eight time-dependent operators of the Liouville basis which satisfy the decomposition Eq. (3). For this end, we introduce a static operator basis for the SU(3) algebra, composed of traceless Gell-Mann operators {λ 1 , ...,λ 8 } and the identity [62]. As required, the time-dependent Hamiltonian is within the algebra and can be expressed in terms of the operator basis: In the next step, we search for a suitable unitary time-dependent transformation that defines the timedependent operator basis: where = α 2 + β 2 , leads to the desired decomposition of the equation of motion, Eq. (3) [52]. This transformation conserves the simple commutation relations between the Gell-Mann matrices while rotating the basis operators with the Hamiltonian. A similar choice was chosen by Eberly and Hioe [52], with the restriction α(t ) ∝ β(t ). The When χ 1 and χ 2 are constant, an exact solution is obtained by diagonalization. For slowly varying χ 1 and χ 2 , the evolution is approximated by the inertial solution Eq. (6), where λ k and F k are eigenvalues and eigenvectors of M, and k = for all k.

D. Adiabatic and inertial STIRAP in Liouville space
In the three-level basis, the STIRAP procedure is a technique to completely transfer population between the states |1 and |3 . This procedure is immune to losses from spontaneous emission originating from the intermediate state, |2 , and is robust under small variations of the experimental parameters [34]. It relies on the two-photon resonance condition, = 21 = 32 , and adiabatic dynamics. This is achieved when the adiabatic parameter is sufficiently small, χ 2 1. A dynamical symmetry viewpoint serves as an elegant approach to understand the STIRAP procedure. Within this framework, the STIRAP is a consequence of the dynamical invariance ofT 8 (t ) when the dynamics are sufficiently slow. In the adiabatic limit, (χ 2 → 0), implying that d dtT H 8 → 0, which means that its expectation value is constant [63]: Here, ρ i j = i|ρ| j , ρ ii is the population of the ith level and ρ i j are the coherences (i = j). Any initial stateρ(0), which is a linear combination ofÎ andT 8 (0), is therefore also a dynamical invariant. This form incorporates the system state during the STIRAP procedure, 013064-6 which maintains this form throughout the dynamics. In such a process, the Rabi frequencies α(t ) and β(t ) determine the boundary conditions. The population transfer, in the adiabatic regime, is obtained by the following protocol: At initial time, α = 0 and β > 0, implying that only the first state is populated,ρ(0) = |1 1|, Eq. (28). This is manifested by T (α = 0, β > 0) ∝ ρ 11 as shown in Eq. (27). During intermediate times α, β > 0, leading to generation of coherences between states |1 and |3 and rise in population of the third level. This can be witnessed by the nonvanishing terms, proportionate to ρ 13 and ρ 31 and ρ 33 , in T 8 , and nonvanishing off-diagonal terms inρ(t ). At the final time t = T , α > 0 and β = 0, completing a transition of the system towardρ(T ) = |3 3|. The form of adiabatic invariant, Eq. (28), rationalizes the counterintuitive order of pulses of the STIRAP driving protocol [34].
The STIRAP technique, described above, is based on the conservation ofT 8 (t ), which applies only in the adiabatic limit (χ 2 → 0). In the following, we show that the STIRAP can be generalized to the inertial regime (only requiringχ 1 ,χ 2 → 0). This technique is based on general dynamical invariants, which incorporate the adiabatic invariants as a special case.
Diagonlizing M leads to two independent dynamical invariants along with the identity operator [64] andD, given in Appendix D. For slowly varying χ 1 and χ 2 the inertial solution applies, implying that Ĉ (t ) and D (t ) are constant.
The operatorĈ serves as a generalization ofT 8 and converges to it in the adiabatic limit (χ 2 → 0)Ĉ. We can employ this property to construct an inertial STIRAP by setting the same boundary condition as the adiabatic process while uplifting the restriction of slow driving at intermediate times. In this procedure, the adiabatic parameter χ 2 can be large as long as the change in χ 2 is sufficiently slow.
The inertial STIRAP is achieved by the following protocol: At initial time, the Rabi frequencies are set as α(0) = 0, β(0) > 0, χ 2 (0) = 0, and χ 1 = 0, and the system is initialized in the |1 state. Under inertial driving, these initial conditions imply that the system remains in the form throughout the dynamics. From an initial stationary state (χ 2 = 0), the driving is accelerated, leading to α, β > 0 and χ 2 ∼ 1 at transient times and decelerated at the final stage, achieving α > 0, β = 0 and χ 2 = 0 at the final time. Such a protocol transfers population between states |1 and |3 . The inertial STIRAP is expected to share the same robustness as the adiabatic STIRAP. As a simple demonstration, we consider a delta-correlated noise in timing of the driving. Such a process is equivalent to adding random noise to the generalized Rabi frequency (t ), Eq. (26) [65,66]. The effective equation of motion becomes where n is proportional to the noise amplitude. In this case, the noise has no effect on the eigenoperators with vanishing eigenvalues (the time-dependent constants of motion [67]). The dynamics of the transient eigenoperatorsF k are accompanied by an additional decay at rate 2 n λ 2 k , while the phase remains unaffected.
The SU(3) framework can be employed to generalize this scheme, by constructing inertial STIRAP protocols for an Nlevel Hamiltonian. This can be achieved by utilizing the same techniques used to generalize the adiabatic STIRAP [34,68].

E. A fully connected three-level system
A similar framework of the SU(3) algebra is employed to describe the dynamics of a fully connected three-level system, Fig. 2. We study a system with the time-dependent Hamiltonian, where all the Hamiltonian parameters are time dependent. We use the transformation generated by Eq. We consider an inertial protocol varying χ and ξ slowly, with an initial state including a superposition of two eigenvectors in Liouville space:ρ(0) = 1 3Î +v 3 (0) +v 6 (0) (v i is the eigenoperator of λ i ). The trajectory is shown in Fig. 5, superimposed on the eigenvlaue surfaces. Figure 6 compares the inertial solution to a numerical integration of the equations of motion. For Hamiltonian Eq. (32), the geometrical phase in Liouville space is identically zero since the Berry connection vanishes, cf. Appendix D.
The demonstrated inertial protocol can be utilized in quantum control, extending the adiabatic protocols studied in Ref. [69], for a landscape of conical intersections.

V. EIGENOPERATORS EXISTENCE PROOF
The form of Eq. (3) seems to restrict the inertial theorem to specific Hamiltonian dynamics. This is not the case, as we will show in the following proof. We claim that for any HamiltonianĤ (t ) with time-dependent analytical parameters, one can find an orthonormal basis in Liouville space (a set of time-dependent orthonormal operators) such that the generator of the dynamics in Liouville space can be decomposed as Proof. Let {|ϕ } be an orthogonal basis of the system's N-dimensional Hilbert space. Then the set of matricesX nm = |ϕ n ϕ m | form a basis for the associated N 2 -dimensional Liouville space. The dynamics of the Liouville vector X = {X 1,1 ,X 2,1 , ...,X N−1,N ,X N,N } is generated by the Heisenberg . When the Hamiltonian is contained in the closed algebra, formed by {X }, the dynamics of X can be expressed as Next, we define a basis of Liouville space given by where M(t ) = −iV (dV † /dt ) + iVGV † . The desired decomposition is now obtained by choosing a suitable basis transformation of the form V = PWOP −1 , where P is an arbitrary constant unitary matrix and W (t ) is determined by the differential equation dW/dt = −iM(t )W, with initial condition W (0) = I, the identity in Liouville space. The matrix O is chosen to be a diagonal (in the {X } basis) time-dependent real matrix. Substituting V into Eq. (34) leads to the desired decomposition M(t ) = PD(t )P −1 .
To conclude, this result shows that for an analytic Hamiltonian, contained within the operator algebra, by choosing a suitable time-dependent operator basis, all the timedependence of the generator can be absorbed into D and the time-dependent basis. Once such a basis is found, the inertial theorem can be used to generalize this specific solution to a broad family of inertial solutions.

VI. COMBINING INERTIAL SOLUTIONS AND CONSTRUCTION OF INERTIAL HAMILTONIANS
Constructing inertial solutions relies on the intuition related to the system's operator algebra. As a result, the solution may appear restrictive. This disadvantage can be overcome by combining inertial solutions to generate a family of inertial protocols for a larger algebra. This construction then generates a new inertial Hamiltonian. The method is based on the insight that deriving the Hamiltonian from a known propagator is relatively simple (in contrast, the reverse procedure, for a time-dependent Hamiltonian, is extremely involved).
We assume a set of known inertial solutions in Liouville space for the set of Hamiltonians {Ĥ (i) (t )}. The solutions are defined by the Liouville propagators {U (i) (t, 0)}; these determine the dynamics of the operators basis set: v H (t ) = U (i) (t, 0) v H (0). Each Liouville propagator has a matching propagator in the standard Hilbert space of wave functionŝ U (i) (t, 0). In turn, each propagator depends on rapid eigenvalues {λ (i) (i) } and slowly varying inertial parameters {χ (i) }. By taking a product of inertial propagators, we obtain a new inertial propagator, which describes the dynamics accurately when all the inertial parameters change slowly. The new inertial Hamiltonian is determined by the Shrödinger equation where the time dependence in the Hamiltonians and propagators is emitted to simplify the presentation. The form of Eq. (36) defines a broad family of inertial Hamiltonians; in each one, the rapid degrees of freedom { } can be set as arbitrary time-dependent continuous functions, while the inertial parameters satisfy inertial condition.

013064-8
This procedure allows deriving inertial Hamiltonians for arbitrary dimensional systems. Moreover, since, Eq. (35) defines an inertial solution, the procedure can be repeated iteratively increasing the controlable degrees of freedom in the final Hamiltonian.

Explicit demonstration
To demonstrate the method, we combine two inertial solutions of the two-level system, Sec. IV B, and derive a three-level inertial Hamiltonian. The first step requires mapping the Liouville propagator to the associated Hilbert space propagator. The transformation is performed by utilizing the vec-ing procedure, which maps the Hilbert space dynamics to Liouville space [35,37]. The density operator maps to a vector in Liouville space r by ordering the columns ofρ, i.e., for an N by N density operator the (i, j) matrix entry ofρ maps to the ( j − 1)N + i entry of r. Operation of operators from both left and right onρ, asÂρB, maps to the operation (B T ⊗Â) r in Liouville space. Therefore, the time evolutionÛ (t, 0)ρ(0)Û † (t, 0) leads to the Liouville propagator U (t, 0) =Û * (t, 0) ⊗Û (t, 0), where * designates the complex conjugation operation.
In the eigenstates representation {|φ k (t ) }, the wavefunction propagator can be expressed aŝ where {α k } are the associated phases. Transforming the dynamics according to the vec-ing procedure in the eigenbasis representation ofÛ (t, 0) leads to U (t, 0) = diag( N k,l e i(α k (t )−α l (t )) ), while r corresponds to the Liouville state vector in the time-dependent operator basis {|φ k (t ) φ l (t )|}. The entries of r(t ) are given by the matrix element ofρ(t ) in the {|φ k (t ) } basis.
In Liouville space, the propagator is diagonal in the eigenoperator basis, which implies that the set of eigenoperators correspond to the set {|φ k (t ) φ l (t )|}. This identification allows deriving the eigenstates ofÛ (t, 0) and the matching phases {α k (t )}, which determine the propagator, Eq. (37).
We demonstrate the reversed mapping by deriving the Hilbert space propagator for the inertial solution of the twolevel system, Sec. IV B. By following a similar methodology, any Liouville propagator can be mapped to the corresponding Hilbert space propagator. First, the eigenoperators of the propagator {F 1 ,F 2 ,F 3 } are obtained, diagonalizing the generatorB, Eq. (23). Formally, the eigenoperators are the operators associated with the rows of the diagonalizing matrix P −1 , satisfying PBP −1 = diag(0,κ, −κ ). We next normalize these Liouville vectors to obtain an orthonormal set of corresponding operators {Ĝ 1 ,Ĝ 2 ,Ĝ 3 }. These are characterized by simple dynamics (Heisenberg picture): Similarly, the expectation value ofĜ 3 (t ) leads to the conjugate equation of Eq. (38). These identities imply the explicit form of the Hilbert space propagator: (39) For different inertial solutions, such as the harmonic oscillator or three-level system Secs. IV A and IV C, the Hilbert space propagator is obtained by following a similar procedure. First the eigenoperators are divided to pairs, each constitute adjoint pairs. Taking a product of the two operators (two options existĜ iĜ j andĜ jĜi ) gives a projection operator of an eigenstate |φ k (t ) . The eigenstates can then be obtained via diagonalization. The two corresponding eigenstates lead to an analogous equation as Eq. (38), which determines the propagator phases α k (t ) and subsequently, from Eq. (37), the Hilbert space propagator.
We have now reached the point where we can combine two TLS propagatorsÛ (1) (t, 0) andÛ (2) (t, 0), Eq. (39), to obtain an inertial Hamiltonian for a three-level systemÛ (t, 0) = U (1) (t, 0)Û (2) (t, 0). We consider two propagators that couple two energy levels with a single common energy state, which is taken to be the second energy state. Each propagator is generated by an Hamiltonian in the form of Eq. (21), with distinct frequencies ε (i) (t ), ω (i) (t ), and inertial parameters μ (i) (t ), with i = 1, 2. Substituting the associated propagators, Eq. (39), into Eq. (36) leads to a new inertial Hamiltonian where the correction term iŝ and [Û (1) ,Ĥ (2) ] is given in Appendix E. The norm of the correction HamiltonianĤ 12 depends on the control parametersÛ (1) andÛ (2) . We compared the norm of the correction to the norms ofĤ (1) andĤ (2) and find that for typical parameter values the norm is at least an order of magnitude smaller (|μ| ∼ 1 and (i) ∼ 10). This suggests a perturbative approach with respect to this correction.

VII. EXTENDING THE INERTIAL THEOREM TO OPEN-SYSTEM DYNAMICS
The inertial solution describes the free propagation of isolated systems. In reality, no system is truly isolated; as a consequence, the environment modifies the system dynamics. By combining the inertial theorem and the nonadiabatic master equation (NAME) [16] a reduced description of the system dynamics can be obtained, where the influence of the bath is treated implicitly.
The crucial step of the derivation includes solving the free propagation, which in turn is used to obtain the system-bath interaction HamiltonianĤ I in the interaction representation. Applying the inertial theorem to expandĤ I in terms of the eigenoperatorsF k , Eq. (6) and following the derivation presented in Ref. [16] leads to a NAME incorporating the effect of a slowly accelerated drive, d χ/dt = dμ/dt > 0 (see Appendix C).
The Master equation in the interaction representation is given: Here,ρ S (t ) is the system's density operator in the interaction representation relative to the free evolution and F j ≡F j (0). The termH LS (t ) is the time-dependent Lambtype shift Hamiltonian. This master equation, Eq. (42), is an explicit time-dependent version of the Markovian Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) Master equation [70,71]. Within the derivation of Eq. (42), the inertial theorem eigenoperators,F k , Eq. (6), are identified as the jump operators of the master equation. These determine the instantaneous attractor of the dynamical map and the decay rates [16]. The decay rates γ (α k ) are related to the Fourier transform of the bath correlation functions with effective frequencies α k (t ). These effective frequencies are the derivative of the accumulated phases, associated with the eigenvalues ofF k . In Appendix C, the construction of Eq. (42) is demonstrated for a driven system weakly coupled to a bath.
The framework of the inertial NAME, Eq. (42), has been employed in an open system control problem: accelerating thermalization [72,73]. In addition, the theory allowed constructing a fully quantum Carnot analog engine [74,75], highlighting quantum signatures in quantum thermodynamics [76,77].

VIII. GEOMETRIC PHASE IN LIOUVILLE SPACE
In 1984, Berry showed that a system transported adiabatically by varying parameters of the Hamiltonian around a circuit acquires an additional geometric phase [78]. Following a similar proof, we show that the operatorF k (θ ) attains a new geometric phase, φ k , when the parameters {χ } are transformed slowly in circuit C in parameter space, cf. Appendix A. The geometric phase has the form where The geometric phase in Liouville space, Eq. (43), has a different physical significance compared to the Berry phase. The Berry's phase is an accumulated phase of the wave function and therefore is nonvanishing only for a closed circuits including a degeneracy of eigenenergies. As a property of the wave function, it can only be witnessed by interference.
In contrast, the geometric phase in Liouville space influences the physical observables directly. Such observables are determined by a linear combination of the eigenoperators, [associated with the vector v H , Eq. (6)] and the initial system state. Moreover, unlike the Berry phase, the eigenvectors F k are uniquely defined by χ . As a result, φ k is nonvanishing for open circuits in the parameter space {χ }. For a closed circuit, φ k vanishes when V k ( χ ) is analytic within the area encompassed by the circuit. When the circuit surrounds a singularity, which may occur in the case of degeneracies n λ n = k λ k , the geometric phase is nonvanishing and will directly affect the physical state.

Geometric phase examples in Liouville space
As a first demonstration, we consider a two-level system in a time-dependent magnetic field. This system is represented by the Hamiltonian where B i (t ) = f (t )b i (t ) are the components of the magnetic field B(t ) and f (t ) = | B| is a time-dependent function. In terms of the spin operator basis {Ŝ x ,Ŝ y ,Ŝ z }, the dynamics in Liouville space are generated by For such a Hamiltonian, the inertial limit coincides with the adiabatic limit under the scaling t → t 0 dt f (t ). Nevertheless, this example is instructive, as it demonstrates the properties of the geometric phase in Liouville space and highlights the distinctions relative to the Berry phase in the wave-function Hilbert space.
For simplicity, the magnetic field is varied while keeping the magnitude of b(t ) = {b x (t ), b y (t ), b z (t )} T constant and equal to unity. For such a protocol, it is natural to express the dynamics in terms of the spherical angles ( b(t ) = b(θ (t ), ϕ(t ))), leading to the eigenoperatorsF k (θ, ϕ) and associated geometric phases in Liouville space φ k , see Appendix F. A rotation of b(t ) toward directionn(θ f , ϕ f ) leads to accumulated geometric phases: φ ± = ± cos (θ f ) φ, where ϕ = ϕ f − ϕ i is the change in the azimuthal angle, associated with eigenoperatorsF ± correspondingly. These induce a direct effect on the system state.
In Fig. 7, we compare the inertial solution, including the geometric phase, to a solution where the phase is omitted and an converged numerical solution. The result demonstrates the importance of the geometric phase for an accurate dynamical description. Moreover, unlike the Berry phases in Hilbert space, φ k influence the dynamics for closed as well as open circuits in the (B x , B y , B z ) space. An addition example demonstrating the geometric phase employs the Hamiltonian of a fully driven three-level system, where all the Hamiltonian parameters can be time dependent and = α 2 + β 2 . The inertial solution for Hamiltonian Eq. (47) is derived in a similar manner as in Sec. IV C. The time-dependent transformationŜ, Eq. (26), leads to the dynamical operator basis in Liouville space and the M matrix. The associated inertial variables are χ = (αβ −αβ )/ 3 , γ = / , η 1 , η 2 and ξ = δ/ .
An excursion on the eigenvalue manifold is generated by varying the parameters η 1 and η 2 . Figure 8 compares a numerical solution of the dynamics with the inertial solution, with and without the geometric phase. We take the initial state aŝ ρ(0) = 1 3Î +v 1 (0) +v 8 (0) and vary η 1 and η 2 cyclically. The errors of the complete inertial solution arise from deviations from the inertial limit.
The overall phase of the inertial solution is of importance in the derivation of open system dynamics cf. Sec. VII. The phase enters into the detailed balance expression and thus is influenced by the geometric phase. As a result, the kinetic rates of the master equation are dependent on the geometric phase.

IX. DISCUSSION
To summarize, the inertial solutions are constructed in Liouville space, employing a time-dependent operator basis. The transformation to this basis incorporates a part of the time dependence while an additional fast timescale is absorbed into the eigenvalues of the propagator { (t )λ( χ )}. This leads to the required decomposition, allowing us to bypass the time-ordering obstacle. The transformation is the key toward As expected, the distance between the converged solution and the solution not including the geometric phase increases. In contrast, the complete inertial solution remains accurate, achieving distances close to unity. All the system observables show a typical behavior as presented in (b). The observable of the complete inertial solution remains close to the numerical solution, while the solution lacking the geometric phase deviates from the expected value as the phase increases. Gaps between the inertial and numerical solution arise from the small deviation from the inertial limit. The model parameters are = 1/ √ 2, = − √ 2, χ = 0.75, ξ = 0.5, α(t ) = − sin ( χt ), β(t ) = cos ( χt ), η 1 = cos (g(t )), and η 2 = sin (g(t )) with g(t ) = (7π/4)(t/t f ) + π/4 for a final time t f a.u. For a complete circle in parameter space, the total accumulated phase is φ tot ≈ 1.04π . obtaining a suitable decomposition, leading to the inertial solution. An appropriate time-dependent operator basis always exists (cf. Sec. V) but its construction requires ingenuity. The symmetries of Lie algebras can serve as a guideline to guess an appropriate transformation to the time-dependent basis.
The solution is characterized by a set of slowly varying inertial variables {χ } and expressed in terms of the eigenoperators {F }, along with dynamical and geometric phases {φ}. The eigenoperators with vanishing eigenvalues correspond to time-dependent constants of motion [17,18], meaning that the expectation values of these operators are invariant under the driven dynamics. In contrast, the eigenoperators with nonvanishing eigenvalues carry a time-dependent phase.
The inertial solutions share common features with the quantum adiabatic solutions. These are a consequence of the quantum adiabatic theorem, which was derived 90 years ago by Born and Fock [79]. It states that for a slowly varying HamiltonianĤ (t ), an eigenstate |n(0) of the initial HamiltonianĤ (0), remains an eigenstate |n(t ) of the instantaneous HamiltonianĤ (t ) throughout the process. By inserting the instantaneous eigenstate solution into the timedependent Schrödinger equation, the validity of the adiabatic approximation is be determined.
Comparing the structures of the two solutions, the eigenstates of the instantaneous Hamiltonian in the adiabatic solution are replaced by the eigenoperators of the propagator in the inertial solution. The dynamics in both cases includes a geometric and dynamical phases, which differ by their physical meaning. Moreover, the validity of the inertial solution is verified in a similar fashion, by substituting the solution into the dynamical equation in Liouville space. In the special case, where a single transformation generates the time-dependent dynamical operator basis from a time-independent basis, the inertial solution can be formulated as an adiabatic solution in the interaction representation in Hilbert space.
The differences between the adiabatic and inertial solutions are highlighted by applying both methods to the same driven system. This is quantified by the inertial and adiabatic parameters, ϒ and μ. The adiabatic parameter restricts the range of validity of the adiabatic approximation, while the inertial parameter limits the inertial solutions. Typically, the inertial parameter is associated with the change in the adiabatic parameter, and therefore, the range of valid protocols is significantly enhanced to include rapid changes associated with a large μ.
Corrections to the adiabatic approximations have been developed. There is some confusion in their terminology. One approach is based on adding a counterdiabatic control to the Hamiltonian, which eliminates transition terms between adiabatic states, thus, maintaining adiabatic evolution [80][81][82]. A different approach termed superadiabatic is based on a perturbative treatment in orders of the adiabatic parameter [81,83]. An alternative is to employ the adiabatic transformation iteratively. For the SU(2) case, this can be carried out to infinite order [84,85]. This approach could be hard to generalize for a larger algebra with more than one adiabatic parameter and is closest in spirit to the inertial procedure.
Generally, similar closed-form solutions serve as a platform to construct control protocols, for example, the adiabatic STIRAP procedure, which has become extremely popular in contemporary physics [34]. As closed-form solutions, the inertial solutions also generate a constructive family of control strategies. As an example, we introduced the inertial STIRAP, which can be incorporated in similar control tasks as the adiabatic STIRAP but allows for faster control.
An additional possible control strategy concerns the geometric and dynamical phases of the inertial solution. These serve as a new template for interference. Overall, since the inertial solution is stable [51], the control protocols, based on the solution, are expected to be robust in the presence of noise [86]. An analysis of the noise has been performed in Eq. (31) for the three-level system. The generalization to other algebras is straightforward. Similar treatment concerning the robustness were carried out for the Landau-Zener scenario [87,88].
The scope of inertial solutions can be extended by combining individual inertial solutions, Sec. VI. For example, in a multilevel Hilbert space inertial solutions with a common levels can be combined, generating a new time-dependent control Hamiltonian.
The inertial solution also has a direct application for opensystem control. Since any quantum system interacts with the environment to some extent, complete control includes taking into account the affect of the environment. In contrast to the typical analysis, the driving has a direct influence on the system-bath interactions [16]. In the weak system-bath coupling limit, the influence of the free dynamics can be incorporated by the NAME, Sec. VII. In such open system control scenarios, the controller affects the system directly through the driving and indirectly through the system-bath coupling. This property enables extending the common unitary quantum control to transformations involving changes in entropy [72].

APPENDIX A: GEOMETRIC PHASE
We derive the geometric phase in Liouville space, assuming a general non-Hermitian generator M [Eq. (3) in the main text]. The derivation follows the original derivation of Berry [78], extending the solution to a non-Hermitian generator. When v H (χ (t )) completes a contour C in the parametric of {χ } (not necessarily a closed), the inertial solution acquires a geometric phase of the form When the matrix M includes three inertial variables χ = {χ 1 , χ 2 , χ 3 } T , the calculation of the geometric phase is simplified by utilizing common vector calculus identities and Stoke's theorem. Following Berry's derivation [78] and identities Eqs. (11) and (12) in the main text lead to the final result: (A2)

APPENDIX B: COMPARISON OF THE INERTIAL AND ADIABATIC SOLUTIONS FOR SLOW DRIVING
We compare the inertial and adiabatic solutions to converged numerical results for slow driving. A linear ramp protocol is considered, for the harmonic oscillator model, the oscillator frequency increases linearly in time, with δ = (ω(t f ) − ω(0))/t f . The two-level system is modified by a similar protocol, varying the Rabi frequency linearly, (t ) =¯ (0) +δ · t. The comparison between the two solutions is presented in Fig. 9, demonstrating the superiority of the inertial approximation over the adiabatic approximation, in both the adiabatic and nonadiabatic regimes.

APPENDIX C: NONADIABATIC MASTER EQUATION
We present a derivation of a master equation for a driven quantum system interacting with a thermal electromagnetic field with temperature T , see Ref. [16] for further analysis.
The composite system is represented by the Hamiltonian whereĤ (t ) is the driven system Hamiltonian, the bath Hamiltonian is composed of all the bath modes of the formĤ B = λ=1,2 k ω kb † λ ( k)b λ ( k), andĤ I is the system bath interaction term. The interaction term under the dipole approximation can be written asĤ I = E · D, where D is the system dipole operator and E is the electromagnetic field operator. Such a field operator obtains the form E = where V is the volume of the field, ε 0 is the electric constant,b λ ( k) andb † λ ( k) are the annihilation and creation operators of a bath mode in the kth direction with a frequency ω k , (k ≡ | k|), and polarization λ.
Following the microscopic derivation [89][90][91], we transform to the interaction picture relative to the free Hamiltonian H (t ) +Ĥ B . We assume the conditions are such that the inertial approximation is valid. In this regime, the dipole operator in the interaction representation can be decomposed in terms of the time-independent eigenoperators {F }, Eqs. (4) and (6) in the main text, asD where n (t ) is given by Here,F n ≡F n (0), a n = tr( D(0)F † n ) and an upscript tilde designates operators in the interaction picture. Utilizing Eq. (C2), the composite Hamiltonian in the interaction picture can be written as (C5) whereρ B is the density operator of the bath. Assuming the bath correlation functions decay fast relative to the external driving we approximate k (t − s) as where . This approximation is justified, as the bath correlation functions decay in a typical timescale which is much smaller than the timescale of the change in the system parameters, namely, the function k (t ). Thus, the contribution to the integral in Eq. (C5) vanishes when the approximation Eq. (C6) is not satisfied, see Ref. [16] for further details.
Gathering Eqs. (C4)-(C6) leads to with the spectral correlation tensor given by We assume i (t ) + j (t ) 1 for i (t ) = − j (t ), and by performing the secular approximation terminate terms in Eq. (C7) which oscillate rapidly. Furthermore, by following a similar derivation as presented in Ref. [89] Part II, Sec. 3.4.1, the spectral correlation tensor i j can be calculated and written as a sum of two terms j (α) ≡ i j = δ i j ( 1 2 γ j (α) + iS j (α)), with γ j (α) = α 3 |a j d| 2 12π 2h ε 0 c 3 (1 + N (α)) , and S j (α) = |a j d| 2 6π 2h ε 0 c 3 P ∞ 0 dω k ω 3 Here, c is the speed of light, P designates the Cauchy principle part, and N (α) is the occupation number of the Bose-Einstein distribution at frequency α.
The final form of the NAME in the interaction picture can be written as d dtρ  k corresponds to the eigenoperatorF k , which is obtained by summing over the product of the coefficients and the basis operators. For F k = { f 1 j , f 2 j , f 3 j } T ,F j = f 1 jĤ + f 2 jL + f 3 jĈ . The eigenvectors that correspond to the eigenoperators of the bottom block 2 × 2 matrix are with the eigenvalues λ + = κ 2 and λ − = − κ 2 . The dynamics of the eigenoperators has an additional scale ω(t ) ω(0) associated with the diagonal of terms of B. For example, the dynamics of the operator associated with F 1 (vanishing eigenvalue) iŝ F H 1 (t ) = ω(t ) ω0F 1 (0), where the upper-script H designates that the operator is in the Heisenberg picture. Explicitly, the dynamics becomeŝ (0)) .
The geometric phase [Eq. (7) in the main text] is obtained by integrating over the Berry connection: i( G k , ∇ χ F k ) = i( G k , ∂ μ F k ). This expression vanishes for the eigenvectors of the harmonic oscillator, therefore, the inertial solution does not contain a geometric phase.

Two-level system
The eigenvectors that correspond to the eigenoperators and eigenvalues of the propagator [Eq. (23) where the invariant isĜ 1 (t ), given in Eq. (E1). In a similar fashion as in the harmonic oscillator model, the Berry connection, associated with the various eigenvectors, vanishes.