Adiabatic landscape and optimal paths in ergodic systems

Whether one is interested in quantum state preparation or in the design of efficient heat engines, adiabatic (reversible) transformations play a pivotal role in minimizing computational complexity and energy losses. Understanding the structure of these transformations and identifying the systems for which such transformations can be performed efficiently and quickly is therefore of primary importance. In this paper we focus on finding optimal paths in the space of couplings controlling the system's Hamiltonian. More specifically, starting from a local Hamiltonian we analyze directions in the space of couplings along which adiabatic transformations can be accurately generated by local operators, which are both realizable in experiments and easy to simulate numerically. We consider a non-integrable 1D Ising model parametrized by two independent couplings, corresponding to longitudinal and transverse magnetic fields. We find regions in the space of couplings characterized by a very strong anisotropy of the variational adiabatic gauge potential (AGP), generating the adiabatic transformations, which allows us to define optimal adiabatic paths. We find that these paths generally terminate at singular points characterized by extensive degeneracies in the energy spectrum, splitting the parameter space into adiabatically disconnected regions. The anisotropy follows from singularities in the AGP, and we identify special robust weakly-thermalizing and non-absorbing many-body"dark"states which are annihilated by the singular part of the AGP and show that their existence extends deep into the ergodic regime.


I. INTRODUCTION
With the rapid progress of quantum technologies, the design of efficient protocols to control and numerical methods to describe quantum systems quickly moved to the forefront of current research. To achieve a better performance, a crucial element is the ability to perform adiabatic transformations, i.e. transformations between states that are adiabatically connected [1]. For example, quantum annealing and adiabatic quantum computation are based on an adiabatic process transforming a simple initial ground state into a final non-trivial eigenstate, and were shown to be a universal tool for quantum computation [2]. Likewise, any quantum gate operation can be designed using an adiabatic protocol [3,4]. The experimental preparation of equilibrium states in isolated or nearly-isolated systems such as cold atoms or NV centers is often achieved by adiabatic transformations of the Hamiltonian, starting from a simple initial state. In some cases, including Floquet-engineered systems [5][6][7], such a procedure is not only convenient but is actually required, since these systems do not naturally thermalize by interacting with their environment. In the context of thermodynamics, adiabatic (reversible) processes are also of crucial importance. They allow one to minimize the dissipative losses associated with an increase of entropy and achieve the maximal possible efficiency of energy conversion, e.g. in heat engines and refrigerators [8,9]. On the theoretical side, adiabatic transformations underly many concepts including the Schrieffer-Wolff transforma-tion [10][11][12], and the dressing of quasiparticles by interactions underlying e.g. Fermi liquid theory [13]. Such transformations not only allow us to theoretically understand the properties of low-energy Hamiltonians, but also provide a convenient tool to greatly improve the efficiency of numerical methods, allowing one to focus on particular subspaces of interest [12].
A standard limitation of our ability to use adiabatic transformations is that they, almost by definition, have to be extremely slow. In many-body interacting systems, unless we are interested in the ground state of a gapped system, the necessary time scales are exponentially large with the system size [1,14]. A similar exponential slowing down is required if we are interested in following a ground state which either crosses a first-order phase transition [15,16], enters a quantum glass regime [17,18], or follows from an annealing protocol solving a hard computational problem [19].
From the computational point of view, strict upper bounds on the rate of parameter change result in heavy numerical costs. On the experimental side, they lead to a very slow state preparation and large energy processing times. Moreover, the necessary long time scales are generally inaccessible in experimental setups. Systems cannot be perfectly isolated from their environment, leading to decoherence and noise which can destroy the state or erase the information that adiabaticity is trying to preserve. Rather recently, it was realized that this problem can be circumvented and adiabatic transformations can be sped up, in principle arbitrarily, by adding an additional term to the Hamiltonian, suppressing all dynamical/diabatic transitions. Such ideas were first introduced in 2003 by M. Demirplak and S. Rice [20] and independently in 2009 by M. Berry [21] and were subsequently termed counterdiabatic (CD) or transitionless driving. The topic of counterdiabatic driving and the related field of shortcuts to adiabaticity has recently gained tremendous attention in both experimental and theoretical literature [22][23][24][25][26][27][28][29][30][31][32].
In counterdiabatic protocols one applies an additional term to the Hamiltonian, proportional to the generator of adiabatic transformations, the so-called adiabatic gauge potential (AGP). This extra term suppresses all diabatic(non-adiabatic) excitations/losses. The main difficulty of this approach is that the AGP is generally highly non-local. Furthermore, the AGP is not only useful in counterdiabatic driving, but also contains a wealth of information on the geometry of eigenstates and diabatic response [1] and serves as a very sensitive probe of quantum chaos [33]. The exact AGP is local only in some special situations, including symmetry transformations or transformations of the ground state of a gapped system [34][35][36]. Fortunately, even if the exact AGP is generally out of reach, it was recently realized that in some specific instances we can find an approximate yet accurate local AGP using a variational minimization [1,26,37]. The resulting local AGP was shown to be highly efficient both in solving computationally difficult problems [38,39] and performing efficient Schrieffer-Wolff transformations [12,40]. Still, several general and unanswered questions remain: (i) When do such local approximations apply? (ii) Which are the optimal protocols for local adiabatic evolution? (iii) Can we learn which states are most heavily affected by diabatic effects from these local approximations and is it possible to identify the states for which dissipation is minimal?
In this work, we first focus on finding an optimal path in the space of system's parameters to design local protocols for adiabatic state preparation. Very often, physical systems are controlled by multiple parameters, e.g. pressure, temperature, chemical potential, external electric and magnetic fields in thermodynamics or single-spin controls and two-spin interactions in quantum control. While the order in which these parameters are changed will not matter if everything happens perfectly adiabatically, the diabatic effects can vary drastically depending on how these parameters are tuned. It is then natural to ask for the optimal path in the space of parameters, minimizing diabatic effects. This will be the focus of this work.
For concreteness, we will consider protocols satisfying the time-dependent Schrödinger equation (we set = 1 throughout the text), where the Hamiltonian depends on a set of timedependent control parameters λ(t) and we initialize the system at t = 0 in a stationary eigenstate of H( λ(0)) (all our results immediately extend to mixed initial states). The question is then how to vary λ(t) such that the state remains close to an instantaneous eigenstate of H( λ(t)).
To answer this question, we analyze the adiabatic landscape of a fairly generic non-integrable 1D Ising model characterized by two independent couplings (cf. Eq. (3)). Specifically, we show that the variational adiabatic gauge potential (VAGP), which gives the best local approximation to the exact AGP (see Eqs. (4) and (5)), forms a two-dimensional vector space, and the directions where the norm of the VAGP is minimal define the optimal paths minimizing diabatic effects. We mainly focus on infinite temperature states, where the equilibrium properties of the system are completely featureless. Nevertheless, the problem of adiabatic continuation remains well defined and highly nontrivial. We find that the evolution along the optimal direction is efficient; that is, eigenstates which are drawn from the middle of the spectrum remain close to the instantaneous eigenstates, maintaining small energy variance. As we will show below (see also Refs. [26,33]) the AGP can be expressed through the long-time limit of non-equal time correlation functions of the operators conjugate to the coupling. Therefore, they cannot be analyzed by the methods of equilibrium statistical mechanics. Our findings thus imply that temperature plays a much smaller role in adiabatic transformations than in equilibrium settings. Let us now introduce the Hamiltonian that we will analyze in this work, describing the quantum Ising model in the presence of a longitudinal and transverse field, as and we introduce a shorthand notation that is convenient for translationally-invariant systems with and so on. Fixing the coupling in front of the Ising interaction ZZ to be unity, J = 1, h and g will be taken as control parameters throughout this paper. The main results of our paper are summarized in Fig. 1: each point represents a choice of couplings defining a Hamiltonian, and the lines show the optimal adiabatic directions presented as a flow diagram. As will be discussed later, this diagram has a very rich structure and is in many respects similar to the standard equilibrium phase diagrams (except that, as already pointed out, it corresponds to an infinite temperature). Let us now summarize the most essential findings reflected in this figure, which will be explained in detail in the paper. , where the optimal direction is approximately the radial one. The norm of the variational gauge potential is highly anisotropic: near the source flows the norm is small and nearly system-size independent along the optimal directions, while increasing drastically in the orthogonal direction and diverging exactly at the points (a) and (c). For the 5-body ansatz an additional singular point appears at (h, g) = (1, 0) [(b)], strongly disrupting the optimal directions in its vicinity.
exponentially-large degeneracies in the energy spectrum, which we call macroscopic degeneracies. These singularities serve as sources/sinks of adiabatic flows and play a similar role to critical points in equilibrium phase diagrams.
• Close to these singularities, the VAGP becomes infinitely anisotropic, with highly-anisotropic regions extending far away from the singularities. This high anisotropy implies that optimal directions, along which local adiabatic transformations are highly efficient, remain well defined. Such optimal directions define paths with minimal dissipation and maximum fidelity for state preparation.
• Near the singular points, there are special manybody "dark" states annihilated by the diverging part of the VAGP. These states exist throughout the entire energy spectrum; similar to the anisotropic regions these are highly robust and extend deep into the ergodic regimes, bearing many parallels with the recently-discovered quantum scars and the eigenstates of constraint models [41][42][43][44][45][46][47][48][49]. While these states form an exponentially-small fraction of the total Hilbert space, their total number can still be exponentially large; as they are immune to the usual dissipation they can be efficiently prepared both numerically and experimentally.
• The optimal adiabatic directions allow us to define adiabatic flows similar to the renormalization group flows, as shown in the figure, and these flows in turn define adiabatically-connected families of Hamiltonians. The norm of the exact AGP is equivalent to the Fubini-Study metric defining the distance between eigenstates of adiabaticallyconnected Hamiltonians [1]. Therefore, these flows can be interpreted as lines approximately minimizing the local distance between eigenstates (or more accurately between energy shells) of different Hamiltonians. Along these flows, both states and operators can be dressed to a very good accuracy under the unitary transformations generated by the local VAGP. In particular, such directions are characterized by the existence of nearly-conserved operators, which are locally-dressed operators conjugate to the coupling along these directions.
• Near the singular points, the VAGP diverges in all directions except for the optimal one. However, the divergent part of the VAGP is a well-defined local operator, implying that a local dressing can be used to efficiently perform adiabatic rotations near these singularities. Combined with the fact that all adiabatic flows terminate at one of such singularities, we arrive at the interesting conclusion that any optimal adiabatic path between two generic points goes through one of these singularities. In other words, the system first has to be brought to the singular point, then a local rotation needs to be performed, before going to the target point along a different flow line. Importantly, such a path can be always found locally by following the optimal direction of the adiabatic flow.
• The optimal directions generally depend on the support/size of the variational ansatz (see the top and bottom halves in Fig. 1), i.e. the support of the operator generating approximate adiabatic transformations. New singularities appear in the higher-order variational ansatz with an increased local support, reflecting higher-order divergences in the perturbative expansion of the AGP. These singularities arise from the degeneracies associated with higher-order interactions and appear at rational couplings, bearing many similarities to the divergences appearing in both KAM theory [50] and locator expansions [51]. The emergence of higherorder singularities indicates that it is not possible to improve local dressing, by either adding additional local terms to the CD protocol or by slowing down the ramping rate in the absence of CD driv-ing, without abruptly altering the path near these new singularities.
• The adiabatic flow diagram remains well-defined even at infinite temperature, where no structure exists in the equilibrium state according to statistical mechanics. Interestingly, many of its features persist at all temperatures, all the way down to the ground state at zero temperature.
We confirm these general findings with numerical simulations for the non-integrable 1D Ising model described by the Hamiltonian (3). Our results can have a broad range of applications in various problems, beyond simply finding optimal paths for annealing or state preparation. In particular, they can be used to find efficient local conservation laws and corresponding "most-integrable" directions, to find the nearest integrable (simple) points that are locally connected to a Hamiltonian of interest, to define most efficient ways of obtaining effective lowenergy theories starting from a noninteracting model, and so on.
This paper is organized as follows: In Sec. II, we introduce the VAGP and define the optimal adiabatic directions. Applications to approximate CD driving and slowest operators are also explained there. Sec. III is the highlight of this paper, where we obtain the flow diagram that defines the optimal directions at each point of the coupling space. We demonstrate that both for conventional adiabatic driving and for the approximate CD protocols state preparation along the optimal paths shows a much better performance than along the orthogonal directions. We explain that the flows terminate/start at special sources/sinks, where the VAGP develops divergencies in the orthogonal directions, becoming infinitely anisotropic, and show how these singularities arise from the perturbative expansion of the exact AGP. We then explain the emergence of special dark states unaffected by the singular part of the VAGP. In Sec. IV, we study how the VAGP depends on the size of the variational ansatz and explain the emergence of new singularities near rational values of h. We then use the VAGP to construct approximate local conserved operators and analyze their life times in Sec. V. Details of the perturbative expansion are given in Sec. VI and Sec. VII is reserved for conclusions.

II. VARIATIONAL ADIABATIC GAUGE POTENTIAL
In this section we will give a brief introduction to the concept of the (variational) adiabatic gauge potential, emphasizing its structure as a vector in a system with multiple controls (tunable parameters). Much of this discussion can be found in earlier papers [1,26,37], but is included here in order to be self-contained and to make an explicit connection of VAGP with slow operators [52,53], operator spreading [54][55][56][57][58][59][60], and emergent conservation laws [61], which will be relevant for the presented flow diagram.

A. Theoretical background
Let us consider a family of Hamiltonians H( λ), where λ specifies the space of available couplings or controls. Any protocol corresponds to a time-dependent choice of λ(t), with an adiabatic protocol corresponding to a vanishing time-derivative |˙ λ(t)|.
The effects of time-dependent couplings are most clearly illustrated in the instantaneous (co-moving) eigenstates of the Hamiltonian |n( λ) , satisfying Any change in the control parameters corresponds to a change in the eigenstates, and one can formally define the adiabatic gauge potential (AGP) as the Hermitian operator A( λ) generating these basis changes [1]: in which ∂ j is the partial derivative w.r.t. λ j . Note that, since eigenstates are only defined up to a phase (or more general rotations in the presence of degeneracies), the AGP is not uniquely defined and supports a gauge freedom.
We will be interested in finding the time evolution (1) of an initial pure state |ψ(t = 0) under time evolution governed by a time-dependent Hamiltonian H( λ(t)), where the only explicit time dependence is through the control parameters [62]. Expanding this state in the comoving basis |ψ(t) = n a n (t)|n( λ(t)) , it is easy to check that the time evolution in this new basis is governed by the moving Hamiltonian Specifically, H nl m (t) = n( λ(t))|H m (t)|l( λ(t)) , which takes the form of a regular matrix representation of the Schrödinger equation, but with time-dependent basis states, which are accounted for by the second term in Eq. (7). In the limit˙ λ → 0 this additional term vanishes such that there are no transitions between instantaneous eigenstates of H( λ). At non-vanishing |˙ λ| the extra term in the moving Hamiltonian, proportional to the AGP, cannot be neglected. Since H( λ) is by construction diagonal in the co-moving frame, all diabatic excitations/losses are generated by the off-diagonal elements of the AGP. Following Ref. [1], Eq. (5) can be recast as an operator equation in which The matrix G j ( A) is diagonal in the eigenbasis of H and its diagonal matrix elements are given by ∂ j n ( λ), the generalized forces conjugate to λ j . In other words, one can view any infinitesimal deformation of the Hamiltonian along the λ j direction ∂ j H as consisting of a spectrum change encoded in G j and an eigenbasis rotation encoded in A j . Eq. (9) remains well-defined in both the classical and thermodynamic limits. However, with the exception of symmetry transformations/integrable systems, the solutions to this equation are generally unstable to infinitesimal perturbations and might not even exist in either of these limits [1,14,33]. Therefore, finding approximate local gauge potentials is essential to circumvent this problem. One goal of this paper is to convey that, even though the exact AGP might be ill-defined, such local approximations can be well-defined and meaningful.
A particularly powerful approach to finding approximate solutions is the variational method. It is based on the observation that Eq. (9) can be interpreted as the minimization condition for the auxiliary action S( A) [37] δS Approximate solutions of Eq. (9) can be found by choosing a specific subset of operators as an ansatz for the AGP and finding the minimum of the action. We call the resulting solution the (local) variational adiabatic gauge potential (VAGP). Also note that the action for the VAGP in Eq. (11) can be interpreted as the action at infinite temperature. In principle, it can be extended to finite temperatures through the introduction of a thermal state exp [−βH] in S (see Ref. [37]), although this strongly complicates the resulting minimization. In this paper we focus on variational manifolds consisting of all local operators with a given support (see Sec. II B for details). One can develop a similar expansion based on nested commutators of ∂ j H and H [26]. We checked that this second expansion leads to very similar conclusions. Despite being an approximate solution, as we discuss below, the local VAGP can be used to determine highly nontrivial properties of the system. Let us mention a few of them.
Approximate counterdiabatic driving. -The notion of counterdiabatic (CD) driving immediately follows from this derivation, since the exact solution of A can be used to completely suppress energy dissipation by evolving a system with the CD Hamiltonian including an additional term˙ λ · A( λ), Representing the evolution in the co-moving frame of H( λ(t)), the additional counterdiabatic term cancels, such that the moving frame Hamiltonian is exactly given by H( λ(t)), which is diagonal and hence does not lead to any excitations or dissipation. Namely, starting from any energy eigenstate |ψ(t = 0) = |n( λ(0)) the state at later times remains an instantaneous eigenstate |n( λ(t)) . In the limit of an infinitely fast rate of change |˙ λ| → ∞, the AGP dominates, and the resulting evolution can be seen as a pure dressing of the initial state. We will refer to a protocol corresponding to H( λ(t)), where no CD term is present, as the unassisted protocol. While the exact AGP generally cannot be realized in many-body systems, the use of local approximations from the variational minimization has already been shown to lead to a significant suppression of transitions [26,[37][38][39]63]. As such, the availability of an accurate local VAGP can also be used to reduce dissipation and design efficient annealing protocols.
Approximate state dressing. -Starting from an initial eigenstate of the instantaneous Hamiltonian, counterdiabatic driving can be interpreted as interpolating between two limits:˙ λ → 0 returns adiabatic state preparation, whereas˙ λ → ∞ dresses the initial state with the (approximate) gauge potential. Namely, in this limit the Schrödinger equation reduces to For an exact AGP, |ψ(t) = |ψ( λ(t)) and this equation reduces to This corresponds to a (quasi-)adiabatic dressing of the initial state [12,34,36,40]. The possibility of such dressing with a (quasi-)local A is a crucial ingredient in classifying topological phases, where all ground states within a given phase can be adiabatically connected using a local dressing.
Operator spreading. -A formal solution to Eq. (9) can be found using the Lehmann's representation as where is the operator conjugate to the parameter λ j , ∂ j H, in the Heisenberg representation w.r.t. the instantaneous Hamiltonian H. For classical Hamiltonian systems, this representation originates from C. Jarzynski [14]. As mentioned before, the exact solution is highly sensitive to the choice of ∂ j H and the limit → 0 will generally diverge in chaotic systems. Keeping finite then corresponds to finding an approximate AGP, which will be local for a local ∂ j H due to the finite support of (∂ j H) (t) at finite times, following recent results on operator spreading (e.g. [58]) and Lieb-Robinson bounds [64]. This representation has also been combined with the variational principle to find an efficient variational ansatz in chaotic many-body systems [26]. Conservation laws and slowest operators. -A local AGP immediately implies an additional local conservation law, since G j ( A) by definition commutes with the Hamiltonian. Minimizing the action then corresponds to obtaining a 'slowest operator' [52], minimizing the commutator with the Hamiltonian (setting the time scale for thermalization), which then becomes an exact conserved quantity if the local VAGP becomes an exact AGP. Interestingly, if we consider the representation of the AGP through Eq. (15) with finite , the corresponding G j ( A) exactly coincides with the approximately-conserved operator obtained by the time-averaging of (∂ j H) (t) introduced in Ref. [61]. In particular, using Eq. (15) it is easy to check that namely, it is the part of ∂ j H that is conserved and does not decay with time.

B. Optimal adiabatic directions
From Eq. (7) it can be seen that all diabatic transitions are induced by the AGP. For a time-dependent change along a certain direction˙ λ = |λ| n λ , for a fixed rate of change |λ| along a direction set by a unit vector n λ , these transitions can be expected to be maximally suppressed along directions where the norm of A λ = n λ · A is minimal. In the same way that the gap between the ground state and the first excited state sets the time scale for quantum annealing, the norm of the AGP along a certain direction sets the scale for the rate of change of the control parameter |λ|: for small ||A λ || the control parameter can be changed rather fast without inducing large diabatic effects, whereas for large ||A λ || even slow deformations of the Hamiltonian immediately lead to diabatic transitions. While the local VAGP is not exact, it contains information about transitions through local interactions, which are often the most damaging because they can lead to a large energy transfer. We will demonstrate below that this is indeed the case.
Given a multi-dimensional space of control parameters, one can thus set the optimal direction as the direction for which the norm of the VAGP is minimal. In principle, one can define different norms, so the minimization procedure is not unique. For example, one could choose norms tailored for particular states, e.g. the ground state. In this work we will use the Fröbenius (L2) tracenorm, equivalent to the common infinite-temperature norm. These norms have the advantage that they can be easily calculated in large systems (including the thermodynamic limit) without any need to diagonalize the Hamiltonian. As such, the actual minimization is particularly straightforward. Remarkably, it was observed in Refs. [26,[37][38][39] that this infinite-temperature norm still provides excellent results even considering e.g. only dissipation from the ground state. Colloquially, using the Fröbenius norm to find the VAGP is similar to optimizing an ice cream recipe inside a very hot oven and then applying this recipe inside a freezer to efficiently prepare the ice cream. Remarkably, this procedure works amazingly well in various systems.
Rather than keeping the discussion maximally general, we will focus on a two-dimensional parameter space with controls set by g and h (see Eq. 3), such that λ = (g, h) and analyze infinitesimal deformations (h + δ cos ϕ, g + δ sin ϕ) with an infinitesimal δ, such that n λ = (cos ϕ, sin ϕ). The generalization of this methodology to more parameters is straightforward. Crucially, the action S defined in Eq. (11) is quadratic in the variational parameters, such that the minimization will give rise to a set of linear equations and the VAGP in an arbitrary direction will be a linear combination of the solutions corresponding to δh (ϕ = 0) and to δg (ϕ = π/2). We can write in which A h and A g minimize the action S h and S g respectively. The Hamiltonian is set by the parameters λ, as denoted in the subscript (where we dropped the vector notation), while the argument denotes the direction in which this Hamiltonian is varied. Defining it can easily be checked that the norm of the VAGP is minimal for ϕ = α ± π/2, α ∈ [−π/4, π/4], and maximal in the orthogonal directions ϕ = α and α + π if ||A g || > ||A h ||, while in the other case the extrema are exchanged (see also Appendix A). We will call these directions optimal and orthogonal respectively. In the following sections, we will analyze the geometric structure of these directions and the resulting anisotropy as a function of (g, h). Note that this also highlights that the directions set by ϕ and ϕ + π are equivalent since they correspond to the same perturbation, only with a different sign (which does not influence the norm of the VAGP).
For translationally-invariant spin-1/2 systems of size L with periodic boundary conditions, like those described by the Hamiltonian (3), we define the k-body operator space H k , k < L, as the zero-momentum space of all operators having support of up to k sites, where we will choose strings of Pauli matrices as basis operators: where the index n stands for the set {s 1 , . . . , s k } and σ s i is one of the Pauli operators {σ x , σ y , σ z , 1} acting on the site i. To avoid double-counting the identity operator is excluded from the right boundary, i.e. s k = 1. We will use a local variational ansatz with a fixed support: We call this the k-body ansatz of the variational calculation, and solve Eq. (11) with the ansatz (21). Since all operators O n are traceless and orthogonal, satisfying Tr(O n O m )/N = Dδ nm , where D = 2 N is the Hilbert space dimension, the minimization problem is straightforward and the solution is formally given by where ad P k HP k A ≡ [P k HP k , A], ad −1 P k HP k is the pseudoinverse of ad P k HP k , and P k is a super-operator which projects an operator onto H k .
In the limit where this operator basis is complete we can consider e.g. projectors on eigenstates as basis operators, which returns the formal solution which can be checked to be equivalent to Eq. (15).

III. ADIABATIC FLOW DIAGRAM OF THE QUANTUM ISING MODEL WITH LOCAL VAGP
In this section we will discuss in detail the flow diagram and the emerging physical implications for a particular, but fairly generic, quantum Ising model, which we introduced earlier in Eq. (3). We will first analyze this diagram using the VAGP obtained within the lowest-order approximation, which already yields non-trivial results. Namely, we will consider a variational manifold with support up to three sites for the VAGP. The motivation for this ansatz is that, as we discuss below, it reproduces the leading-order behavior and the most important singularities of the exact AGP near the strongest macroscopic degeneracy points. These singularities underly several key properties of the adiabatic flows and allow us to reveal the origin of special dark weakly-thermalizing states similar to those found in e.g. Ref. [40]. In the next section, we will then show how the results of this section are affected by adding terms with a larger support into the variational manifold. Before discussing our findings, let us mention a few properties of the Ising model that will be relevant later in the paper.
• There are two integrable lines corresponding to i) g = 0: the so-called classical Ising model with strictly local integrals of motion (z-magnetization for each spin) and ii) h = 0: the transverse field Ising model, which maps to free fermions through the Jordan-Wigner transformation and which has quasi-local integrals of motion constructed from fermion bilinears [65,66]. There is an additional trivially-integrable point corresponding to h 2 + g 2 → ∞, which describes noninteracting spins. Away from these points the model is believed to be chaotic, satisfying the eigenstate thermalization hypothesis (ETH) [67].
• The ground state of the Ising model undergoes a quantum phase transition from an anti-ferromagnet corresponding to small magnetic field to a paramagnet at large magnetic field [68]. On the integrable lines, the critical line separating the two phases terminates at the points (h, g) = (2, 0) and (0, 1). We note that changing the sign of the ZZ coupling moves this phase transition line from the ground state to the most excited state. Therefore, this sign does not affect our "infinite temperature" flow diagram.
• The "classical Ising" line g = 0 additionally contains macroscopic (exponential) degeneracies of the spectrum at any rational value of the longitudinal field h. In particular, at h = 0 and H = ZZ, any configuration with the same number of domain walls has the same energy, e.g. |. . . ↑↑↓ . . . and |. . . ↑↓↓ . . . . At h = 2 and H = ZZ +2Z, any local spin flip from a local "down" to "up" state that creates two domain walls does not change the energy of the system, e.g. |. . . ↓↓↓ . . . and |. . . ↓↑↓ . . . are degenerate. In a similar way, at other rational points of h one can always find many combinations of spin flips leaving the energy of the system invariant. Finally, the h → ∞ point is also macroscopically degenerate: the energy does not change under arbitrary spin flips preserving total magnetization.

A. Flow diagram for the 3-body variational ansatz
As mentioned in Section II, one can systematically define the adiabatic flow diagram by following the directions of the minimal norm of the VAGP. The resulting diagram with respect to the couplings (h, g) as obtained within the 3-body variational ansatz for the VAGP is shown in the bottom half of Fig. 1 as well as in Fig. 2. Note that, on the one hand, the representation of the diagram on a sphere is more natural since all the Hamiltonians with large magnetic field are equivalent to each other up to trivial spin rotation and correspond to the same point in Fig. 1. On the other hand, the "Cartesian" representation shown in Fig. 2 is easier to visualize in the most interesting regime where neither h nor g are too large.
One can observe that the optimal flows form radial patterns centered around singularities at (h, g) = (0, 0) and at (2, 0) (as well as near h 2 + g 2 → ∞ in the spherical representation). Interestingly, these singularities lie at the endpoints of any adiabatic flow: if we start at any generic point (h, g) and follow the optimal adiabatic direction, we will end up in one of these singularities. Likewise, these singular points are good starting points for quantum state preparation in e.g. quantum annealing protocols, because any point of the control space (h, g) can be reached by starting at either of these singularities. At first sight this result seems surprising: these singular points are clearly the points corresponding to large macroscopic degeneracies, where adiabatic transformations are ill-defined. Indeed, our common understanding of adiabatic transformations suggests that one should avoid situations with closing gaps between eigenstates. Thus, naively, one should generally avoid such singular points. As we will show, this reasoning only applies to the orthogonal azimuthal directions, where the norm of the VAGP becomes divergent and strong diabatic effects come into play. However, such divergences remain suppressed in the radial directions. Let us also point out that at the singular points the Hamiltonian splits into a sum of mutually commuting terms, such that its eigenstates are factorisable and thus easy to prepare.
The radial flow near h = 0 implies that the optimal deformation of the Hamiltonian is along the instantaneous magnetic field, (δh, δg) ∝ (h, g). Intuitively, one can understand this result using the domain wall picture: at small magnetic fields one can think about the Ising model as a weakly interacting gas of domain walls separating regions of positive and negative magnetization. The number of the domain walls is conserved by the ZZinteraction. In this manifold of states the Z-magnetic field plays the role of an effective linear potential and the X-magnetic field plays the role of the domain-wall hopping amplitude. The two terms can be combined into an effective non-interacting Hamiltonian describing these domain walls. The radial deformation of h and g then amounts to a simultaneous rescaling of these two parameters of the effective Hamiltonian, which does not induce diabatic transitions between the eigenstates. Similar considerations apply to the other singularity at (2, 0), where the effective Hamiltonian becomes the PXP model [69] with h − 2 playing the role of the potential and g playing the role of the magnetic field. At the third degenerate point, at infinite magnetic field, the radial deformation is trivially the most adiabatic direction, since it simply amounts to rescaling the full Hamiltonian. We emphasize that, while this intuition can generally be justified by considering low-energy effective Hamiltonians, the optimal directions remain well-defined for all eigenstates. We justify this conclusion below by analytically constructing the VAGP near these points, where the radial directions are explicitly shown to be non-singular.

B. State preparation along the optimal flow directions
Before discussing the emergent features of the flow diagram in more detail, let us immediately analyze its implications for quantum state preparation. All calculations and presented diagrams hold at the operator level, so it is natural to first ask if the (operator) flows for the VAGP are representative of similar flows in the context of quantum state preparation, where only a single eigenstate is relevant. Second, if such state preparation is assisted with local counterdiabatic driving using the VAGP, a follow-up question is if the optimal directions using the VAGP are also the ones where the approximate counterdiabatic driving is maximally effective. Here, we will present numerical evidence that suggests a positive answer to these two questions.
Since we are not necessarily interested in the ground state and will consider excited states, a good measure for the proximity of any prepared state |ψ to an eigenstate of the instantaneous Hamiltonian is the energy variance δE: where |ψ is the state prepared according to a protocol following a particular, e.g. optimal, path. If the system is prepared in an exact eigenstate of H( λ), this en-  Fig. 4). System size is L = 12.
In the unassisted protocol, the energy variance is already much smaller along the optimal directions than along the orthogonal directions. Applying local 2-body CD driving, the energy variance drastically reduces even more along the optimal direction, while it only gradually decreases in the orthogonal direction. At infiniteλ the state along the optimal direction can similarly be accurately approximated by a 2-body dressing of the initial state, whereas the accuracy along the orthogonal direction only gradually increases.
ergy variance clearly reduces to zero, whereas a non-zero value indicates how strongly this state has mixed with different-energy eigenstates. We will consider unassisted state preparation protocols and approximate CDD protocols, where the adiabatic evolution is assisted by the strictly local VAGP. In both cases we will compare different paths in control space. For the CDD protocols we numerically solve the Schrödinger equation using the Hamiltonian (12) along a given path λ(t) with A( λ) replaced by its variationally-obtained approximation. The initial state is chosen to be one of the eigenstates of the initial Hamiltonian H near the middle of the spectrum, and we then compute the energy variance at the final value of λ according to Eq. (24). While the results are presented for a single (generic) eigenstate, we checked that these are representative for most eigenstates (exceptions will be discussed in Sec. III D).
All protocols are characterized by the total time duration T , where the limit of large T corresponds to adiabatic evolution, while the limit of small T corresponds to the instantaneous quench for the unassisted protocol and to a dressing of the initial state with the VAGP for the CDD protocol. We choose a smooth protocol to help eliminate diabatic effects at the protocol boundaries [1] interpolating from λ(0) = 0 to λ(T ) = 1, where we set the total protocol duration T = 2 for concreteness, and take (h(t), g(t)) = (h(0), g(0)) + λ(t)(h(T ), g(T )). However, we checked that all the presented results remain qualitatively similar for other time dependences. In Fig. 3 we present the resulting energy variance of the final state for different preparation protocols with the same final Hamiltonian but different initial Hamiltonians, corresponding to different directions of state preparation. For the optimal protocol, the initial point is chosen as (h, g) = (0 + , 0 + ), with a small = 0.01 lifting the degeneracies of the eigenstates, which is then linearly evolved along the radial direction to the final point (h, g) = (0.5, 0.5) (cf. green line in the inset of Fig. 4). This can be contrasted with the state preparation protocol along the orthogonal direction, taking as initial control parameters (1 − , ) and again linearly deforming the Hamiltonian to the same final point (0.5, 0.5) (cf. red line). We checked that starting from another point along the orthogonal direction, namely ( , 1 − ), leads to similar results (cf. Fig. 12).
In order to compare the unassisted protocols, we consider a linear rampλ = 1/T and present the final energy variance for different ramp rates along different directions in Fig. 4. It is clear that the protocol along the optimal direction generally has an energy variance that is orders of magnitude smaller than the energy variance along the orthogonal direction. Even more, when increasing T (nearing adiabaticity), the energy variance for the optimal path decreases much faster, as indicated by the steeper slope in the log-log scale. Interestingly, for evolution along the sub-optimal direction starting at (1 − , ), the energy variance does not decrease in the interval 0.01 ≤ 1/T ≤ 0.1, indicating a complicated landscape of energy level crossings. A similar situation occurs, for example, in Floquet systems [70]. Still, we checked that eventually the energy variance starts decreasing again for 1/T ≤ 0.005.
Using the calculated VAGP for approximate local CDD (see Eq. (12)) to improve on the unassisted protocol, panels (a) and (c) in Fig. 3 show the energy variance for the CDD protocols along the optimal direction with either finite duration T = 2 (a) or infinitely fast T → 0 (c), which effectively corresponds to dressing the initial state with the VAGP. Different colors correspond to a different size of the variational ansatz for the VAGP, with the unassisted protocol included as reference. Panels (b) and (d) show related results for state preparation along the orthogonal direction. Again, it is clear from the plot that the energy variance is generally smaller for state preparation along the optimal direction. Even more, including (approximate) local counterdiabatic terms can be used to drastically reduce the energy variance along the optimal direction. We note that along the optimal direction the VAGP for 1-body ansatz is found to be exactly zero; therefore the results for k = 0 and 1 completely overlap each other. While including the approximate counterdiabatic term along the orthogonal direction also systematically reduces the energy variance with increasing ansatz size, its effect is not as pronounced as along the optimal direction.

C. Asymptotic behavior of the VAGP near singular points
From the structure of adiabatic flows shown in Figs. 1 and 2, it is clear that the points (0, 0) and (2, 0) play a special role, serving as sources/sinks of these flows. As already mentioned, these points also correspond to Hamiltonians with macroscopic (exponential) degeneracies in their energy spectrum. As will be discussed in this section, these points control many important properties of the AGP, including the large anisotropy between optimal and orthogonal directions and the existence of special dark/non-thermal states far from the edges of the spectrum.
In order to understand these properties, we consider perturbative expansions of the exact AGP near these two singular points. The full formalism will be developed in Sec. VI, and here we will focus on the leading-order terms only. Near (h, g) = (0, 0), the dominant term in the perturbative expansion is given by with r = | λ| = g 2 + h 2 . Here the angle θ characterizes the magnetic field in the Hamiltonian H(h, g) through (h, g) = (r cos θ, r sin θ), whereas the angle ϕ characterizes the direction in which this magnetic field is perturbed (δh, δg) ∝ (cos ϕ, sin ϕ).
From Eq. (26) it is clear that the AGP diverges at (0, 0) for a general ϕ. However, along the radial direction ϕ = θ the singular term exactly vanishes, indicating that the radial direction is the optimal one. It is also evident that the anisotropy between the optimal and orthogonal directions diverges near this singularity. This perturbative expansion also highlights that the variational ansatz for the VAGP minimally requires 3-body terms in order to correctly capture the singularity and the corresponding anisotropy. The increasing anisotropy as the magnetic field goes to zero is clearly visible in the 3-body VAGP, as illustrated in Fig. 5. In this plot we show the norm of the 3-body VAGP along the optimal and orthogonal directions as a function of r at a fixed angle θ = arctan(0.2), such that g = 0.2 h. The lines are the fits to the constant (optimal) and 1/r (orthogonal) asymptotes expected from perturbation theory. Interestingly, the perturbative scaling of the norm of VAGP extends up to a relatively large value of the coupling r = 0.4, such that the effects from the singular point can remain important deep into the ergodic regime of the flow diagram. In Appendix B, the individual weights of the terms in the expansion are compared with the scalings from perturbation theory, and it is confirmed that the dominant terms are of the form (26).
The operator divergence can immediately be connected to the eigenstate structure of the Hamiltonian at (0, 0). As already noted, the energy of the model only depends on the number of domain walls, leading to macroscopic degeneracies in the eigenspectrum. The operator Y − ZY Z can be seen as a 'dressed' version of the spin flip operator Y , which however only creates a spin flip if it does not change the number of domain walls, connecting the degenerate eigenstates. This macroscopic degeneracies in H and their splitting by the perturbation effectively dominate the perturbative AGP and lead to well-defined local terms.
A very similar structure emerges near the second singularity (2, 0), where the perturbative expansion of the exact AGP yields (see again Sec. VI) where now r = (h − 2) 2 + g 2 and ϕ is again the angle characterizing the deformation δλ. We introduced the notation P for the projector on the down state of the spin along the z-direction. In the extended notation, the P Y P term reads Same as near the (0, 0)-singularity, the AGP diverges as r → 0 except in the radial direction φ = θ. Therefore the AGP again becomes infinitely anisotropic in the limit r → 0. This singularity is precisely reflected in the flow diagram indicating that the optimal directions are radial. Interestingly, and not accidentally, the leading-order singularity of the AGP is nothing but the generator of spin rotations of the effective low-energy P XP model emerging near the (2, 0) point [69]. This model was already shown to satisfy highly unusual properties, including the existence of weakly thermalizing quantum scar states [69] and the existence of nearby integrable deformations of the Hamiltonian [71]. In the next section, we will show that some (and probably all) unusual properties of this model are encoded in the exact AGP and can be observed in its local variational approximation.
Since it was recently noted that the AGP generates the effective Schrieffer-Wolff Hamiltonian, it is also worthwhile to note that the effective PXP Hamiltonian can be obtained by performing the Schrieffer-Wolff transformation using the VAGP [12].

D. Many-body dark states
The local structure of singularities of the AGP near the macroscopically degenerate points also allows for the existence of special states that are simultaneously eigenstates of the Hamiltonian and are annihilated by (or are possibly other eigenstates of) the leading divergent part of the AGP. From Eq. (7) it is clear that such states should be largely immune to any time-dependent protocols λ(t). They are thus approximately dark states.
Let us start by analyzing such states near the singularity at (2, 0). From Eq. (26) it follows that the divergent part of the AGP in any direction except the radial one scales as We can readily see that A s has many zero eigenstates that are simultaneous eigenstates of H at r = 0. An example of such a state is There are (exponentially) many other such dark states, which can e.g. be obtained by increasing the length of the domains of | ↑ spins. From Eqs. (7) and (8), the time evolution of such a | ψ 1 in the co-moving basis under an arbitrary time-dependent protocol is given by where A n is the remaining non-divergent part of the AGP as We see that the state |ψ 1 is unaffected by the term A s , the main source of diabatic excitations in general states, and is thus only weakly excited. Because this statement is general and is not tuned to the details of the protocol λ(t), this state approximately behaves as a many-body dark state. The remaining non-divergent terms A n entering Eq. (30) can be further suppressed by means of local . Even for the unassisted protocol, the final energy variance is much smaller for | ψ1 than for | ψ2 . Introducing a local counterdiabatic term rapidly decreases the final energy variance in the dark states, whereas the energy variance remains largely unchanged in the Néel state (see insets).
CDD. As we show in Sec. VI, A n has a well defined expansion in terms of local operators and thus the dark states only acquire local dressing near singularities and remain highly nonthermal (with e.g. low entanglement entropy) even far from the singularity, in the ergodic regime. To demonstrate the advantage of the many-body dark state in the context of quantum state preparation, we consider a CD protocol with the VAGP, starting at the singular point (2, 0) and subsequently increasing the transverse magnetic field up to the point (2, 0.5). We consider two scenarios, starting with two different initial states, a dark state | ψ 1 and a bright (non-dark) Néel state that is not annihilated by the singular part of the AGP: We choose the protocol given by Eq. (25) with protocol duration T = 1. The results of the simulations are shown in Fig. 6. In Appendix D, we analyze the dressing of less symmetric dark and bright initial states, and show that they exhibit a very similar qualitative behavior. Even for the unassisted protocol (blue lines), we can already see in the figure that the energy variance of the dressed dark state is a factor of 20 smaller than that of the bright state. This ratio quickly increases if we increase the protocol duration. The difference between the dark and bright states becomes even more pronounced in the presence of the local CD term. We see that the energy variance of the bright Néel state | ψ 2 is almost unaffected by the counterdiabatic term, only decreasing from 1.496 to 1.468 as we go from the unassisted protocol to the CDD with the 3-body ansatz. On the other hand, the energy variance of the prepared dark state reduces from 0.085 (unassisted) in the unassisted protocol to 0.001 (3-body CDD) for the dark state. Such a small energy variance implies that the prepared state is very close to an eigenstate of the system. The fact that this state is prepared in a short time T = 1 using a local CD Hamiltonian also implies that this state is nonthermal, e.g. it exhibits area law entanglement.
It is easy to check that the dark states, i.e. the zeroenergy eigenstates of the P Y P Hamiltonian, are simultaneously the zero-energy eigenstates of the low-energy effective P XP Hamiltonian. Interestingly, the AGP allows us to find these special states without prior knowledge of the effective Hamiltonian.
One can similarly analyze the structure of the AGP near the other singularity at (0, 0). From Eq. (26) it follows that the divergent part of the AGP is given by This operator clearly annihilates two pairs of states: i) fully-polarized states | ↑↑ . . . ↑↑ and | ↓↓ . . . ↓↓ and ii) the two Néel states | ↑↓↑↓ . . . ↑↓ and | ↓↑↓↑ . . . ↓↑ . The two Néel states are clearly the degenerate ground states, such that it is not surprising that they can be efficiently dressed locally as we introduce a nonzero finite magnetic field. The two ferromagnetic states are the most excited states, i.e. the states with maximal energy. As we turn on the Z-magnetic field, one of the polarized states remains the most excited state -it is again not surprising that this state can be locally dressed. However, the second polarized state quickly enters the energy continuum and yet, because it is annihilated by A s , it only weakly hybridizes with other states and remains highly non-thermal. This dark state was recently discovered in Ref. [40] (cf. Fig. 4 there) as a state with anomalously low entanglement. Interestingly, in this case the ground and most excited states can be immediately determined as the zero states of the AGP, without any need to diagonalize the full Hamiltonian.

IV. FLOW DIAGRAM WITH THE HIGHER ORDER VARIATIONAL ANSATZ
A. Scaling of the VAGP norm with the ansatz size Having analyzed the emerging adiabatic flow diagram within the 3-body variational ansatz, we now consider what happens when we increase the support of the VAGP to more than three sites. First, we study how the norm of the VAGP changes with the increasing ansatz size k. A slow increase of ||A λ || with k indicates that increasing the support of the ansatz only has a small effect on the VAGP, such that its local approximation is stable and accurate. Conversely, a fast increase of ||A λ || with ansatz size would indicate that the exact AGP is highly non-local and the local variational ansatz is not very stable. In Fig. 7, we analyze the norms of the VAGP in the optimal (blue) and orthogonal (red) directions at two different sets of couplings: (5/3, 1/10) (top) and (1/3, 1/3) (bottom). The first point is close to the g = 0 classical Ising line and relatively far from the singular points, whose structure is explained below. The second point is dominated by its proximity to the (0, 0) singularity, but it is not too close to it. In both cases we observe a large anisotropy between the optimal and orthogonal directions starting from k = 3. In particular, we see that the norm of the VAGP in the orthogonal direction rapidly increases to a large value as k reaches 3 and then remains relatively flat for the first set of couplings, and increases more gradually with k for the second set of couplings. In  , g) and the arrows denote the optimal direction for deformations (δh, δg).
The color now represents the logarithm of the ratio of the norm in the optimal direction over that in the orthogonal direction, ranging from blue (nearly anisotropic) to yellow (highly isotropic). Source flows are clearly visible not just at (h, g) = (0, 0) and (2, 0), but also at (1, 0) and (2/3, 0).
both cases the AGP norm in the optimal direction increases slowly with k. As we will show below, when we keep increasing the ansatz size, new singularities affecting the VAGP start to emerge. These singularities can discontinuously change the optimal direction, at the same time drastically reducing the anisotropy of the AGP.
Hence, the optimal direction is no longer radial, except exactly at the singularity, where θ = π/2. Since the operator part of the diverging contribution to the AGP contains four-body operators, this singularity will only manifest in the VAGP if we use a 4-body ansatz or higher. This is exactly what is shown in Fig. 1, where the flow diagram for the 5-body ansatz contains sources/sinks at both the 3-body singularities (0, 0) and (2, 0) and the additional singularity (1, 0). Increasing the support of the ansatz will lead to additional singularities, which can be captured in higherorder terms in the perturbative expansion. As such, higher-order singularities will become even more suppressed in orders of r, such that they will manifest themselves only some distance away from the degenerate g = 0 line. In Fig. 8, we show the flow diagram for the 8-body variational ansatz. The arrows again indicate the optimal directions, and the color now represents the anisotropy, i.e. the ratio of the VAGP norm along the optimal and the orthogonal directions, with yellow indicating a higher anisotropy. New singularities at h = 1 and h = 2/3 become visible in this plot, accompanied by additional, non-radial, structures around them.
Clearly, the leading-order singular term can be singled out either perturbatively or variationally. As argued above, the corresponding operators should connect states that are exactly degenerate at the corresponding singular point. In Table I we summarize these leadingorder operators and illustrate how they connect degenerate states through correlated spin flips, inducing both the macroscopic degeneracies in the eigenspectrum and the divergences in the VAGP.
We note an interesting feature following from Fig. 8: as we increase the size of the variational ansatz, in some regions the optimal direction can switch. This is most clearly visible near the point h = 1 and small g. Within the 3-body ansatz, the optimal direction is nearly horizontal (cf. Fig. 2), while in the higher-body ansatz (k > 4) the optimal direction is nearly vertical (cf. Fig. 8). This discontinuity indicates that it is impossible to improve the accuracy of the VAGP in the horizontal direction by increasing the support of the variational ansatz: the new singularity prevents us from doing so. The only way to continue improving local state preparation is to change the direction. It is clear that such a sudden change should introduce some ambiguity in finding the optimal path in the space of couplings in the vicinity of the singularity. Indeed, we see that regions of small anisotropy surround the singularity at (1, 0) -in such regions the difference between the optimal and the orthogonal directions is less pronounced.

V. VAGP AND APPROXIMATELY CONSERVED OPERATORS
As we discussed above, the VAGP for deformations along the direction λ j is found by minimizing the norm of the operator G j (cf. Eqs. (10) and (11)). If the VAGP is exact, then G j is a conserved operator conjugate to the direction λ j . However, for an approximate VAGP, G j is only approximately conserved because it has a nonzero commutator with the Hamiltonian. It is clear that the norm of the commutator [G j , H] is a measure for the accuracy of this approximate conservation law: the smaller the norm, the better the conservation law. In some sense, this norm serves as a proxy to the magnitude of the difference between the exact and the local variational AGP. If this difference is small, we can simultaneously implement accurate local counterdiabatic driving and construct a local nearly-conserved operator. These qualitative considerations are indeed correct, as we show below by analyzing the accuracy of such conservation laws in the optimal directions at different couplings and different ansatz sizes.
A more convenient and physical measure characterizing the accuracy of the conservation law is the lifetime of G j measured in an eigenstate |n of the Hamiltonian. The latter can be computed from the short-time expansion of the connected non-equal time correlation function [52]: . (34) From this expansion one can define a state-averaged normalized decay rate (inverse lifetime) for the operator G j as A small decay rate indicates that the operator G j is nearly conserved, at least up to times of the order 1/Γ j . For the exact AGP, obviously, Γ j = 0. In Fig. 9 we show the lifetimes of the operators G j computed in the optimal direction, i.e. the direction shown by arrows in Fig. 8, as a function of i) h at fixed g = 0.2 (panel a); ii) as a function of g at fixed h = 0.15, and iii) as a function of the total magnetic field along the diagonal direction h = g. Different lines on each panel refer to different ansatz sizes. In all the cases we chose the direction λ j to be the optimal one for the corresponding ansatz size. In the panel (a), showing Γ j as a function of h at a fixed small value of g, we see several characteristic features. First of all, it is clear that increasing the ansatz size increases the lifetime of the nearly-conserved operators. Furthermore, the decay rate exhibits non-monotonic peaks at the singular points of the AGP. As we increase the ansatz size Γ j becomes more sensitive to the higher-order singularities. Thus the effect of the singularity near h = 2 is very strong at k = 3, i.e. at the 3-body ansatz level but becomes very small for larger k. This picture is consistent with our previous analysis, suggesting that the divergent contributions to the VAGP, corresponding to leading singularities, are local and as such can be eliminated by the local VAGP. Higher-order singularities then require a VAGP with increasing support. Another very interesting feature emerges if we analyze the dependence of Γ j on g at fixed small h = 0. 15 [panel (b)]. Namely, the decay rate exhibits a clear maximum near g = 1, corresponding to the quantum critical point at zero temperature [68]. Interestingly, the maximum in Γ is clearly pronounced despite the fact that we analyze the operator lifetimes at infinite temperature, where static observables do not exhibit any signatures associated with criticality, consistent with recent results from Ref. [40]. At h = 0, i.e. in the limit of the integrable transverse field Ising model, this result is known from prior work [1,16]. The plot shown in Fig. 9 suggests that, even if the integrability is broken, the maximum of Γ j remains well-defined and again highlights how temperature plays a much less important role when we define quantum criticality through the diabatic response encoded in the AGP.

VI. PERTURBATIVE EXPANSION
In this final section we present a derivation of the divergences appearing in the VAGP by developing a perturbative expansion of the exact AGP in small g near g = 0, i.e. near the classical Ising limit, using the integral representation of the AGP given by Eq. (15). We will outline only a sketch of the derivation here and provide some key results, further details of all calculations can be found in Appendix F.
We will denote the Hamiltonian at the solvable point g = 0 as H 0 = ZZ+hZ and find a perturbative expansion for A λ at H = H 0 + gX in powers of g for general ∂ λ H, The first-order contribution can be found by setting g = 0 in Eq. (15), where any time-dependence is taken to be in the interaction picture, (∂ λ H (0) )(t) ≡ e iH0t (∂ λ H)e −iH0t . The next order can be found by taking the derivative w.r.t. g in Eq. (15), where In order to simplify the notations we use X(t) instead of X (0) (t). Higher-order terms can be found by taking higher-order derivatives of Eq. (15), leading to an iterative evaluation scheme. We will only analyze the first two orders here. We will separately calculate the dominant terms for ∂ λ H = X and ∂ λ H = Z, yielding A g and A h correspondingly. Given a general perturbation (δh, δg) ∝ (cos ϕ, sin ϕ), we can write A λ (ϕ) = cos ϕA h + sin ϕA g .
Given that H 0 = ZZ + hZ, in the interaction picture Z (0) (t) = Z is time-independent, and hence A (0) h = 0. For A g , we need to first evaluate X(t), which can be done analytically (see Eq. (F2) and Appendix F). It represents a sum of eight different independent operators with support up to k = 3 with time-dependent coefficients. For h = 0, 2 the integral of X(t) is well behaved in the limit → 0 and we can find This expression clearly diverges at h = 0 and h = 2. Collecting the diverging terms near these singularities, we recover the expressions quoted earlier (Eqs. (26) and (27)) in the limit ϕ → π/2 and θ → 0. Exactly at the singular points the divergent terms commutes with the Hamiltonian H 0 and can be subtracted from the AGP. This sudden discontinuity is not accidental, since the direction along g becomes exactly radial at the singular points, which is optimal. The cancellation of divergences also follows from Eq. (F2) and arises from the fact that the limits → 0 and h → 0, 2 do not commute. An explicit evaluation of Eq. (37) at h = 0 returns: similarly at h = 2 we find The first non-vanishing contribution to A h is A h , which can be immediately obtained from Eq. (38) (see again Appendix F for details): In a similar fashion, one can compute an exact analytic expression for A (1) g , showing the emergence of the new singularity at h = 1. This expression is rather long, so is is only explicitly given in the Appendix F.

Interestingly, while formally A
(1) h is obtained as a higher-order term than A (0) g , it contains the same type of singularities at h = 0 and h = 2. Moreover, it also only contains terms with support of up to three sites: both these terms will appear in, e.g. the 3-body variational ansatz. Physically, A (1) h plays the same role as A (0) g because both appear as the leading non-vanishing contributions to the AGP in the perturbative expansion. For this reason it suffices to analyze the following "leading order" perturbative AGP: As we will show next, the AGP in this form allows us to understand key features of the adiabatic flows near the singularities at (0, 0) and (2, 0). Using Eqs. (40) and (43), we can minimize the norm of the perturbative AGP (44) with respect to ϕ and find the optimal direction as The corresponding perturbative flow diagram is shown in Fig. 10. It is clearly highly similar to the variational flow diagram obtained for the 3-body variational ansatz (cf. Fig. 2), confirming how the (numerically straightforward) variational approach is able to identify the most important local contributions to the AGP. It is easy to check that from Eq. (44) we can recover the asymptotic behavior of the AGP close to the singularities at (0, 0) and (2, 0) (cf. Eqs. (26) and (27)).

VII. CONCLUSIONS
We developed a general approach for analyzing the adiabatic landscape in systems described by a family of Hamiltonians characterized by several controls (couplings). This approach is based on minimizing the norm of the local variational adiabatic gauge potential, which serves as the local generator of adiabatic transformations.
We applied this method to a one-dimensional Ising model in the presence of both transverse and longitudinal fields. In this model we determined the optimal directions as those where the norm of the adiabatic gauge potential is minimal, which can be used to immediately define continuous paths along which diabatic effects are suppressed (cf. Fig. 1). Along these optimal paths one can design highly efficient local and experimentally feasible counterdiabatic driving protocols. These paths are also useful for various other applications, including finding local nearly-conserved operators, dressed elementary excitations such as quasiparticles or domain walls (see also Refs. [12,40]), constructing effective Hamiltonians via the Schrieffer-Wolff transformation, numerically computing approximate eigenstates using efficient numerical methods such as the DMRG-X algorithm [72], designing optimal paths for quantum annealing protocols, suppressing dissipative losses in thermal machines and more. Interestingly, finding these optimal paths does not require diagonalizing the Hamiltonian of the system either exactly or approximately and can be done even in the thermodynamic limit.
We found that these optimal paths always start/terminate at the points corresponding to Hamiltonians exhibiting macroscopic degeneracies of the spectrum, which play a role similar to the role of quantum critical points in equilibrium phase diagrams. As we approach these singularities, the anisotropy between the optimal and the orthogonal directions diverges. The most divergent contributions to the adiabatic gauge potential are local and can be singled out either perturbatively or variationally. Increasing the support of the variational gauge potential, additional (weaker) divergences start to emerge, strongly affecting the flow diagram in their vicinity. Close to these singularities we can identify special dark states: mutual eigenstates of the Hamiltonian at the singular point and the divergent part of the adiabatic gauge potential. These dark states are highly robust against various time-dependent perturbations and can be efficiently locally dressed by the non-divergent part of the VAGP. They persist deep in the ergodic regime extending far away from the singularities. Physically, these dressed dark states correspond to spin configurations that can remain non-thermal for extremely long times. Our method provides a general prescription of finding such non-thermal states in interacting systems.
Finally, we showed that the optimal directions are associated with the existence of local nearly-conserved operators. Thus there is an interesting and direct connection between our ability to perform efficient local adiabatic transformations along particular directions and the existence of long-lived operators, which are locally dressed deformations of the Hamiltonian along these optimal directions.
FIG. 11. Scaling with r of each term in the VAGP along the orthogonal (a) and optimal (b) direction. All the norms of the non-vanishing terms in the 3-body ansatz are shown. The horizontal axis is r = h 2 + g 2 and the vertical one is the norm of each term |cn| 2 (markers). In the orthogonal direction, the dominant term O(1/r) is given by Y − ZY Z, whereas the dominant term along the optimal direction is O(1) and given by Y Z + ZY . All other terms exhibit higher-order scaling, where selected fits (full lines) are included as a guide to the eye. In the Appendix, we present a short argument for why the radial direction is generally the optimal one near singular points. Let us assume that the Hamiltonian of the system can be written as where H 0 is a Hamiltonian describing a macroscopically-degenerate point and V is some generic perturbation. The exact AGP can be represented in the degenerate eigenbasis of the instantaneous Hamiltonian as [1] A λ = i where H | n = n | n . It is clear from this expression that the AGP generally diverges in the limit → 0 as 1/ε, since in first-order perturbation theory n − m ≈ n|V |n − m|V |m . Note that in the proper eigenstates the matrix V is approximately diagonal within each degenerate manifold. This divergence is, however, canceled if the deformation ∂ λ H is diagonal (to the first order) in ε in this basis. In particular, this is the case when λ = ε, which precisely corresponds to radial deformations defining the optimal directions (cf. Sec. III C).
be readily done as Its commutators with X(t) and Z follow as Note that the even terms in Y do not contribute to the AGP: as follows from e.g. Eq. (E2) the AGP is explicitly imaginary for a real Hamiltonian (see also Ref. [37]). These terms, however, will contribute to higher order corrections to the AGP. Using these expressions together with Eqs. (37) we can immediately recover the leading-order contributions to the AGP shown in Sec. VI. Likewise, using (38) we can find the first subleading corrections, which we will show below for different values of h: • h = 0, 1, 2 From the expansion (F6) we recover both the singularities close h = 0 and h = 2, as discussed in the main text, as well as the emergence of a new singularity close to h = 1, which one can check is proportional to P (XY + Y X)P . One can further check that a similarly divergent term, also proportional to P (XY + Y X)P , appears in the second order correction to A h , which we only show for completeness away from singularities, i.e. h = 0, 1, 2: Collecting the terms that will be singular at h = 1 appearing in the expressions for A (1) g and A (2) h , we can obtain Eq. (33).