## Provably Efficient Adiabatic Learning for Quantum-Classical Dynamics

Changnan Peng, 1 Jin-Peng Liu, 2, 3 Gia-Wei Chern, 4 and Di Luo 2, 5, 6, ∗

1 Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 2 Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 3 Simons Institute and Department of Mathematics, University of California, Berkeley, CA 94704, USA 4 Department of Physics, University of Virginia, Charlottesville, VA 22904, USA 5 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions 6 Cambridge, MA 02138, USA

Department of Physics, Harvard University, (Dated: August 9, 2024)

Quantum-classical hybrid dynamics is crucial for accurately simulating complex systems where both quantum and classical behaviors need to be considered. However, coupling between classical and quantum degrees of freedom and the exponential growth of the Hilbert space present significant challenges. Current machine learning approaches for predicting such dynamics, while promising, remain unknown in their error bounds, sample complexity, and generalizability. In this work, we establish a generic theoretical framework for analyzing quantum-classical adiabatic dynamics with learning algorithms. Based on quantum information theory, we develop a provably efficient adiabatic learning (PEAL) algorithm with logarithmic system size sampling complexity and favorable time scaling properties. We benchmark PEAL on the Holstein model, and demonstrate its accuracy in predicting single-path dynamics and ensemble dynamics observables as well as transfer learning over a family of Hamiltonians. Our framework and algorithm open up new avenues for reliable and efficient learning of quantum-classical dynamics.

Introduction. - Efficient simulation of quantumclassical hybrid dynamics is crucial to multi-scale modelings of a wide range of physical systems, opening new avenues for advancements in material science, chemistry, and drug discovery by providing a more comprehensive understanding of molecular interactions [1-7]. A common hybrid dynamics approach relies on the adiabatic approximation, where two well-separated timescales of a system allow one to treat the slow dynamics classically while quantum calculations are used to solve the fast, often electronic, degrees of freedom that adiabatically follow the classical dynamics. A well-known example is the Born-Oppenheimer approximation which is widely used in ab initio molecular dynamics [8]. The significance of quantum-classical hybrid dynamics lies in its potential to revolutionize how we model and predict the behavior of complex systems, ultimately pushing the frontiers of both fundamental research and practical applications.

fields facilitate the simulation of large intricate systems by accurately capturing the essential quantum mechanical effects while maintaining computational efficiency. Despite intensive studies and wide applications of ML force field models over past decades, critical questions on the error bound, sample complexity and generalizability of the ML methods have remained unresolved.

The simulations of quantum-classical dynamics, however, is computationally challenging due to not only an exponentially large Hilbert space of quantum subsystems and repeated time-consuming quantum calculations at every time step, but also the nonlinear differential equation coupled both the quantum and the classical variables. In the past two decades, machine learning (ML) has emerged as a powerful tool in developing force fields and inter-atomic potentials for ab initio molecular dynamics [9-20]. ML force-field approaches have recently been generalized to enable large-scale dynamical simulations of condensed-matter lattice systems [21-25]. This approach leverages massive datasets of quantum mechanical results to train models that can predict the potential energy surfaces with high precision. ML-enhanced force

In this work, we establish a generic theoretical framework for analyzing quantum-classical adiabatic dynamics with learning algorithms. We start with the analysis of the approximately constant linear model, derive the error bounded condition for the non-linear model, and introduce the relaxation method to check the error bounded property for a generic model, which lays down a solid foundation for the reliability of learning algorithms in quantum-classical adiabatic dynamics. Inspired by the recent development of quantum information theoretic learning theory for quantum many-body systems [26-39], we develop a provably efficient adiabatic learning (PEAL) algorithm for quantum-classical dynamics, which offers a sample complexity scaling logarithmically with the system size and favorable scaling of evolution time. We benchmark our algorithm on the Holstein model and demonstrate accurate prediction of the single path dynamics and ensemble dynamics observables, as well as transfer learning across different couplings between quantum and classical degrees of freedom.

Adiabatic Quantum-Classical Dynamics Learning. -We consider a general quantum-classical Hamiltonian:

$$\int _ { \mathrm e } ^ { \vdash } \hat { H } = \hat { H } _ { q } + \sum _ { \alpha, i } g _ { \alpha } \hat { O } _ { \alpha, i } G _ { \alpha, i } ( \vec { P }, \vec { Q } ) + H _ { c l } ( \vec { P }, \vec { Q } ), \quad ( 1 )$$

FIG. 1: Schematic diagram for PEAL.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-f78ddced-2e72-41eb-b477-92ea350b0c87.jpg)

where ˆ H q is the Hamiltonian operator for the quantum degrees of freedom, g α is the quantum-classical coupling coefficient, α is the index of different types of couplings, i ∈ { 1 , . . . , L } is the label of local regions (e.g. lattice sites), ˆ O α,i is the quantum operator that enters the coupling, G α,i ( ⃗ , Q P ⃗ ) is a general function of the classical degree of freedom ⃗ Q = ( Q , . . . , Q 1 L ) and its canonical momentum ⃗ P = ( P , . . . , P 1 L ), and H cl ( ⃗ , Q P ⃗ ) is the classical Hamiltonian.

We study the adiabatic evolution dynamics driven by the quantum-classical Hamiltonian in Eq. (1). In this adiabatic limit, similar to the Born-Oppenheimer approximation in quantum molecular dynamics, the quantum subsystem is assumed to quickly relax to the ground state of the total Hamiltonian ˆ , H while the classical degrees of freedom follow the Hamilton's equations of motion (EOM):

$$\frac { d } { d t } Q _ { j } = \sum _ { \alpha, i } g _ { \alpha } \langle \hat { O } _ { \alpha, i } \rangle \frac { \partial } { \partial P _ { j } } G _ { \alpha, i } + \frac { \partial } { \partial P _ { j } } H _ { c l }, \quad \ \ ( 2 ) \, \sum _ { \alpha \, \text{act si} } ^ { \, \text{mom} } \, \text{act si} \, \text{ } \, \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{                                                                                                                                                                                                         }                                                                                                                                                                                                       
    \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{
    \text{                                                                                                                                                                                                     }   \text{ } \text{                                                                                                                                                                                                 }       }                                                                                                                                                                                                 }      
    \text{                                                                                                                                                                                                 }    }                                                                                                                                                                                                    }    }   }                                                                                                                                                                                                     }   }                                                                                                                                                                                                 }    }  
    \text{                                                                                                                                                                                               }         }                                                                                                                                                                                               }        
    \text{                                                                                                                                                                                           }             }                                                                                                                                                                                           }            
    \text{                                                                                                                                                                                         }               }                                                                                                                                                                                         }              
    \text{                                                                                                                                                                                     }$$

$$\frac { d } { d t } P _ { j } = - \sum _ { \alpha, i } g _ { \alpha } \langle \hat { O } _ { \alpha, i } \rangle \frac { \partial } { \partial Q _ { j } } G _ { \alpha, i } - \frac { \partial } { \partial Q _ { j } } H _ { c l } - \gamma P _ { j }, \ ( 3 ) \quad \text{satisf}.$$

where j ∈ { 1 , . . . , L } , ⟨·⟩ is the ground state expectation, and γ &gt; 0 is the damping coefficient due to dissipation.

The conventional way to solve the above equations is through an iterative scheme, where one first updates { Q j ( ) t , P j ( ) t } based on the classical EOM, and then solve the ground state from ˆ . H The updated ground state is used to compute expectation values ⟨ ˆ O α,i ⟩ which determine the driving terms of the EOM for the next step. However, repeated ground-state calculations of ˆ H at every time step could be computationally costly using quantum state solvers (QSS), such as exact diagonalization (ED) [40], density functional theory [41], DMRG [42], neural network [43-47], and quantum computers [48-51] .

We consider a general quantum-classical ML model for learning such dynamics, which we call Adiabatic Dynamics Model Learning (ADML). ADML consists of two components, which uses machine learning to predict the quantum observables and evolves the classical observables using classical numerical schemes. The ML force- field models widely used in quantum molecular dynamics can be viewed as special classes of ADML [9-20]. Our goal here is to predict the adiabatic dynamics of Eq. (1) with a learning-based approach. Given access to a dataset ⋃ N s s =1 { ( g α , P ⃗ ( ) t , Q t , ⃗ ( ) ⟨ ˆ O ⟩ ( )) t s : t ∈ T } s , where T s is the set of sampled time steps, N s is the number of initial conditions { ( g α , P ⃗ ( t init ) , Q t ⃗ ( init )) s } N s s =1 sampled from a distribution D init , and ˆ stands for O ˆ O α,i in Eq. (1) or other operators of insterest but not in the Hamiltonian, the task is to design ADML for predicting the dynamics starting with other { ( g α , P ⃗ ( t init ) , Q t ⃗ ( init )) } from D init . In the following, we analyze and derive the error bounded conditions for ADML.

̸

(i) Approximately constant linear model. To serve as a starting point, we consider a simple example of Eq. (1), ˆ H = ˆ H q + ∑ i gO Q ˆ i i + ∑ ( i 1 2 M P 2 i + 1 2 kQ 2 i ) , where we only consider one type of quantum-classical coupling with G α,i ( ⃗ , Q P ⃗ ) = Q i , and the classical Hamiltonian H cl ( ⃗ , Q P ⃗ ) is for simple harmonic oscillators with mass M and spring constant k . Further assuming that during the dynamical process we are interested in, the response ∂ O ⟨ ˆ i ⟩ ∂Q i is approximately a constant, and the offdiagonal response ∂ O ⟨ ˆ i ⟩ ∂Q j ( i = j ) is approximately zero, the EOM of the system is then reduced to that of independent simple harmonic oscillators.

Focusing on the classical degree of freedom Q i , we can view the quantum-classical coupling gO Q ˆ i i as a driving force on the oscillator. Since the EOM in Eq. (2) and (3) are approximately linear in this case, the accumulated momentum and position errors between ML and the exact simulation using QSS, p t ( ) = P i, ML ( ) t -P i, Exact ( ), t q t ( ) = Q i, ML ( ) t -Q i, Exact ( ) t ( i index suppressed), also satisfy a similar EOM:

$$\frac { d } { d t } q ( t ) & = \frac { 1 } { M } p ( t ), & ( 4 ) \\ \frac { d } { d t } p ( t ) & = F ( t ) - K q ( t ) - \gamma p ( t ) + o ( q ( t ) ), & ( 5 )$$

$$\frac { a } { d t } p ( t ) = F ( t ) - K q ( t ) - \gamma p ( t )$$

where we define the error force F t ( ) = -gδ O ⟨ ˆ i ⟩ ( ) t and the error stiffness K = g ∂ O ⟨ ˆ i ⟩ ∂Q i + k . δ O ⟨ ˆ i ⟩ ( ) t = ⟨ ˆ O i ⟩ ML ( ⃗ Q ML ( )) t - ⟨ ˆ O i ⟩ Exact ( ⃗ Q ML ( )) t is the ML prediction error at the t -th time step. o q t ( ( )) means higher order terms in q t ( ) (See Supplemental Material for the derivation). Given a bounded F t ( ), even if it is tuned to drive the oscillator optimally, as long as K,γ &gt; 0, the oscillator cannot be driven to infinite amplitude, i.e. the accumulated momentum and position errors are bounded. For F t ( ) to be bounded, it suffices to have δ O ⟨ ˆ i ⟩ ( ) bounded, which motivates us to define an Error t Bounded Property such that when δ O ⟨ ˆ i ⟩ ( ) is t bounded, p t ( ) and q t ( ) are also bounded. More precisely, we define the following

Definition 1 (Error Bounded Property) . A model satis-

fies the Error Bounded Property with respect to δ O ⟨ ˆ i ⟩ ( ) t for t ∈ [ t init , t end ] , if and only if the following claim is true: ∃ C , C q p &gt; 0 such that ∀ ϵ &gt; 0 , if ∀ t ∈ [ t init , t end ] , | δ O ⟨ ˆ i ⟩ ( ) t | 2 ≤ ϵ , then there are | q t ( ) | ≤ C q √ ϵ and | p t ( ) | ≤ C p √ ϵ , ∀ t ∈ [ t init , t end ] .

Proposition 1. The approximately constant linear model satisfies the Error Bounded Property if K &gt; 0 .

The proof of Prop. 1 is in Supplemental Material.

(ii) Non-linear model. We can generalize the approximately constant linear model to allow non-linearity. We drop the assumption that the response ∂ O ⟨ ˆ i ⟩ ∂Q i is approximately a constant, allow G α,i ( ⃗ , Q P ⃗ ) to be a non-linear function of Q i , and allow the potential 1 2 kQ 2 i to include non-quadratic component in Q i , as long as we can Taylor expand the potential at its minimum. These generalizations can be absorbed by a redefined F t ( ) and a timedependent error stiffness K t ( ). Unlike the approximately constant linear model, the oscillator could have infinite amplitude even if K t ( ) &gt; 0 for all the time. We present a condition in Supplemental Material which guarantees a bounded bounded in the worst case scenario. Summarized as an informal theorem, we have:

Theorem 1 (Error Bounded Condition for non-linear model (Informal)) . If K t ( ) &gt; M γ/ ( 2) 2 + C with a positive constant C for all t , and the error stiffness K t ( ) 's fluctuation, as well as G α,i ( ⃗ , Q P ⃗ ) and its first derivatives are bounded, then the non-linear model satisfies the Error Bounded Property in Def. 1.

The proof for Thm. 1 is provided in Supplemental Material. We can further apply it to the Hamiltonian in normal mode with the quadratic momentum under Fourier transform, which could appear in a more general setup.

(iii) General Relaxation Method. We now come back to the most general ADML. We allow arbitrary G α,i ( ⃗ , Q P ⃗ ) and H cl ( ⃗ , Q P ⃗ ) in Eq. (1), and we make no assumption on the response ∂ O ⟨ ˆ i ⟩ ∂Q j . Because in the most general case the classical degrees of freedom are no longer decoupled, we restore the i index of the accumulated momentum and position errors, p i ( ) = t P i, ML ( ) t -P i, Exact ( ), t q i ( ) = t Q i, ML ( ) t -Q i, Exact ( ). t

For a general ADML, the EOM for the errors q i ( ), t p i ( ) can in general be derived from Eq. (2) and (3): t

$$\frac { d } { d t } q _ { i } = \sum _ { j } \left ( \mathcal { K } _ { q _ { j }, q _ { i } } q _ { j } + \mathcal { K } _ { p _ { j }, q _ { i } } p _ { j } \right ) + \mathcal { F } _ { q _ { i } } + o ( q, p ), \quad ( 6 ) \quad \text{symm} \\ \text{ple cco} \\ \text{follow}$$

$$\frac { d } { d t } p _ { i } = \sum _ { j } \left ( \mathcal { K } _ { q _ { j }, p _ { i } } q _ { j } + \mathcal { K } _ { p _ { j }, p _ { i } } p _ { j } \right ) + \mathcal { F } _ { p _ { i } } - \gamma p _ { i } + o ( q, p ), \quad \text{The} \epsilon$$

(7)

where the error stiffness matrix K q p ( ) j ,q ( p ) i and the error force vector F q p ( ) i depend on ( ⃗ , Q P ⃗ ), and F q p ( ) i is linear in δ O ⟨ ˆ α,i ⟩ ( ) (see Supplemental Material for details). t

While it is difficult to write down an error bounded condition for a general ADML, we propose a relaxation method to provide a sufficient (but not necessary) condition to check the Error Bounded Property in Def. 1. The idea is to consider the worst case scenario. If in the worst case scenario the error still converges, then it is safe to use ADML to accelerate our dynamical simulation. The relaxation method is as follows. First, we require K elements, as well as G α,i ( ⃗ , Q P ⃗ ) and its first derivatives, are bounded. Second, assuming | δ O ⟨ ˆ α,i ⟩ ( ) t | 2 ≤ ϵ , we identify possible upper and lower bounds for the elements in K and F . The bounds just only need to be effective during the time range of the simulation. Third, to achieve the worst case scenario, we want to maximize dq /dt i when q i is positive, and minimize when negative (the same for p i ). Therefore, we insert the upper bound of K q p ( ) j ,q ( p ) i when q p ( ) i and q p ( ) j have the same sign, and the lower bound if the opposite sign. F q p ( ) i are adjusted to their upper or lower bounds accordingly. Fourth, we perform a classical simulation of the EOM in Eq. (6) and (7), with the worst case scenario stated above. Finally, if the simulation shows there exist constants C q i and C p i such that for any ϵ &gt; 0, there are | q i ( ) t | ≤ C q i √ ϵ and | p i ( ) t | ≤ C p i √ ϵ during the time range of interest, then the Error Bounded Property is verified with the relaxation method.

Provably Efficient Adiabatic Learning. - Next, we present a provably efficient learning algorithm for the above ADML models based on quantum information theory, which we call Provable Efficient Adiabatic Learning (PEAL). The PEAL algorithm, equipped with a learning model M and a classical ordinary differential equation (ODE) solver, consists of the following steps.

↦

(i) Data collection for training. We sample a set of N s initial conditions { ( g α , P ⃗ ( t init ) , Q t ⃗ ( init )) s } N s s =1 from a distribution D init . We evolve the system with QSS and ODE solver to get ( ⃗ P t , ⃗ ( ) Q t , ( ) ⟨ ˆ O ⟩ ( )) t s . For each s , we sample a set of time steps T s uniformly from [ t init , t end ]. (ii) Model training. We use the dataset ⋃ N s s =1 { ( g α , P ⃗ ( ) t , Q t , ⃗ ( ) ⟨ ˆ O ⟩ ( )) t s : t ∈ T } s to train a model M : ( g α , P , Q ⃗ ⃗ ) → ⟨ ˆ O ⟩ with the learning algorithm developed in [27]. (iii) Prediction. For any unseen new initial condition ( g α , P ⃗ ( t init ) , Q t ⃗ ( init )) new ∼ D init , PEAL outputs the dynamical trajectory ( ⃗ P t , ⃗ ( ) Q t , ( ) ⟨ ˆ O ⟩ ( )) t new for t ∈ [ t init , t end ], by alternatively updating ⟨ ˆ O ⟩ ( ) with t M and ⃗ P t , ⃗ ( ) Q t ( ) with ODE solver, integrated with our symmetry-preserving techniques shown later. The sample complexity and error bounds are summarized in the following informal theorem.

Theorem 2 (Provably Efficient Adiabatic Learning Theorem (Informal)) . When the Error Bounded Property in Def. 1 is satisfied, for T time steps quantum-classical adiabatic dynamics of an n -qubit gapped system, with sample complexity O (log( n )) , PEAL gives rise to controllable accumulated errors of classical variables and all k-local, bounded quantum observables scaling as (i) O ( √ T ) for

FIG. 2: PEAL (dashed) vs. exact simulation (solid) in single-path prediction. Blue curves correspond to standard learning with g = 1 4, and orange curves correspond to transfer learning with . g = 1 39. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-b026d74b-d997-44d6-bca5-6453cde363d3.jpg)

generic model M (ii) O ( √ log T ) for sub-Gaussian M 's prediction error (iii) independent on T for bounded M 's prediction error.

We leave the proofs for the above theorem in Supplemental Material. We note that the computational time for PEAL's prediction under a fixed g α is O (min { nN ,c n s ( ) } T ), where c n ( ) is the ODE solver per time step complexity and the nN s factor comes from the model M [27]. For Theorem 2, it can also be applied to an n -qubit system with ground states of exponentialdecay correlation functions using recent results [29]. We highlight that PEAL works for unseen couplings g α, new / ∈ { g α,s } N s s =1 , demonstrating transfer learning over a family of Hamiltonian. Furthermore, PEAL can predict any k -local observable ˆ O even if it does not appear in the Hamiltonian.

density wave (CDW) in zero temperature. In Supplemental Material, we derive the relation between phonon amplitude and CDW response function, showing the error stiffness K is always positive when it's close to a checkerboard configuration. We also numerically measure the error stiffness K t ( ) during the time range of interest, confirming the Error Bounded Property (see Supplemental Material). Therefore, we can apply PEAL to the 1D adiabatic Holstein model with a bounded error guarantee during the dynamics.

Numerical Experiments. - We consider the Holstein model [52-54] which describes the electron phonon interaction as follows:

$$H = - \sum _ { i, j } c _ { i } ^ { \dagger } c _ { j } - g \sum _ { i } \left ( c _ { i } ^ { \dagger } c _ { i } - \frac { 1 } { 2 } \right ) Q _ { i } + \sum _ { i } \left ( \frac { P _ { i } ^ { 2 } } { 2 M } + \frac { k Q _ { i } ^ { 2 } } { 2 } \right ), \quad \quad 6 \times \\ \intertext { m h o m e } \times \in \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ m h o m e \ sigma } \times 4 \times$$

where c i is related to the fermionic degree of freedom, while Q i and P i are the position and momentum of the phonon as classical degrees of freedom.

In the numerical experiment, we study such model on a 1D periodic chain with L = 50 sites. We set M = k = 1 and the damping coefficient γ = 0 1. . The electronic degree of freedom at time t will be the ground state of H e ( { Q t i ( ) } ) due to the fast relaxation, showing a 1D strong Anderson localization [55]. Therefore, PEAL's requirement of a ground state with exponential-decay correlation functions is satisfied. Analytically, the 1D adiabatic Holstein model always cools down to a charge

Since Eq. (8) has the U (1) and translation symmetry, we develop a symmetry-preserving PEAL. The U (1) global symmetry is respected by conserving the total electron density, and the translation symmetry by applying the same model on all sites. The U (1) symmetry preserving is a new feature in our PEAL that does not exist in previous literature and we have shown that the

FIG. 3: Test errors for standard learning (blue) and transfer learning (orange). Inner pannels: n PEAL vs. n Exact for standard learning with g ∈ G SL (Left) and transfer learning with g ∈ G TL (Right).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-b0991c50-2237-49cc-a878-4551f181ca85.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-65d952f9-9983-4aa3-9ab5-74a3cc70add0.jpg)

-

-

FIG. 4: Consistent agreement between PEAL and the exact simulation for the ensemble Q t Q i ( ) j ( ) correlation t at t = 100 of g = 1 4 (Left) and . g = 1 7 (Right). .

symmetry-preserving PEAL maintains provably efficient error bound in Supplemental Material.

We begin with the single-path prediction task by training a model using some initial conditions and certain g values. The goal is to predict observable dynamics from different initial conditions under both the training g values (standard learning) and unseen g values (transfer learning). Here, we choose g ∈ G SL = { 1 3 . , 1 32 1 34 . , . , 1 36 . , 1 38 1 4 . , . } for training. In Figure 2, we demonstrate the single-path prediction by PEAL. The solid curves are the exact simulation using exact diagonalization (ED) for QSS and RK4 for classical ODE solver [56], and the dashed curves are with PEAL. Blue curves are for g = 1 4 (standard learning), and orange . curves for g = 1 39 . (transfer learning). We present the time evolution of six different observables during the dynamics, which are the total charge density wave (CDW) ∑ i ( -1) i n i , the electron density at the first site n , the phonon amplitude at the first site Q , the phonon momentum at the first site P , the hopping term ⟨ c † 0 c 1 ⟩ , and the next-nearest-neighbor correlation ⟨ c † 0 c 2 ⟩ . We note that CDW is a sum of local observables, ⟨ c † 0 c 1 ⟩ does not enter classical EOM, and ⟨ c † 0 c 2 ⟩ does not even appear in the Hamiltonian. Nevertheless, all observables in the PEAL curves agree well with their corresponding exact simulations, demonstrating that PEAL guarantees a controllable error for k -local observable and well performs in transfer learning.

In Figure 3, we present the sample complexity of PEAL. We use the same training data in the singlepath prediction, build the test set with g ∈ G SL but with different initial conditions, and demonstrate transfer learning with g ∈ G TL = { 1 31 1 33 1 35 . , . , . , 1 37 . , 1 39 . } . For illustration, we consider the root mean square of | n PEAL -n Exact | as the test error. The test errors for both standard learning and transfer learning show a logarithmic scaling law with the number of samples. The inner panels also demonstrate that n PEAL in standard learning and transfer learning match well with n Exact .

Generically, the ensemble dynamics (i.e. the dynamics average over an ensemble of paths from different initial conditions) can provide more information of the underlying physics. In Figure 4, we compare the ensemble averages of the correlation function Q t Q i ( ) j ( ) t at time t = 100 for g = 1 4 and . g = 1 7. Here, . g = 1 4 prediction . is based on PEAL training in g ∈ G SL , and g = 1 7 pre-. diction is based on another model from PEAL training in g ∈ { 1 6 . , 1 62 . , 1 64 1 66 1 68 . , . , . , 1 7 . } where long-lasting domain walls exist. For g = 1 4, the system has cooled . down to a pure checkerboard configuration, and shows an oscillating correlation function. For g = 1 7, since the . system has domain walls, it shows a decaying correlation function and a large variance. PEAL provides a good agreement with the exact simulations for both the mean and the variance over different g values.

Conclusion. - In this work, we establish a comprehensive theoretical framework for analyzing quantumclassical adiabatic dynamics using learning algorithms. We provide a systematic analysis for the error bounded properties of the approximately constant linear model, the non-linear model, and generic models, ensuring the reliability of ADML. We develop a provably efficient adiabatic learning algorithm PEAL, demonstrating logarithmic scaling of sample complexity with system size and favorable evolution time scaling. Benchmarking PEAL on the Holstein model, we achieve accurate predictions of single-path dynamics and ensemble dynamics observables, with effective transfer learning across various quantum-classical coupling strengths. Our framework and algorithm opens up new directions for efficient learning in quantum-classical dynamics, including applying PEAL to quantum experiments, integrating advanced machine learning techniques, and extending to generic non-equilbrium processes.

Acknowledgements - The authors acknowledge helpful discussion with Max Metlitski, Lingyu Yang and Zhi Ren. DL acknowledges support from the NSF AI Institute for Artificial Intelligence and Fundamental Interactions (IAIFI). JPL acknowledges support by the National Science Foundation (PHY-1818914, CCF-1729369), the NSF Quantum Leap Challenge Institute (QLCI) program (OMA-2016245, OMA-2120757), and a Simons Foundation award (No. 825053). GWC acknowledges the support of the US Department of Energy Basic Energy Sciences under Contract No. DE-SC0020330. The authors acknowledge the MIT SuperCloud for providing HPC resources that have contributed to the research results reported within this paper.

∗ correspondence: diluo@mit.edu

[1] R. Kapral and G. Ciccotti, Mixed quantum-classical dynamics, The Journal of chemical physics 110 , 8919 (1999).

[2] R. Kapral, Progress in the theory of mixed quantumclassical dynamics, Annu. Rev. Phys. Chem. 57 , 129 (2006).

[3] H. Lin and D. G. Truhlar, Qm/mm: what have we

- learned, where are we, and where do we go from here?, Theoretical Chemistry Accounts 117 , 185 (2007).
- [4] M. W. van der Kamp and A. J. Mulholland, Combined quantum mechanics/molecular mechanics (qm/mm) methods in computational enzymology, Biochemistry 52 , 2708 (2013).
- [5] B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer, Hybrid quantum-classical approach to correlated materials, Phys. Rev. X 6 , 031045 (2016).
- [6] M. C. Melo, R. C. Bernardi, T. Rudack, M. Scheurer, C. Riplinger, J. C. Phillips, J. D. Maia, G. B. Rocha, J. V. Ribeiro, J. E. Stone, et al. , Namd goes quantum: an integrative suite for hybrid simulations, Nature methods 15 , 351 (2018).
- [7] O. T. Unke, M. St¨ ohr, S. Ganscha, T. Unterthiner, H. Maennel, S. Kashubin, D. Ahlin, M. Gastegger, L. Medrano Sandonas, J. T. Berryman, et al. , Biomolecular dynamics with machine-learned quantum-mechanical force fields trained on diverse chemical fragments, Science Advances 10 , eadn4397 (2024).
- [8] D. Marx and J. Hutter, Ab initio molecular dynamics: basic theory and advanced methods (Cambridge University Press, 2009).
- [9] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Phys. Rev. Lett. 98 , 146401 (2007).
- [10] A. P. Bart´k, M. C. Payne, R. Kondor, and G. Cs´nyi, o a Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett. 104 , 136403 (2010).
- [11] Z. Li, J. R. Kermode, and A. De Vita, Molecular dynamics with on-the-fly machine learning of quantummechanical forces, Phys. Rev. Lett. 114 , 096405 (2015).
- [12] A. V. Shapeev, Moment tensor potentials: A class of systematically improvable interatomic potentials, Multiscale Modeling &amp; Simulation 14 , 1153 (2016).
- [13] J. Behler, Perspective: Machine learning potentials for atomistic simulations, The Journal of Chemical Physics 145 , 170901 (2016).
- [14] V. Botu, R. Batra, J. Chapman, and R. Ramprasad, Machine learning force fields: Construction, validation, and outlook, The Journal of Physical Chemistry C 121 , 511 (2017).
- [15] L. Zhang, J. Han, H. Wang, R. Car, and W. E, Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics, Phys. Rev. Lett. 120 , 143001 (2018).
- [16] R. T. McGibbon, A. G. Taube, A. G. Donchev, K. Siva, F. Hernandez, C. Hargus, K.-H. Law, J. L. Klepeis, and D. E. Shaw, Improving the accuracy of m¨ller-plesset pero turbation theory with neural networks, The Journal of Chemical Physics 147 , 161725 (2017).
- [17] H. Suwa, J. S. Smith, N. Lubbers, C. D. Batista, G.-W. Chern, and K. Barros, Machine learning for molecular dynamics with strongly correlated electrons, Phys. Rev. B 99 , 161107 (2019).
- [18] S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Sch¨tt, and K.-R. M¨ller, Machine learning of acu u curate energy-conserving molecular force fields, Science Advances 3 , e1603015 (2017).
- [19] S. Chmiela, H. E. Sauceda, K.-R. M¨ uller, and A. Tkatchenko, Towards exact molecular dynamics simulations with machine-learned force fields, Nature Communications 9 , 3887 (2018).
- [20] H. E. Sauceda, M. Gastegger, S. Chmiela, K.-R. M¨ller, u and A. Tkatchenko, Molecular force fields with gradientdomain machine learning (GDML): Comparison and synergies with classical force fields, The Journal of Chemical Physics 153 , 124109 (2020).
- [21] P. Zhang and G.-W. Chern, Arrested phase separation in double-exchange models: Large-scale simulation enabled by machine learning, Phys. Rev. Lett. 127 , 146401 (2021).
- [22] S. Zhang, P. Zhang, and G.-W. Chern, Anomalous phase separation in a correlated electron system: Machinelearning enabled large-scale kinetic monte carlo simulations, Proceedings of the National Academy of Sciences 119 , e2119957119 (2022).
- [23] P. Zhang and G.-W. Chern, Machine learning nonequilibrium electron forces for spin dynamics of itinerant magnets, npj Computational Materials 9 , 32 (2023).
- [24] C. Cheng, S. Zhang, and G.-W. Chern, Machine learning for phase ordering dynamics of charge density waves, Phys. Rev. B 108 , 014301 (2023).
- [25] X. Cheng, S. Zhang, P. C. H. Nguyen, S. Azarfar, G.-W. Chern, and S. S. Baek, Convolutional neural networks for large-scale dynamical modeling of itinerant magnets, Phys. Rev. Research 5 , 033188 (2023).
- [26] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, and J. Preskill, Provably efficient machine learning for quantum many-body problems, Science 377 , eabk3333 (2022), https://www.science.org/doi/pdf/10.1126/science.abk3333.
- [27] L. Lewis, H.-Y. Huang, V. T. Tran, S. Lehner, R. Kueng, and J. Preskill, Improved machine learning algorithm for predicting ground state properties, Nature Communications 15 , 895 (2024).
- [28] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16 , 1050 (2020).
- [29] E. Onorati, C. Rouz´, D. S. Fran¸a, and J. D. Watson, Efe c ficient learning of ground &amp; thermal states within phases of matter, arXiv preprint arXiv:2301.12946 (2023).
- [30] J. Haah, R. Kothari, and E. Tang, Learning quantum hamiltonians from high-temperature gibbs states and real-time evolutions, Nature Physics , 1 (2024).
- [31] H.-Y. Huang, S. Chen, and J. Preskill, Learning to predict arbitrary quantum processes, PRX Quantum 4 , 040337 (2023).
- [32] M. C. Caro, H.-Y. Huang, M. Cerezo, K. Sharma, A. Sornborger, L. Cincio, and P. J. Coles, Generalization in quantum machine learning from few training data, Nature communications 13 , 4919 (2022).
- [33] J. Gibbs, Z. Holmes, M. C. Caro, N. Ezzell, H.-Y. Huang, L. Cincio, A. T. Sornborger, and P. J. Coles, Dynamical simulation via quantum machine learning with provable generalization, Physical Review Research 6 , 013241 (2024).
- [34] M. Fanizza, Y. Quek, and M. Rosati, Learning quantum processes without input control, PRX Quantum 5 , 020367 (2024).
- [35] R. Levy, D. Luo, and B. K. Clark, Classical shadows for quantum process tomography on near-term quantum computers, Physical Review Research 6 , 013029 (2024).
- [36] M. C. Caro, H.-Y. Huang, N. Ezzell, J. Gibbs, A. T. Sornborger, L. Cincio, P. J. Coles, and Z. Holmes, Outof-distribution generalization for learning quantum dynamics, Nature Communications 14 , 3751 (2023).
- [37] A. Bakshi, A. Liu, A. Moitra, and E. Tang, Structure

learning of hamiltonians from real-time evolution, arXiv preprint arXiv:2405.00082 (2024).

- [38] E. Onorati, C. Rouz´ e, D. S. Fran¸a, and J. D. Watson, c Provably efficient learning of phases of matter via dissipative evolutions, arXiv preprint arXiv:2311.07506 (2023).
- [39] H.-Y. Huang, Y. Tong, D. Fang, and Y. Su, Learning many-body hamiltonians with heisenberg-limited scaling, Physical Review Letters 130 , 200403 (2023).
- [40] A. Weiße and H. Fehske, Exact diagonalization techniques, in Computational Many-Particle Physics , edited by H. Fehske, R. Schneider, and A. Weiße (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 529-544.
- [41] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Physical review 136 , B864 (1964).
- [42] S. R. White, Density-matrix algorithms for quantum renormalization groups, Physical review b 48 , 10345 (1993).
- [43] D. Luo and B. K. Clark, Backflow transformations via neural networks for quantum many-body wave functions, Physical review letters 122 , 226401 (2019).
- [44] G. Carleo and M. Troyer, Solving the quantum manybody problem with artificial neural networks, Science 355 , 602 (2017).
- [45] D. Luo, G. Carleo, B. K. Clark, and J. Stokes, Gauge equivariant neural networks for quantum lattice gauge theories, Physical review letters 127 , 276402 (2021).
- [46] D. Luo, Z. Chen, K. Hu, Z. Zhao, V. M. Hur, and B. K. Clark, Gauge-invariant and anyonic-symmetric autoregressive neural network for quantum lattice models, Physical Review Research 5 , 013216 (2023).
- [47] D. Luo, Z. Chen, J. Carrasquilla, and B. K. Clark, Autoregressive neural network for simulating open quantum systems via a probabilistic formulation, Physical review letters 128 , 090501 (2022).
- [48] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Physical Review A 92 , 042303 (2015).
- [49] J. R. McClean, J. Romero, R. Babbush, and A. AspuruGuzik, The theory of variational hybrid quantumclassical algorithms, New Journal of Physics 18 , 023023 (2016).
- [50] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'brien, A variational eigenvalue solver on a photonic quantum processor, Nature communications 5 , 4213 (2014).
- [51] U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Optimal quantum phase estimation, Physical review letters 102 , 040403 (2009).
- [52] T. Holstein, Studies of polaron motion: Part i. the molecular-crystal model, Annals of Physics 8 , 325 (1959).
- [53] R. M. Noack, D. J. Scalapino, and R. T. Scalettar, Charge-density-wave and pairing susceptibilities in a twodimensional electron-phonon model, Phys. Rev. Lett. 66 , 778 (1991).
- [54] J. Bonˇ ca, S. A. Trugman, and I. Batisti´, Holstein poc laron, Phys. Rev. B 60 , 1633 (1999).
- [55] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109 , 1492 (1958).
- [56] P. L. DeVries and J. E. Hasbun, A first course in computational physics (Jones &amp; Bartlett Publishers, 2011).
- [57] E. W. Weisstein, Complete elliptic integral of the first kind (2024).
- [58] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel,

B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and ´ douard Duchesnay, E Scikit-learn: Machine learning in python, Journal of Machine Learning Research 12 , 2825 (2011).

## Supplemental Material for Provably Efficient Machine Learning for Adiabatic Quantum-Classical Dynamics

## I. Dimensionless Model

Consider the standard spinless Holstein model:

$$H = - t _ { n n } \sum _ { \langle i j \rangle } c _ { i } ^ { \dagger } c _ { j } - g \sum _ { i } \left ( c _ { i } ^ { \dagger } c _ { i } - \frac { 1 } { 2 } \right ) Q _ { i } + \sum _ { i } \left ( \frac { P _ { i } ^ { 2 } } { 2 M } + \frac { k Q _ { i } ^ { 2 } } { 2 } \right ).$$

And the equation of motion for the phonons is

$$\frac { d Q _ { i } } { d t } = \frac { P _ { i } } { M }, \quad \frac { d P _ { i } } { d t } = g n _ { i } - k Q _ { i },$$

where n i = ⟨ c c † i i ⟩ is the on-site fermion number. The mass m and elastic constant k are related by the familiar formula,

$$\omega & = \sqrt { \frac { k } { M } },$$

The inverse ω -1 gives a characteristic time scale for the dynamical problem. Next, one can introduce a 'length scale' Q 0 for the displacement of the simple harmonic oscillator. The energy related to Q at a given site is

$$E ( Q ) = - g n _ { i } Q _ { i } + \frac { k Q ^ { 2 } } { 2 }.$$

Assuming electron number n ∼ 1, minimization with respect to Q gives ∂E/∂Q | Q 0 = 0:

$$Q _ { 0 } & \sim \frac { g } { k }.$$

From this one can then introduce a scale for the momentum via the relation dQ/dt = P/M

$$\omega Q _ { 0 } \sim \frac { P _ { 0 } } { M } \implies P _ { 0 } = M \omega Q _ { 0 } = \frac { M \omega g } { k }.$$

We can now define the dimensionless time, displacement and momentum as

$$\tilde { t } = \omega t, \ \tilde { Q } _ { i } = \frac { Q _ { i } } { Q _ { 0 } }, \ \tilde { P } _ { i } = \frac { P _ { i } } { P _ { 0 } }.$$

In terms of dimensionless quantities, the equation of motion is then simplified to

$$\frac { d \tilde { Q } _ { i } } { d \tilde { t } } = \tilde { P } _ { i }, \quad \frac { d \tilde { P } _ { i } } { d \tilde { t } } = n _ { i } - \tilde { Q } _ { i }.$$

Next, we consider the tight-binding Hamiltonian for the fermions for a given { Q i } configuration. We factor out the nearest-neighbor hopping constant t nn and use it as the unit for energy. Also, we use the dimensionless ˜ : Q i

$$H _ { T B } = t _ { n n } \left [ - \sum _ { \langle i j \rangle } c _ { i } ^ { \dagger } c _ { j } + \frac { g Q _ { 0 } } { t _ { n n } } \sum _ { i } \tilde { Q } _ { i } c _ { i } ^ { \dagger } c _ { i } \right ].$$

The coefficient of the second term above gives an important dimensionless parameter for Holstein model. Instead of t nn , we can introduce the bandwidth of the tight-binding model: W = 4 t nn for the 1D model. We then define a dimensionless electron-phonon coupling

$$\lambda = \frac { g Q _ { 0 } } { W } = \frac { g ^ { 2 } } { k W }.$$

The dimensionless tight-binding Hamiltonian then becomes

$$\tilde { H } _ { T B } = - \sum _ { \langle i j \rangle } c _ { i } ^ { \dagger } c _ { j } + 4 \lambda \sum _ { i } \tilde { Q } _ { i } c _ { i } ^ { \dagger } c _ { i }.$$

One can see that, using these dimensionless quantities, the only adjustable parameter of the adiabatic dynamics of the Holstein model is this dimensionless λ . In general, for real materials λ ≲ 1. For example, we can set it to λ = 0 5 . or 0 25 in the simulations. .

## II. Holstein model charge density wave response analysis

Consider a tight binding model on a 1D lattice:

$$H = - t \sum _ { i } ( c _ { i } ^ { \dagger } c _ { i + 1 } + h. c. ) - g \sum _ { i } ( - 1 ) ^ { i } Q c _ { i } ^ { \dagger } c _ { i },$$

where i ∈ { 0 1 , , . . . , L - } 1 , L is the lattice size and is even. The lattice has periodic boundary condition. The effective potential on the lattice is staggered: [ -gQ, + gQ, -gQ, + gQ,.. . ]. The unit cell consists of 2 lattice sites.

We can solve the single-particle wavefunctions using the ansatz

$$| \psi _ { k } \rangle = ( a, b e ^ { i k }, a e ^ { 2 i k }, b e ^ { 3 i k }, \dots ) ^ { \mathrm T }, \ \ k \in \{ 0, \frac { 2 \pi } { L }, \dots, \frac { 2 \pi } { L } ( \frac { L } { 2 } - 1 ) \}.$$

Note that the range of k is halved because the unit cell is doubled.

The eigenvalue equation H ψ | k ⟩ = E ψ | k ⟩ becomes

$$- g Q a - 2 ( \cos k ) b & = E a, & ( \mathbb { S } 1 4 ) \\ - ( \mathfrak { P } \infty \in k \lambda \perp \alpha \Delta k & = E h & ( \mathbb { S } 1 \mathfrak { K } )$$

$$- ( 2 \cos k ) a + g Q b = E b,$$

which requires E = ± √ ( gQ ) 2 +(2cos k ) 2 to have nontrivial solutions.

The solutions form two bands. At half filling, all the states in the lower band, which has negative E , are filled with a particle. The charge density wave amplitude for | ψ k ⟩ is n k = ( a 2 -b 2 ) / 2. From the eigenvalue equation and the normalization ( L/ 2)( a 2 + b 2 ) = 1, we can solve the CDW amplitude for | ψ k ⟩ :

$$n _ { k } = \frac { \text{sgn} ( g Q ) } { L \sqrt { ( \frac { 2 \cos k } { g Q } ) ^ { 2 } + 1 } },$$

where sgn( gQ ) is the sign of gQ .

The total charge density wave at half filling is

$$n = \sum _ { k } n _ { k } = \sum _ { k } \frac { \text{sgn} ( g Q ) } { L \sqrt { ( \frac { 2 \cos k } { g Q } ) ^ { 2 } + 1 } }.$$

The range of k has been described above.

In the infinite lattice limit, L →∞ , the sum over k turns into an integral

$$n \to \int _ { 0 } ^ { \pi } \frac { d k } { 2 \pi } \frac { \text{sgn} ( g Q ) } { \sqrt { ( \frac { 2 \cos k } { g Q } ) ^ { 2 } + 1 } } = \text{sgn} ( g Q ) \frac { \text{EllipticK} [ - ( \frac { 2 } { g Q } ) ^ { 2 } ] } { \pi },$$

where EllipticK[ m ] ≡ π 2 2 F 1 ( 1 2 , 1 2 ; 1; m ) is the complete elliptic integral of the first kind, 2 F 1 ( a, b ; c ; x ) is the hypergeometric function [57].

When gQ → 0, the result says n → 0. However, the derivative ∂n ∂ gQ ( ) diverges at zero as log ( 1 gQ ) . More precisely, n ∼ -gQ log( gQ ) for small gQ .

The property of n as a function of gQ is important, especially for gQ close to zero. This is because when the system stabilizes, the equation of motion tells us k spring Q = gn , i.e. the forces are balanced. This means a straight

FIG. S1: The charge density wave amplitude n as a function of the staggering potential gQ . For L = 4 N , the function has a discontinuity at zero: There is a constant CDW = 1 L even when g → 0. For L = 4 N +2, the function is continuous and has a finite slope at zero: There is a critical slope for which a straight line cross the origin could have other crossing points with the function. i.e. There exists a critical g value, below which the system does not have a stable CDW configuration. For L →∞ (dashed curve), the function is continuous and has a logarithmically diverging slope at zero: There is always CDW, but the amplitude is exponentially small for small g .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-8601cb36-bf95-4dfd-a1e7-52b02366ad44.jpg)

line n = ( k spring /g 2 ) gQ that crosses the origin. The number of crossing points between this straight line and the n vs. gQ function determines whether the system have a stable CDW in the long time. Figure S1 shows the function curves for different lattice size.

For example, in the infinite lattice limit, the derivative diverges at zero. Therefore, any straight line crossing the origin with a finite slope will cross the function at some other points. This means that there always exists a stable CDW for the infinite lattice. However, we can estimate how large the CDW amplitude should be for a small g . Combine the straight line n = ( k spring /g 2 ) gQ and the asymptotic behavior n ∼ -gQ log( gQ ), we can solve a nontrivial crossing point at the phonon amplitude Q ∼ e -k spring /g 2 /g . We see that the amplitude decays exponentially when g → 0.

However, for a finite system size L , the situation is different. For L = 4 N , the function n has a finite limit ± 1 L when gQ → 0, and has a discontinuity at 0. Therefore, a straight line crossing the origin always has non-trivial intersections with the function curve. i.e. There exists CDW with amplitude at least 1 L for all values of g . For L = 4 N +2, the function n goes to zero when gQ → 0, and has a finite slope at 0. When g is small, ( k spring /g 2 ) exceeds this slope, the straight line only has intersection with the function curve at zero. When g is large, ( k spring /g 2 ) smaller than this slope, there are non-trivial intersections. Therefore, there is a phase transition due to the finite system size: When g &lt; g crit , no CDW; when g &gt; g crit , there is CDW. The system size we study, L = 50, is in this situation.

We can analytically solve g crit . We can take the derivative of Eq. S18 at the limit to zero:

$$\frac { \partial n } { \partial ( g Q ) } \Big | _ { g Q \to 0 ^ { + } } = \sum _ { k } \frac { 1 } { 2 L | \cos k | }.$$

By identifying ( k spring /g 2 crit ) with this slope, we get

$$g _ { c r i t } = \sqrt { \frac { 2 k _ { \text{spring} } L } { \sum _ { k \ \overline { | \cos k | } } ^ { 1 } } }.$$

Some numerical values for k spring = 1 are listed in Table. S1

For our choice L = 50, g crit = 1 1524. . The domain wall formation g value for L = 50 is around 1 6. . Therefore, we pick g ∈ [1 3 . , 1 4] in our numerical experiment, to avoid the finite size effect and the domain walls. .

In the main text, we see the bounded error requires k spring -g ∂n ∂Q &gt; 0. We can prove this is true when the system

TABLE S1: Table of slope (Eq. S19) and g crit values (Eq. S20) for different L . We choose k spring = 1.

| L      | 2      | 6        | 10       | 22       | 50       | 102      | 1002     | L →∞        |
|--------|--------|----------|----------|----------|----------|----------|----------|-------------|
| slope  | 0 . 25 | 0 . 4167 | 0 . 4972 | 0 . 6224 | 0 . 7529 | 0 . 8664 | 1 . 2300 | ∼ log( L )  |
| g crit | 2      | 1 . 5492 | 1 . 4182 | 1 . 2676 | 1 . 1524 | 1 . 0743 | 0 . 9017 | ∼ 1 √ log L |

̸

is close to the stable CDW configuration. We apply the self-consistency equation k spring = gn/Q . From the concavity of the n vs. gQ function when gQ &gt; 0, and the convexity when gQ &lt; 0, we see n gQ &gt; ∂n ∂ gQ ( ) as long as gQ = 0. Therefore, k spring -g ∂n ∂Q = g 2 ( n gQ -∂n ∂ gQ ( ) ) &gt; 0 almost surely, because the gQ = 0 case has measure 0. This completes our proof.

## III. From the EOM of Q P , to the EOM of q , p

## Approximately constant linear model

In the main text, we introduced the 'approximately constant linear model', which has the Hamiltonian:

$$\hat { H } = \hat { H } _ { q } + \sum _ { i } g \hat { O } _ { i } Q _ { i } + \sum _ { i } \left ( \frac { 1 } { 2 M } P _ { i } ^ { 2 } + \frac { 1 } { 2 } k Q _ { i } ^ { 2 } \right ).$$

It has the EOM for Q i and P i :

$$\frac { d } { d t } Q _ { i } ( t ) = \frac { 1 } { M } P _ { i } ( t ),$$

$$\frac { \overset { \cdot \cdot } { d } } { d t } P _ { i } ( t ) = - g \langle \hat { O } _ { i } \rangle - k Q _ { i } ( t ) - \gamma P _ { i } ( t ),$$

where the quantum observable ⟨ ˆ O i ⟩ is computed based on the configuration ⃗ Q t ( ).

The EOM above is equally valid for both the exact dynamical path and the ML dynamical path:

$$\frac { d } { d t } Q _ { i, \text{Exact} } ( t ) = \frac { 1 } { M } P _ { i, \text{Exact} } ( t ),$$

$$\frac { d } { d t } Q _ { i, \text{Exact} } ( t ) & = \frac { 1 } { M } P _ { i, \text{Exact} } ( t ), \\ \frac { d } { d t } P _ { i, \text{Exact} } ( t ) & = - g ( \hat { O } _ { i } ) _ { \text{Exact} } - k Q _ { i, \text{Exact} } ( t ) - \gamma P _ { i, \text{Exact} } ( t ),$$

and

$$\frac { d } { d t } Q _ { i, \text{ML} } ( t ) = \frac { 1 } { M } P _ { i, \text{ML} } ( t ),$$

$$\frac { \cdots } { d t } P _ { i, M L } ( t ) = - g \langle \hat { O } _ { i } \rangle _ { M L } - k Q _ { i, M L } ( t ) - \gamma P _ { i, M L } ( t ).$$

Take the difference between the ML EOM and the exact EOM, define the accumulated position and momentum errors q t ( ) = Q i, ML ( ) t -Q i, Exact ( ), t p t ( ) = P i, ML ( ) t -P i, Exact ( ), we get t

d

1

dt q t

( ) =

M p t ,

( )

(S28)

$$\frac { \ddots } { d t } p ( t ) & = - g ( \langle \hat { O } _ { i } \rangle _ { \text{ML} } - \langle \hat { O } _ { i } \rangle _ { \text{Exact} } ) - k q ( t ) - \gamma p ( t ).$$

However, here the observables depend on different configurations, ⟨ ˆ O i ⟩ ML = ⟨ ˆ O i ⟩ ML ( ⃗ Q ML ( )), t and ⟨ ˆ O i ⟩ Exact = ⟨ ˆ O i ⟩ Exact ( ⃗ Q Exact ( )). t When we take their difference, we not only need to take care about the difference between the prediction methods, we also need to take care about the difference between the configurations ⃗ Q ML and ⃗ Q Exact .

We do a Taylor expansion and apply the assumption that the off-diagonal response ∂ O ⟨ ˆ i ⟩ ∂Q j ( i = ) is approximately j

̸

zero. We can make the difference between the observables into two terms:

$$\langle \hat { O } _ { i } \rangle _ { \text{ML} } ( \vec { Q } _ { \text{ML} } ( t ) ) - \langle \hat { O } _ { i } \rangle _ { \text{Exact} } ( \vec { Q } _ { \text{Exact} } ( t ) ) & = \langle \hat { O } _ { i } \rangle _ { \text{ML} } ( \vec { Q } _ { \text{ML} } ( t ) ) - \langle \hat { O } _ { i } \rangle _ { \text{Exact} } ( \vec { Q } _ { \text{ML} } ( t ) ) + \frac { \partial \langle \hat { O } _ { i } \rangle } { \partial Q _ { i } } q ( t ) + o ( q ( t ) ) \\ & = \delta \langle \hat { O } _ { i } \rangle ( t ) + \frac { \partial \langle \hat { O } _ { i } \rangle } { \partial Q _ { i } } q ( t ) + o ( q ( t ) ),$$

where δ O ⟨ ˆ i ⟩ ( ) = t ⟨ ˆ O i ⟩ ML ( ⃗ Q ML ( )) t -⟨ ˆ O i ⟩ Exact ( ⃗ Q ML ( )) is the single-step prediction error. t Insert the two terms back to the EOM, we get the same result as in the main text:

$$\frac { d } { d t } q ( t ) = \frac { 1 } { M } p ( t ),$$

$$\frac { d } { d t } q ( t ) & = \frac { 1 } { M } p ( t ), \\ \frac { d } { d t } p ( t ) & = - g \delta \langle \hat { O } _ { i } \rangle ( t ) - g \frac { \partial \langle \hat { O } _ { i } \rangle } { \partial Q _ { i } } q ( t ) - k q ( t ) - \gamma p ( t ) + o ( q ( t ) ).$$

## Non-linear model

In the main text, we then generalized a bit to the non-linear model, which has the Hamiltonian:

$$\hat { H } = \hat { H } _ { q } + \sum _ { i } g \hat { O } _ { i } G _ { i } ( Q _ { i } ) + \sum _ { i } \left ( \frac { 1 } { 2 M } P _ { i } ^ { 2 } + V _ { i } ( Q _ { i } ) \right ).$$

It has the EOM for Q i and P i :

and

$$\frac { d } { d t } Q _ { i, M L } ( t ) & = \frac { 1 } { M } P _ { i, M L } ( t ), \\ \frac { d } { d t } P _ { i, M L } ( t ) & = - g \langle \hat { O } _ { i } \rangle _ { M L } G ^ { ^ { \prime } } _ { i } ( Q _ { i, M L } ) - V ^ { ^ { \prime } } _ { i } ( Q _ { i, M L } ) - \gamma P _ { i, M L } ( t ).$$

$$\frac { a } { d t } P _ { i, \text{ML} } ( t ) = - g \langle \hat { O } _ { i } \rangle _ { \text{ML} } G ^ { ^ { \prime } } _ { i } ( Q _ { i, \text{ML} } ) - V ^ { ^ { \prime } } _ { i } ( Q _ { i, \text{ML} } ) - \gamma P _ { i, \text{ML} } ( t$$

Take the difference between the ML EOM and the exact EOM, define the accumulated position and momentum errors q t ( ) = Q i, ML ( ) t -Q i, Exact ( ), t p t ( ) = P i, ML ( ) t -P i, Exact ( ), expand to first order of t q t ( ), we get:

$$\frac { d } { d t } q ( t ) & = \frac { 1 } { M } p ( t ), \\ \frac { d } { d t } p ( t ) & = - g ( \langle \hat { O } _ { i } \rangle _ { M L } - \langle \hat { O } _ { i } \rangle _ { \text{Exact} } ) G _ { i } ^ { ^ { \prime } } ( Q _ { i, M L } ) - g ( \hat { O } _ { i } \rangle _ { \text{Exact} } G _ { i } ^ { ^ { \prime \prime } } ( Q _ { i, \text{Exact} } ) q ( t ) - V _ { i } ^ { ^ { \prime \prime } } ( Q _ { i, \text{Exact} } ) q ( t ) - \gamma p ( t ) + o ( q ( t ) ).$$

We apply the same technique as in the approximately constant linear model. Insert

$$\langle \hat { O } _ { i } \rangle _ { M L } - \langle \hat { O } _ { i } \rangle _ { \text{Exact} } = \delta \langle \hat { O } _ { i } \rangle ( t ) + \frac { \partial \langle \hat { O } _ { i } \rangle } { \partial Q _ { i } } q ( t ) + o ( q ( t ) ),$$

$$\frac { d } { d t } Q _ { i } ( t ) = \frac { 1 } { M } P _ { i } ( t ),$$

$$\frac { \overset { \cdot \cdot } { d } } { d t } P _ { i } ( t ) = - g ( \hat { O } _ { i } ) G _ { i } ^ { ^ { \prime } } ( Q _ { i } ) - V _ { i } ^ { ^ { \prime } } ( Q _ { i } ) - \gamma P _ { i } ( t ),$$

where the quantum observable ˆ O i is still computed based on the configuration ⃗ Q t

The EOM above is equally valid for both the exact dynamical path and the ML dynamical path:

⟨ ⟩ ( ).

$$\frac { d } { d t } Q _ { i, \text{Exact} } ( t ) & = \frac { 1 } { M } P _ { i, \text{Exact} } ( t ), \\ \frac { d } { d t } P _ { i, \text{Exact} } ( t ) & = - g ( \hat { O } _ { i } ) _ { \text{Exact} } G _ { i } ^ { ^ { \prime } } ( Q _ { i, \text{Exact} } ) - V _ { i } ^ { ^ { \prime } } ( Q _ { i, \text{Exact} } ) - \gamma P _ { i, \text{Exact} } ( t ),$$

we get

$$\frac { d } { d t } q ( t ) = \frac { 1 } { M } p ( t ), \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$\frac { u u } { d t } p ( t ) & = - g \delta ( \hat { O } _ { i } ) ( t ) G _ { i } ^ { \prime } ( Q _ { i, M } ) - g \frac { \partial ( \hat { O } _ { i } ) } { \partial Q _ { i } } G _ { i } ^ { \prime } ( Q _ { i }, M _ { i } ) q ( t ) - g ( \hat { O } _ { i } ) _ { \text{Exact} } G _ { i } ^ { \prime \prime } ( Q _ { i, \text{Exact} } ) q ( t ) - V _ { i } ^ { \prime \prime } ( Q _ { i, \text{Exact} } ) q ( t ) - \gamma p ( t ) + o ( q ( t ) ). \\ & = F ( t ) - K ( t ) q ( t ) - \gamma p ( t ) + o ( q ( t ) ),$$

$$t ) - K ( t ) q ( t ) - \gamma _ { 1 }$$

where we define the effective driving force and the effective spring constant as

$$F ( t ) = - g \delta \langle \hat { O } _ { i } \rangle ( t ) G _ { i } ^ { ^ { \prime } } ( Q _ { i, \text{ML} } ( t ) ), \\ \alpha / \hat { \wedge } \, \real$$

$$K ( t ) = g \frac { \partial ( \hat { O } _ { i } ) } { \partial Q _ { i } } G _ { i } ^ { ^ { \prime } } ( Q _ { i, \text{ML} } ( t ) ) + g ( \hat { O } _ { i } ) _ { \text{Exact} } G _ { i } ^ { ^ { \prime \prime } } ( Q _ { i, \text{Exact} } ( t ) ) + V _ { i } ^ { ^ { \prime \prime } } ( Q _ { i, \text{Exact} } ( t ) ).$$

This EOM is equivalent to a damped spring with time dependent driving force and spring constant. We can prove the amplitude of the spring will not diverge if F t ( ) is bounded and K t ( ) does not fluctuate largely.

## Generic model

The generic model has the Hamiltonian

$$\hat { H } = \hat { H } _ { q } + \sum _ { \alpha, i } g _ { \alpha } \hat { O } _ { \alpha, i } G _ { \alpha, i } ( \vec { P }, \vec { Q } ) + H _ { c l } ( \vec { P }, \vec { Q } ),$$

and the EOM

$$\frac { d } { d t } Q _ { j } = \sum _ { \alpha, i } g _ { \alpha } \langle \hat { O } _ { \alpha, i } \rangle \frac { \partial } { \partial P _ { j } } G _ { \alpha, i } + \frac { \partial } { \partial P _ { j } } H _ { c l },$$

$$\frac { d } { d t } P _ { j } = - \sum _ { \alpha, i } g _ { \alpha } \langle \hat { O } _ { \alpha, i } \rangle \frac { \partial } { \partial Q _ { j } } G _ { \alpha, i } - \frac { \partial } { \partial Q _ { j } } H _ { c l } - \gamma P _ { j }.$$

We can repeat the techniques above: writing down the EOM for ML and the exact simulation, making a difference of the EOM, defining q i and p i (here the index i cannot be suppressed), doing Taylor expansion, and combining terms up to the first order. After similar derivations, we can get the result

$$\frac { d } { d t } q _ { i } = \sum _ { j } \left [ q _ { j } \mathcal { K } _ { q _ { j }, q _ { i } } ( \vec { P }, \vec { Q } ) + p _ { j } \mathcal { K } _ { p _ { j }, q _ { i } } ( \vec { P }, \vec { Q } ) \right ] + \mathcal { F } _ { q _ { i } } + o ( q, p ),$$

$$\frac { d } { d t } p _ { i } = \sum _ { j } ^ { \cdot } \left [ q _ { j } \mathcal { K } _ { q _ { j }, p _ { i } } ( \vec { P }, \vec { Q } ) + p _ { j } \mathcal { K } _ { p _ { j }, p _ { i } } ( \vec { P }, \vec { Q } ) \right ] + \mathcal { F } _ { p _ { i } } - \gamma p _ { i } + o ( q, p ),$$

where

$$\mathcal { K } _ { q _ { j }, q _ { i } } = \sum _ { \alpha, k } g _ { \alpha } \left ( \frac { \partial \langle \hat { O } _ { \alpha, k } \rangle } { \partial Q _ { j } } \frac { \partial G _ { \alpha, k } } { \partial P _ { i } } + \langle \hat { O } _ { \alpha, k } \rangle \frac { \partial ^ { 2 } G _ { \alpha, k } } { \partial Q _ { j } \partial P _ { i } } \right ) + \frac { \partial ^ { 2 } H _ { c l } } { \partial Q _ { j } \partial P _ { i } },$$

$$\mathcal { K } _ { p _ { j }, q _ { i } } = \sum _ { \alpha, k } g _ { \alpha } \left ( \frac { \partial \langle \hat { O } _ { \alpha, k } \rangle } { \partial P _ { j } } \frac { \partial G _ { \alpha, k } } { \partial P _ { i } } + \langle \hat { O } _ { \alpha, k } \rangle \frac { \partial ^ { 2 } G _ { \alpha, k } } { \partial P _ { j } \partial P _ { i } } \right ) + \frac { \partial ^ { 2 } H _ { c l } } { \partial P _ { j } \partial P _ { i } },$$

$$\mathcal { K } _ { q _ { j }, p _ { i } } = - \sum _ { \alpha, k } g _ { \alpha } \left ( \frac { \partial \langle \hat { O } _ { \alpha, k } \rangle } { \partial Q _ { j } } \frac { \partial G _ { \alpha, k } } { \partial Q _ { i } } + \langle \hat { O } _ { \alpha, k } \rangle \frac { \partial ^ { 2 } G _ { \alpha, k } } { \partial Q _ { j } \partial Q _ { i } } \right ) - \frac { \partial ^ { 2 } H _ { c l } } { \partial Q _ { j } \partial Q _ { i } },$$

$$\mathcal { K } _ { p _ { j }, p _ { i } } = - \sum _ { \alpha, k } g _ { \alpha } \left ( \frac { \partial \langle \hat { O } _ { \alpha, k } \rangle } { \partial P _ { j } } \frac { \partial G _ { \alpha, k } } { \partial Q _ { i } } + \langle \hat { O } _ { \alpha, k } \rangle \frac { \partial ^ { 2 } G _ { \alpha, k } } { \partial P _ { j } \partial Q _ { i } } \right ) - \frac { \partial ^ { 2 } H _ { c l } } { \partial P _ { j } \partial Q _ { i } },$$

$$\mathcal { F } _ { q _ { i } } = \sum _ { \alpha, j } g _ { \alpha } \delta \langle \hat { O } _ { \alpha, j } \rangle ( t ) \frac { \partial } { \partial P _ { i } } G _ { \alpha, j } ( \vec { P }, \vec { Q } ),$$

$$\mathcal { F } _ { p _ { i } } = - \sum _ { \alpha, j } g _ { \alpha } \delta \langle \hat { O } _ { \alpha, j } \rangle ( t ) \frac { \partial } { \partial Q _ { i } } G _ { \alpha, j } ( \vec { P }, \vec { Q } ).$$

## IV. Proof of Proposition 1

Here we prove that the approximately constant linear model satisfies the Error Bounded Property in Def. 1 if K &gt; 0. We use the analogy to a damped harmonic oscillating spring, phrasing the error force F t ( ) in the main text as the 'driving force' and the error stiffness K in the main text as the 'spring constant'. For a damped spring, given a bounded driving force, even if the force is tuned to drive the spring optimally, as long as K,γ &gt; 0, the spring cannot be driven to infinite amplitude. This can be seen by the following worst-case analysis. In a worst case scenario, the driving force is set to be the maximal value in the direction of the spring movement. We know a constant driving force means a shift of the reference point of the spring. Therefore, by shifting the reference point back and forth, the driving force at most linearly increases the amplitude with the motion cycles of the spring. However, a finite damping decreases the amplitude of the spring by a constant factor in each cycle. Therefore, for a large enough initial amplitude, the decrement of the amplitude due to the damping must exceed the increment of the amplitude due to the driving force in the cycle, and thus the amplitude in the next cycle must be smaller than the initial one, which tells us that the damped spring cannot be driven to infinite amplitude.

Now we analyze the relation between the maximal amplitude ¯, the maximal momentum ¯, and the maximal driving q p force ¯ , using the dimensional analysis. F At t init , there are no accumulated errors, i.e. q t ( init ) = p t ( init ) = 0. Therefore, the only physical quantity that carries the dimension of length is the maximal driving force [ ¯ ] = [ F M L T ][ ][ ] -2 , where [ M ], [ L ], and [ T ] are the dimension of mass, length, and time, respectively. Note that the maximal amplitude [¯] = [ q L ] and the maximal momentum ¯ = [ p L T ][ ] -1 both carries power one of the length dimension [ L ]. Therefore, by dimensional analysis, they are both proportional to the maximal driving force ¯ . F i.e. ∃ ¯ C , C q ¯ p such that ¯ = q ¯ C F q ¯ , ¯ = p ¯ C F p ¯ .

For any ϵ &gt; 0, if | δ O ⟨ ˆ i ⟩ ( ) t | 2 ≤ ϵ , by the definition F t ( ) = -gδ O ⟨ ˆ i ⟩ ( ), we have t | F t ( ) | ≤ | g | √ ϵ . By the analogy above, it means that the maximal force ¯ = F | g | √ ϵ . Therefore, by the meaning of the maximal amplitude and the maximal momentum, there are | q t ( ) | ≤ ¯ = q ¯ C F q ¯ = ¯ C q | g | √ ϵ and | p t ( ) | ≤ √ ¯ = p ¯ C F p ¯ = ¯ C p | g | √ ϵ . Define two new constants C q = | g C | ¯ q and C p = | g C | ¯ p , we get | q t ( ) | ≤ C q √ ϵ and | p t ( ) | ≤ C p ϵ , which is what we want to show. This completes the proof of Prop. 1.

## V. Proof of Theorem 1 (Error Bounded Condition for non-linear model)

Here we present the formal version of Thm. 1 in the main text. We use the analogy to a damped harmonic oscillating spring, phrasing the error force F t ( ) in the main text as the 'driving force' and the error stiffness K t ( ) in the main text as the 'spring constant'.

Theorem 3 (Error Bounded Condition for non-linear model) . Consider a damped harmonic oscillator with the following EOM:

$$\frac { d } { d t } q ( t ) = \frac { 1 } { M } p ( t ),$$

$$\frac { d } { d t } p ( t ) = F ( t ) - K ( t ) q ( t ) - \gamma p ( t ) + o ( q ( t ) ),$$

where q is the classical coordinate, p is the classical momentum, t is time, M is the mass of the oscillator, F t ( ) is a time-dependent driving force, K t ( ) is a time-dependent spring constant, γ is the damping coefficient.The initial condition is q t ( init ) = p t ( init ) = 0 .

Given K max &gt; K min &gt; M γ/ ( 2) 2 satisfying the following inequality:

$$\frac { K _ { \max } } { K _ { \min } } < \exp \left [ 2 \left ( \frac { a r c t a n \, \omega _ { \min } } { \omega _ { \min } } + \frac { \pi - a r c t a n \, \omega _ { \max } } { \omega _ { \max } } \right ) \right ],$$

where ω max = √ K max M γ/ ( 2) 2 -1 , ω min = √ K min M γ/ ( 2) 2 -1 .

If K t ( ) ∈ [ K min , K max ] for all t , and G α,i ( ⃗ , Q P ⃗ ) and its first derivatives are bounded, then the Error Bounded Property in Def. 1 is satisfied, i.e. ∃ C , C q p &gt; 0 such that ∀ ϵ &gt; 0 , if | δ O ⟨ ˆ i ⟩ ( ) t | 2 ≤ ϵ for all t , then there are

$$| q ( t ) | & \leq C _ { q } \sqrt { \epsilon }, \quad | p ( t ) | \leq C _ { p } \sqrt { \epsilon }, \quad \text{for all } t.$$

Proof. Consider a damped harmonic oscillator with the following EOM:

$$\frac { d } { d t } q ( t ) = \frac { 1 } { M } p ( t ),$$

$$\frac { d } { d t } p ( t ) = F ( t ) - K ( t ) q ( t ) - \gamma p ( t ) + o ( q ( t ) ),$$

where q is the classical coordinate, p is the classical momentum, t is time, M is the mass of the oscillator, F t ( ) is a time-dependent driving force, K t ( ) is a time-dependent spring constant, γ &gt; 0 is the damping coefficient. We assume F t ( ) ∈ [ -¯ F,F ¯ ], K t ( ) ∈ [ K min , K max ], where K max &gt; K min &gt; M γ/ ( 2) 2 . Note that with Eq. S45, F t ( ) = -gδ O ⟨ ˆ i ⟩ ( ) t G Q ′ i ( i, ML ( )), the assumption t | δ O ⟨ ˆ i ⟩ ( ) t | 2 ≤ ϵ implies ¯ F ≤ O ( √ ϵ ) given bounded G ′ i .

We consider the worst case scenario. The mass on the spring starts from one side with a large amplitude. As the mass moving towards the other side, the driving force always keeps the maximal value towards the other side. Before the mass crosses the origin, the spring constant is set as K max to maximize the drag towards the other side. After the mass crosses the origin, the spring constant is set as K min to minimize the burden for its moving as far as possible. In this worst case scenario, the mass will stop at its largest amplitude on the other side. If this amplitude on the other side is smaller than the one it started with, than the spring cannot have diverging amplitude, and therefore the error converges.

We now analyze this worst case scenario. Assume the mass starts at q ( -t 0 ) = -q 0 , ˙( q -t 0 ) = 0. At time t = 0 it crosses the origin, q (0) = 0, with a velocity ˙(0) = q v . At time t = t 1 it stops on the other side at q t ( 1 ) = q 1 , ˙( q t 1 ) = 0.

Before the mass crosses the origin, the spring constant is K max . The EOM is

$$\tilde { q } ( t ) = - \gamma \dot { q } ( t ) - \frac { K _ { \max } } { M } q ( t ) + \frac { \bar { F } } { M }, \ \ t \in [ - t _ { 0 }, 0 ].$$

The EOM has a general solution

$$q ( t ) = A e ^ { - \frac { \gamma } { 2 } t } \sin \Omega _ { \max } t + B e ^ { - \frac { \gamma } { 2 } t } \cos \Omega _ { \max } t + \frac { \bar { F } } { K _ { \max } },$$

where Ω max = √ K max M -( γ 2 ) 2 .

Inserting the boundary condition at t = 0, we get equations

$$q ( 0 ) = 0 = B + \frac { \bar { F } } { K _ { \max _ { \varepsilon } } },$$

$$\stackrel { \dots } { q } ( 0 ) = v = \Omega _ { \max } A - \frac { \gamma } { 2 } B,$$

$$A = \frac { 1 } { \Omega _ { \max } } \left ( v - \frac { \gamma } { 2 } \frac { \bar { F } } { K _ { \max } } \right ),$$

$$B = - \frac { \stackrel { \dots \dots } { F } } { K _ { \max } }.$$

Inserting the boundary condition at t = -t 0 , we get equations

$$q ( - t _ { 0 } ) = - q _ { 0 } = - A e ^ { \frac { 3 } { 2 } t _ { 0 } } \sin \Omega _ { \max } t _ { 0 } + B e ^ { \frac { 3 } { 2 } t _ { 0 } } \cos \Omega _ { \max } t _ { 0 } + \frac { \bar { F } } { K _ { \max } }, \\ \Gamma _ { \sim \, \cdot \, \cdot \, \cdot } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \lambda } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \alpha } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \lambda } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \phi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xeta } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \lambda } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xeta } \Gamma _ {\xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ {\xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xeta } \Gamma _{ \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _{ \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xeta } \Gamma } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \axis } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \chi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \sigma } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \zeta } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gamma _ { \xi } \Gam$$

$$\dot { q } ( - t _ { 0 } ) = 0 = A \left ( \frac { \gamma } { 2 } e ^ { \frac { \vec { 2 } t _ { 0 } } { 2 } } \sin \Omega _ { \max } t _ { 0 } + \Omega _ { \max } e ^ { \frac { \vec { 2 } t _ { 0 } } { 2 } } \cos \Omega _ { \max } t _ { 0 } \right ) + B \left ( \frac { \gamma } { 2 } e ^ { \frac { \vec { 2 } t _ { 0 } } { 2 } } \cos \Omega _ { \max } t _ { 0 } + \Omega _ { \max } e ^ { \frac { \vec { 2 } t _ { 0 } } { 2 } } \sin \Omega _ { \max } t _ { 0 } \right ).$$

We can solve t 0 from the second equation:

$$t _ { 0 } & = \frac { 1 } { \Omega _ { \max } } \arccot \frac { ( \gamma / 2 ) A + \Omega _ { \max } B } { - \Omega _ { \max } A + ( \gamma / 2 ) B } \\ & = \frac { 1 } { \Omega _ { \max } } \arccot \left ( - \frac { ( \gamma / 2 ) v - \bar { F } / M } { \Omega _ { \max } v } \right ),$$

where we have inserted the solution of A and B . We can then insert everything into the first equation to solve q 0 , but we do not do it now.

Let's also take a look at the time after the mass crosses the origin. The spring constant changes to K min . The EOM is

$$\ddot { q } ( t ) = - \gamma \dot { q } ( t ) - \frac { K _ { \min } } { M } q ( t ) + \frac { \bar { F } } { M }, \ \ t \in [ 0, t _ { 1 } ].$$

The EOM has a general solution from which we solve

from which we solve

$$q ( t ) = C e ^ { - \frac { 2 } { 2 } t } \sin \Omega _ { \min } t + D e ^ { - \frac { 2 } { 2 } t } \cos \Omega _ { \min } t + \frac { \bar { F } } { K _ { \min } },$$

where Ω min = √ K min M -( γ 2 ) 2 .

Inserting the boundary condition at t = 0, we get equations

$$q ( 0 ) = 0 = D + \frac { \bar { F } } { K _ { \min } },$$

$$\dot { q } ( 0 ) = v = \overset { \dots } { \Omega } _ { \min } C - \frac { \gamma } { 2 } D,$$

$$C = \frac { 1 } { \Omega _ { \min } } \left ( v - \frac { \gamma } { 2 } \frac { \bar { F } } { K _ { \min } } \right ),$$

$$D = - \frac { \overset { \dots } { F } } { K _ { \min } }.$$

Inserting the boundary condition at t = t 1 , we get equations

$$& \nearrow \quad \cdot \quad \cdot \quad \cdot \quad \cdot \\ & q ( t _ { 1 } ) = q _ { 1 } = C e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \sin \Omega _ { \min } t _ { 1 } + D e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \cos \Omega _ { \min } t _ { 1 } + \frac { \tilde { F } } { K _ { \min } }, \\ & \dot { q } ( t _ { 1 } ) = 0 = C \left ( - \frac { \gamma } { 2 } e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \sin \Omega _ { \min } t _ { 1 } + \Omega _ { \min } e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \cos \Omega _ { \min } t _ { 1 } \right ) + D \left ( - \frac { \gamma } { 2 } e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \cos \Omega _ { \min } t _ { 1 } - \Omega _ { \min } e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \sin \Omega _ { \min } t _ { 1 } \right ).$$

We can solve t 1 from the second equation:

$$t _ { 1 } & = \frac { 1 } { \Omega _ { \min } } \arccot \frac { ( \gamma / 2 ) C + \Omega _ { \min } D } { \Omega _ { \min } C - ( \gamma / 2 ) D } \\ & = \frac { 1 } { \Omega _ { \min } } \arccot \frac { ( \gamma / 2 ) v - \bar { F } / M } { \Omega _ { \min } v },$$

where we have inserted the solution of C and D . Note that comparing to the case before, there is no minus sign in the arccot function.

To have a bounded error, we want that for large enough q 0 , there is q 1 q 0 &lt; 1, where q 0 is the starting amplitude and q 1 is the stopping amplitude. i.e. The amplitude never diverges. Iteratively we consider q 1 as the next starting amplitude, and we can see the sequence of amplitudes converges in the long time. Note that the amplitude has the dimension of length, and among all given parameters, only the maximal force ¯ carries the dimension of length. F Using dimension analysis, we see that the maximal amplitude the system could reach will be proportional to the maximal force ¯ . F The maximal momentum the system could reach will also be proportional to ¯ F from dimension analysis. i.e. q t ( ) ≤ O F ( ¯ ), p t ( ) ≤ O F ( ¯ ). Combining with the assumption ¯ F ≤ O ( √ ϵ ), it leads to q t ( ) ≤ O ( √ ϵ ), p t ( ) ≤ O ( √ ϵ ), which is what we want to show.

Now the final step is to simplify the condition q 1 q 0 &lt; 1. Note that for arbitrarily large q 0 , the velocity v can also be arbitrarily large. Therefore, we can consider the limit ¯ F/v → 0. Under this limit,

$$A / v \rightarrow \frac { 1 } { \Omega _ { \max } }, \quad \dots.$$

$$B / v \rightarrow 0, \begin{matrix} B & 1 \\ 1 \end{matrix}$$

$$C / v \rightarrow \frac { 1 } { \Omega _ { \min } }, \quad \dots.$$

$$D / v \to 0, \quad & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & && & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ; & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & - & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & _ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & + & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & = & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &, & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ( & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ^ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & b & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & | & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &$$

$$t _ { 0 } \rightarrow \frac { 1 } { \Omega _ { \max } } \arccot \left ( - \frac { ( \gamma / 2 ) } { \Omega _ { \max } } \right ) = \frac { \pi - \arctan ( 2 \Omega _ { \max } / \gamma ) } { \Omega _ { \max } },$$

$$t _ { 1 } \rightarrow \frac { 1 } { \Omega _ { \min } } \arccot \frac { ( \gamma / 2 ) } { \Omega _ { \min } } = \frac { \arctan ( 2 \Omega _ { \min } / \gamma ) } { \Omega _ { \min } }.$$

Finally,

$$\frac { q _ { 1 } } { q _ { 0 } } & = \frac { C e ^ { - \frac { \gamma } { 2 } t _ { 1 } } } { q _ { 0 } } \sin \Omega _ { \min } t _ { 1 } + D e ^ { - \frac { \gamma } { 2 } t _ { 1 } } \cos \Omega _ { \min } t _ { 1 } + \bar { F } / K _ { \min } } { A e ^ { \frac { \gamma } { 2 } t _ { 0 } } \sin \Omega _ { \max } t _ { 0 } - B e ^ { \frac { \gamma } { 2 } t _ { 0 } } \cos \Omega _ { \max } t _ { 0 } - \bar { F } / K _ { \max } } \\ & \to \frac { ( 1 / \Omega _ { \min } ) e ^ { - \frac { \gamma } { 2 } t _ { 1 } ^ { * } } ( \Omega _ { \min } / \sqrt { K _ { \min } / M } ) } { ( 1 / \Omega _ { \max } ) e ^ { \frac { \gamma } { 2 } t _ { 0 } ^ { * } } ( \Omega _ { \max } / \sqrt { K _ { \max } / M } ) } \\ & = \sqrt { \frac { K _ { \max } } { K _ { \min } } } e ^ { - \frac { \gamma } { 2 } ( t _ { 0 } ^ { * } + t _ { 1 } ^ { * } ) }, \\.. \quad. \quad. \quad. \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cd.. \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cd: \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cdots \ \cd< \quad........................................................................................................................................................................................................$$

where we have used the equality: sin [ arccot ( -( γ/ 2) Ω max )] = Ω max √ K max /M , sin [ arccot ( γ/ 2) Ω min ] = Ω min √ K min /M . Here t ∗ 0 and t ∗ 1 are the limiting values in Eq. S86 and S87, respectively.

$$\text{We define $\omega_{max}=\frac{\Omega_{max}} { ( 7 / 2 ) } = \sqrt { \frac { K _ { \max } } { M ( 7 / 2 ) ^ { 2 } } - 1 } \text{ and $\omega_{min}=\frac{\Omega_{min}} { ( 7 / 2 ) } = \sqrt { \frac { K _ { \min } } { M ( 7 / 2 ) ^ { 2 } } - 1 }. \text{ After rearranging the terms, we}$$

FIG. S2: Numerical measurement of the error stiffness K t ( ) = k -g ∂n ∂Q during the time range of interest. g = 1 4 is . used in the experiment. The highest point in the plot reaches 0.7104, and the lowest is 0.5011. Inserting these values into Eq. S60, we get LHS = 1 418 . &lt; . 1 501 = RHS .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00276v2-figure-3999daa1-024d-45cc-8ac9-0e9084bb1e19.jpg)

simplifies the condition q 1 q 0 &lt; 1 into

$$\frac { K _ { \max } } { K _ { \min } } < \exp \left [ 2 \left ( \frac { \arctan \omega _ { \min } } { \omega _ { \min } } + \frac { \pi - \arctan \omega _ { \max } } { \omega _ { \max } } \right ) \right ],$$

which is the inequality we see in Eq. S60 in Thm. 3.

This completes the proof of Thm. 3.

Note that the inequality in Eq. S60 in Thm. 3 is satisfied in the limit K max -K min K min → 0. We can see this by take a log on both sides of the inequality and get

$$\log \left ( 1 + \frac { K _ { \max } - K _ { \min } } { K _ { \min } } \right ) < 2 \left ( \frac { \arctan \omega _ { \min } } { \omega _ { \min } } + \frac { \pi - \arctan \omega _ { \max } } { \omega _ { \max } } \right ),$$

whose left hand side is O ( K max -K min K min ) and right hand side is O (1) because ω min ≈ ω max . Therefore, in the limit K max -K min K min → 0, the inequality always holds. This agrees with the Error Bounded Property of the approximately constant linear model stated in Prop. 1. (Actually the approximately constant linear model satisfies the Error Bounded Property under a more general condition: instead of K &gt; M γ/ ( 2) 2 we only need K &gt; 0 for the approximately constant linear model.)

If we transform the coordinates of the system to their normal modes, because the normal modes are also described by damped harmonic oscillators, all the analysis works the same for normal modes. This addresses the comment under Thm. 1 in the main text.

In Figure S2, we show an example of the fluctuation of the error stiffness K t ( ) = k -g ∂n ∂Q during the Holstein model dynamics. We see that the error stiffness K t ( ) &gt; M γ/ ( 2) 2 = 0 0025, where . M = 1 and γ = 0 1. . The highest and lowest points in the plot are 0.7104 and 0.5011. We insert these numbers into Eq. S60. The left hand side of the condition results in 1.418, and the right hand side 1.501. Therefore, we confirm that the inequality in Eq. S60 in Thm. 3 is satisfied during the time range we are interested in, and thus the Error Bounded Property is satisfied in the Holstein model.

## VI. Provably efficient adiabatic learning

Here we provide more details about how the PEAL algorithm is performed.

Step I: Collecting data. To be able to train a learning model, we first need to collect data samples from the distribution we are interested in. For example, one can sample a few different initial conditions, and use either classical solver or quantum computer to obtain a few dynamical evolution paths. The 'configuration'-'observable' pairs sampled from those paths can serve as the data set we need. If one wants to predict other k-local observable that is not involved in the dynamics, one can compute that observable along the paths and collect the data for the next training step as well.

Step II: ML training. A provable efficient ML algorithm for predicting ground state properties was presented recently [27]. Here we apply this algorithm to adiabatic dynamics. On the data set collected in Step I, a nonlinear feature map with geometrically local region information is performed. A model is trained with an l 1 -regularized regression (LASSO) on the features. Hyperparameters in the model are properly chosen. One can train a single model for one type of observable if the system has translation symmetry. The same training process can be applied to the data set of observable that is not involved in the dynamics. Moreover, the learning model can be transferred to unseen physical parameter g α in the Hamiltonian. An example is shown in our numerical experiment, also illustrated in Fig. 2 in the main text.

Step III: ML prediction. For dynamical evolution paths with unseen initial conditions, we can use the PEAL prediction to speed up the simulation, instead of repeatedly using costly classical solver or quantum computer. The dynamical simulation is realized by iteratively updating the classical degree of freedom and the quantum degree of freedom. The classical updating step is done by classical ODE solver. The quantum updating step is done by ML prediction. In each iteration, we can add an extra correction step to achieve symmetry-preserving PEAL (see details in the following section). The provably efficiency and controllable error of the PEAL algorithm will be presented later in Thm. 6.

Before we formally present Thm. 6, we would like to mention two other theorems introduced in other works. In Ref. [27], the authors prove two theorems:

Theorem 4. Consider any family of n-qubit geometrically-local Hamiltonians { H x ( ) : x ∈ -[ 1 1] , m } in a finite spatial dimension, such that each local term in H x ( ) depends smoothly on x , and the smallest eigenvalue and the next smallest eigenvalues have a constant gap ≥ Ω(1) between them. Then the ground state properties can be efficiently predicted.

Theorem 5. Given n, δ &gt; 0 , 1 e &gt; ϵ &gt; 0 and a training data set { x , y t t } N t =1 of size

$$N = \log ( n / \delta ) 2 ^ { p o l y l o g ( 1 / \epsilon ) },$$

where x t is sampled from an unknown distribution D and | y t -Tr ( Oρ x ( t )) | ≤ ϵ for any observable O with eigenvalues between -1 and 1 that can be written as a sum of geometrically local observables. With a proper choice of the efficiently computable hyperparameters δ 1 , δ 2 , and B , the learned function h ∗ ( x ) = w ∗ · ϕ z ( ) satisfies

$$\mathbb { E } _ { x \sim D } \left | h ^ { * } ( z ) - T r ( O \rho ( z ) ) \right | ^ { 2 } \leq \epsilon$$

with probability at least 1 -δ . The training and prediction time of the classical ML model are bounded by O ( nN ) = n log( n/δ )2 polylog (1 /ϵ ) .

The request in the first one, that the Hamiltonian has a constant gap, can be soften into that the correlation length has a finite upper bound.

Now we present the formal version of Thm. 2 in the main text.

Theorem 6 (Provably Efficient Adiabatic Learning (PEAL) Theorem) . Given n, δ, η &gt; 0 , 1 e &gt; ϵ &gt; 0 and a training data set { ( g α , P , Q ⃗ ⃗ ) l , ⟨ ˆ O α,i ⟩ } l N l =1 of size

$$N = \log ( n / \delta ) 2 ^ { p o l y l o g ( 1 / \epsilon ) },$$

where l ∈ { 1 2 , , . . . , N } is the index of data points, ( g α l ) is the coupling constant used when collecting the data, ( ⃗ , Q P ⃗ ) l is classical variables sampled from exact simulation data, ˆ O α,i is an observable with eigenvalues in [ -1 1] , that can be written as a sum of geometrically local observables, ⟨ ˆ O α,i ⟩ l is the ground state expectation value with an n -qubit gapped geometrically local Hamiltonian ˆ (( H ⃗ , Q P ⃗ ) ) l . Apply an ML predicted model M learned with a proper choice of the efficiently computable hyperparameters. When the Error Bounded Property in Def. 1 is satisfied, a T -step PEAL prediction Q i, PEAL ( T ) comparing to the exact dynamical process Q i, Exact ( T ) , i.e. the accumulated error, has an error bound

$$| Q _ { i, P E A L } ( T ) - Q _ { i, E x a c t } ( T ) | \leq O ( \sqrt { T \epsilon / \eta } )$$

with probability P ≥ 1 -η -δ .

If further assume the learning error of model M is sub-Gaussian distributed, the error bound can be improved into

$$| Q _ { i, P E A L } ( T ) - Q _ { i, E x a c t } ( T ) | \leq O ( \sqrt { 2 \epsilon \log ( 2 T / \eta ) } )$$

with probability P ≥ 1 -η -δ .

If further more assume the learning error of model M is bounded by √ ϵ almost surely, the error bound can be further improved into

$$| Q _ { i, P E A L } ( T ) - Q _ { i, E x a c t } ( T ) | \leq O ( \sqrt { \epsilon } )$$

with probability P ≥ 1 -δ , i.e. the error is bounded by a constant.

The accumulated error of classical variables P i and all k -local, bounded quantum observables also have the same scaling as the accumulated error of Q i .

Proof. To be able to apply the Error Bounded Property in Def. 1, we request that the prediction errors | δ O ⟨ ˆ ⟩| = |⟨ ˆ O ⟩ PEAL ( ⃗ P PEAL , Q ⃗ PEAL ) -⟨ ˆ O ⟩ Exact ( ⃗ P PEAL , Q ⃗ PEAL ) | are upper bounded for all T number of time steps. From the Thm. 5 above, we learn

$$\mathbb { E } | \delta \langle \hat { O } \rangle | ^ { 2 } \leq \epsilon,$$

from which we can estimate the probability of a single-shot prediction error going beyond a threshold:

$$\text{ite the probability of a single-shot prediction error going beyond a threshold: } \\ \mathbb { P } \left [ | \delta \langle \hat { O } \rangle | \geq A \sqrt { \epsilon } \right ] & = \mathbb { P } \left [ | \delta \langle \hat { O } \rangle | ^ { 2 } \geq A ^ { 2 } \epsilon \right ] \\ & \leq \frac { \mathbb { E } | \delta \langle \hat { O } \rangle | ^ { 2 } } { A ^ { 2 } \epsilon } \\ & \leq \frac { 1 } { A ^ { 2 } }, & & ( S 9 6 ) \\ \text{nt factor we choose to describe the threshold} \, \text{ in the second line } we anniv the Markov's$$

where A ≥ 0 is a constant factor we choose to describe the threshold. In the second line, we apply the Markov's inequality because | δ O ⟨ ˆ ⟩| 2 is non-negative.

For a dynamical simulation process with T number of prediction stpes, the probability of any prediction error going beyond the threshold is controlled by the union bound:

$$\mathbb { P } \left [ \bigcup _ { i = 1 } ^ { T } \left \{ | \delta \langle \hat { O } \rangle ( t _ { i } ) | \geq A \sqrt { \epsilon } \right \} \right ] \leq \sum _ { i = 1 } ^ { T } \mathbb { P } \left [ | \delta \langle \hat { O } \rangle ( t _ { i } ) | \geq A \sqrt { \epsilon } \right ] \leq \frac { T } { A ^ { 2 } }.$$

Taking into account that the trained model could have at most δ probability of failure, we get the probability of not being able to apply our PEAL method is at most ( T A 2 + ) by another union bound. δ

When we are able to control all the prediction errors under the threshold A √ ϵ , we have

$$| q ( t ) | = | Q _ { i, \text{PEAL} } ( T ) - Q _ { i, \text{Exact} } ( T ) | \leq O ( A \sqrt { \epsilon } )$$

by the Error Bounded Property in Def. 1. From the analysis above, we see this is applicable with a probability at least (1 -T A 2 -δ ).

For given η &gt; 0, we choose the threshold factor A = √ T/η , such that η = T A 2 . Therefore, we get our most general bound

$$| Q _ { i, \text{PELAL} } ( T ) - Q _ { i, \text{Exact} } ( T ) | \leq O ( \sqrt { T \epsilon / \eta } )$$

with probability P ≥ 1 -η -δ .

This most general bound is showing a √ T behavior, which happens to be the same as a diffusion model.

This (1 -η -δ ) behavior in the probability is expected, because we have an adiabatic hybrid algorithm with both quantum and classical updates. The δ failure comes from approximating the quantum algorithm, and the η failure comes from the adiabatic quantum-classical dynamics.

We can improve this bound if we can acquire more knowledge on the distribution of the prediction error δ O ⟨ ˆ ⟩ . For example, if δ O ⟨ ˆ ⟩ has a sub-Gaussian distribution with a moment-generating function (MGF)

$$M _ { \delta \langle \hat { O } \rangle } ( s ) = \mathbb { E } \left [ e ^ { s \delta \langle \hat { O } \rangle } \right ] \leq \exp \frac { \epsilon s ^ { 2 } } { 2 },$$

then we can apply the Chernoff bound and get

$$\ddot { \kappa } _ { \kappa } & \quad \lambda \\ \mathbb { P } \left [ | \delta ( \hat { O } ) | \geq B \sqrt { \epsilon } \right ] \leq 2 \exp \left ( - \frac { ( B \sqrt { \epsilon } ) ^ { 2 } } { 2 \epsilon } \right ) = 2 \exp \left ( - B ^ { 2 } / 2 \right ),$$

where B is also a threshold factor.

Repeat the union bound argument, we can get | Q i, PEAL ( T ) -Q i, Exact ( T ) | ≤ O B ( √ ϵ ) with a probability at least (1 -2 T exp ( -B / 2 2 ) -δ ). For given η &gt; 0, we choose B = √ 2 log(2 T/η ), such that η = 2 T exp ( -B / 2 2 . ) Therefore, we get our bound with assuming sub-Gaussian prediction error

$$| Q _ { i, \text{PEAL} } ( T ) - Q _ { i, \text{Exact} } ( T ) | \leq O ( \sqrt { 2 \epsilon \log ( 2 T / \eta ) } )$$

with probability P ≥ 1 -η - √ δ . This improves the general T behavior to a better √ log( T ) behavior.

We can improve this bound even further if we add even stronger assumption on the distribution of the prediction error δ O ⟨ ˆ ⟩ . If we assume there exists a constant factor C such that | δ O ⟨ ˆ ⟩| ≤ C √ ϵ almost surely, then we can directly apply our Thm. 3 and get

$$| Q _ { i, \text{PAL} } ( T ) - Q _ { i, \text{Exact} } ( T ) | \leq O ( \sqrt { \epsilon } )$$

with probability P ≥ 1 -δ .

This is a constant bound which is independent on T . i.e. The accumulated position error of the dynamical simulation is bounded by a constant.

Because we also have a controlled | p t ( ) | in Thm. 3, all the arguments above are also true if replacing Q by P . With the proof of the following lemma, we complete the proof of Thm. 6.

Lemma 1. All k-local, bounded observables have the same provably efficient bound under PEAL, no matter whether they are involved in the dynamics or not. For such an observable ˆ Ω , given a training data set { ( ⃗ , Q P ⃗ ) l , ⟨ ˆ Ω ⟩ } l N l =1 of size N same as in Thm. 6, under various assumptions in Thm. 6, the accumulated error |⟨ ˆ Ω ⟩ PEAL ( ⃗ P PEAL ( T , ⃗ ) Q PEAL ( T )) -⟨ ˆ Ω ⟩ Exact ( ⃗ P Exact ( T , ⃗ ) Q Exact ( T )) | holds similar bounds as in Eq. S92, S93, and S94.

Proof. The accumulated error for an observable is defined as

$$| \langle \hat { \Omega } \rangle _ { P E A l } ( \vec { P } _ { P E A l } ( T ), \vec { Q } _ { P E A l } ( T ) ) = \langle \hat { \Omega } \rangle _ { E x a c t } ( \vec { P } _ { E x a c t } ( T ), \vec { Q } _ { E x a c t } ( T ) ) |.$$

There are two source of error in this operator learning: one from the PEAL path ( ⃗ P PEAL ( T , ⃗ ) Q PEAL ( T )) deviate from the exact simulation ( ⃗ P Exact ( T , ⃗ ) Q Exact ( T )), the other from the single-step ML prediction error. In terms of formula, we have

$$& | ( \hat { \Omega } ) _ { \text{PEAL} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) - \langle \hat { \Omega } \rangle _ { \text{Assact} } ( \vec { P } _ { \text{Assact} } ( T ), \vec { Q } _ { \text{Assact} } ( T ) ) | \\ & \leq | \langle \hat { \Omega } \rangle _ { \text{PEAL} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) - \langle \hat { \Omega } \rangle _ { \text{Assact} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) | \\ & \quad + | \langle \hat { \Omega } \rangle _ { \text{Assact} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) - \langle \hat { \Omega } \rangle _ { \text{Assact} } ( \vec { P } _ { \text{Assact} } ( T ), \vec { Q } _ { \text{Assact} } ( T ) ) |.$$

$$\leq & \, | \langle \hat { \Omega } \rangle _ { \text{PEAL} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) - \langle \hat { \Omega } \rangle _ { \text{Exact} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) | \\ & \, + | \langle \hat { \Omega } \rangle _ { \text{Exact} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) - \langle \hat { \Omega } \rangle _ { \text{Exact} } ( \vec { P } _ { \text{Exact} } ( T ), \vec { Q } _ { \text{Exact} } ( T ) ) |. \\. & \,. & \,. & \,. & \,. & \,. & \,. & \,. & \,. & \,. & \,.$$

We first analyze the second term by the perturbation theory. We have

$$\langle \hat { \Omega } \rangle _ { \text{Exact} } ( \vec { P } _ { \text{PEAL} } ( T ), \vec { Q } _ { \text{PEAL} } ( T ) ) = \langle \text{GS} _ { \text{PEAL} } | \hat { \Omega } | \text{GS} _ { \text{PEAL} } \rangle, & & ( \text{S106} ) \\ \int _ { \hat { \bigwedge } } \vec { n } \quad \vec { n } \quad \vec { \tau } \quad \vec { \tau } \quad \vec { \tau } \quad \cdots \quad \cdots$$

$$\langle \hat { \Omega } \rangle _ { \text{Exact} }$$

where | GS PEAL/Exact ⟩ = GS( | ⃗ P PEAL/Exact ( T , ⃗ ) Q PEAL/Exact ( T )) ⟩ is the ground state wavefunction for the Hamiltonian with classical parameters ( ⃗ P PEAL/Exact ( T , ⃗ ) Q PEAL/Exact ( T )), respectively.

Note that ( ⃗ , Q P ⃗ ) are served as parameters of the Hamiltonian. For slightly different ( ⃗ P PEAL ( T , ⃗ ) Q PEAL ( T )) and ( ⃗ P Exact ( T , ⃗ ) Q Exact ( T )), the corresponding Hamiltonians are slightly different, and thus the difference between the ground states can be analyzed by first order perturbation

$$| G S _ { P E A L } \rangle - | G S _ { E x a$$

where | n Exact ⟩ is the n -th excited state for the Hamiltonian with classical parameters ( ⃗ P Exact ( T , ⃗ ) Q Exact ( T )). E GS,Exact and E n, Exact are energies of | GS PEAL ⟩ and | n Exact ⟩ , respectively. We define the accumulated errors of classical parameters ⃗ p = ⃗ P PEAL ( T ) -⃗ P Exact ( T ), ⃗ q = ⃗ Q PEAL ( T ) -⃗ Q Exact ( T ). ˆ V is the first order Taylor expansion of the Hamiltonian

$$\hat { V } = \sum _ { i } \left ( \frac { \partial \hat$$

Note the the dependence on ( ⃗, ⃗ p q ) in Eq. S108 only appears in ˆ , V which is linear in ( ⃗, ⃗ p q ). Moreover, Eq. S108 is bounded by looking at its left hand side. Therefore, we see that overall | GS PEAL ⟩ - | GS Exact ⟩ is of order O ⃗ p, ⃗ ( q ), and is controlled by the various bounds in Thm. 6.

Therefore, the second term can be controlled by

$$\begin{array} {$$

where ( GS | PEAL ⟩ - | GS Exact ⟩ ) and ( GS ⟨ PEAL | - ⟨ GS Exact | ) is of order O ⃗ p, ⃗ ( q ) and everything else is O (1). Therefore the last line in Eq. S110 is of order O ⃗ p, ⃗ ( q ) because both terms are of order O ⃗ p, ⃗ ( q ). Hence, the second error term satisfies the same bounds as Q does in Thm. 6, which in the proof we analyze the probability of all T steps' predictions having errors under a certain threshold.

The first term is the single-step prediction error δ ⟨ ˆ Ω ⟩ at T -th step, which is analyzed in the proof of Thm. 6. Combining the first term and the second term is equivalent to requiring ( T +1) steps' prediction having errors under a certain threshold, which can be analyzed with the same technique we used in the proof of Thm. 6. More precisely, we have ( T +1) terms in the union bound in Eq. S97. This is equivalent to having results in Thm. 6 by changing T into ( T +1). However, such change does not affect the scaling behaviors of the bounds with T .

Therefore, the first and second terms combined also satisfy the same bounds as Q does in Thm. 6. This completes the proof of Lemma 1, that the accumulated error of any k -local observable satisfy the same bounds in Thm. 6. With this lemma proved, we complete the entire proof of Thm. 6.

## VII. PEAL Implementation

Here we provide more details about how we implement the PEAL algorithm in our numerical experiment on the Holstein model.

PEAL Step I: Collecting data. For each of the six training values g ∈ { 1 3 . , 1 32 1 34 1 36 . , . , . , 1 38 . , 1 4 . } , we randomly sampled 18 independent initial conditions ⃗ Q t ( = 0) iid ∼ N (0 , Q v ) ⊗ L with the standard deviation Q v = 0 2. . L = 50 is the system size. We set ⃗ P t ( = 0) = 0. For each initial condition, We iteratively perform quantum steps and classical steps, for 10,000 epochs, to generate a dynamical evolution path. We set the simulation time step dt = 0 01 and . the total time is 100. In the quantum steps, we use exact diagonalization (ED) to obtain the quantum observables. In the classical steps, we use fourth order Runge-Kutta method. For each path, we randomly sampled 500 pairs of ( ⃗ Q shifti ( ) t , n i ( )) data, where t t is a random time step and i is a random site. ⃗ Q shifti is the ⃗ Q vector with elements rolled by a shift of i , so that the i -th element is put on the first place. Overall, we have 6 × 18 × 500 = 54 000 data , pairs, which effectively serve as 1,080 samples on the size L = 50 system.

The purpose to use this kind of site-shifting data pairs is to respect the translation symmetry of the system. With periodic boundary condition, the system has translation symmetry. Therefore, we can use the same model to predict n i on any location in the system, as long as we also translate the ⃗ Q vector accordingly.

PEAL Step II: ML training. We use the provably efficient ML algorithm introduced in [27]. For each length-50 vector ⃗ Q , we generate 50 local regions, each local region having radius 1 and consisting of 3 elements. We perform a random Fourier feature map on the local regions, generate R number of cosine features and R number of sine features for each local region. R is a hyperparameter chosen from { 5 10 20 40 80 160 , , , , , } . The frequency factor γ ω for the random Fourier feature map is another hyperparameter chosen from { 0 3 . , 0 6 . , 1 2 3 6 10 20 , , , , , } . We trained an l 1 -regularized regression (LASSO) on the non-linear features. The strength of regularization α is determined by

LassoLars method with 4-fold cross validation [58]. For 1,024 samples, the grid-search best hyperparameters are R = 20, γ ω = 6, and α = 4 98 . × 10 -6 . The best model has 316 non-zero linear coefficients after the LASSO feature selection.

PEAL Step III: ML prediction. To speed up the dynamical simulation, we replace the ED steps by the prediction from the learning model. We take care of the U (1) global symmetry of the system. After each ML boosted quantum step, we subtract the electron density n i by the mean of their excess over 1 / 2. Therefore, we can make the total electron number conserved at half filling, and respect the U (1) global symmetry of the system. See more details in the following section.

## VIII. Details of Symmetry-preserving PEAL

In the application of PEAL to specific systems with symmetries, we would like to preserve the symmetries during our numerical simulation. Here we present symmetry-preserving PEAL for two types of symmetries: The U (1) global symmetry and the translation symmetry.

To preserve the U (1) global symmetry, we would like to have the total charge conservation at each time step during the dynamics. When we sum up the model's prediction of the electron charge density at each site of the system, often case the result is not the same as the total charge in the previous time step. We can make a correction on the electron charge density at each site during each time step of the dynamics, in order to conserve the total charge.

There are two possible ways to do the correction. The first one is to uniformly add or subtract a constant on the electron charge density at each site. The second is to uniformly multiply a factor to the electron charge density at each site. However, the second way does not treat an electron (charge density equals to 1) and a hole (charge density equals to 0) equivalently. Therefore, we use the first way to do the correction.

Suppose the learning model M predicts the electron charge density n i = ⟨ c c † i i ⟩ on each site i . In order to preserve the U (1) global symmetry, we would like to have the corrected electron charge density ˜ n i such that ∑ i ˜ n i = L/ 2, assuming at half filling for a system with L sites. It is straightforward to see that we can achieve this charge conservation by doing the correction n i → ˜ n i = n i +∆ n with ∆ n = 1 L ∑ i ( 1 2 -n i ).

In the proof of PEAL error bounds, we analyzed all situations by turning the problem into the calculation of the probability of | δ O ⟨ ˆ ⟩| ≤ D √ ϵ , where the constant D stands for A B , , or C in the proof of Thm. 6. Note that when all | δn i | ≤ D √ ϵ , we have

Now we prove the error bounds for our PEAL algorithm are still valid with this correction. Suppose each n i has a prediction error δn i = n i -n ∗ i , where n ∗ i is the electron charge density computed with QSS, which satisfies ∑ i n ∗ i = L/ 2. We can see ∆ n = -1 L ∑ i δn i , and the error of the corrected electron charge density is δn ˜ i = ˜ n i -n ∗ i = δn i +∆ = n δn i -1 L ∑ i δn i .

$$| \delta \tilde { n } _ { i } | = | \delta n _$$

which means that we are safe to apply every error bounds after we raise the corresponding constant D by a factor of 2. This completes the proof.

When the system has a translation symmetry, such as the one for the periodic chain we studied in the main text, every site i is on equal footing. Therefore, we will get the identical error distribution when we apply a model M i trained for n i to another site j to predict n j . This means that we can save our cost by applying one model repeatedly on every site i instead of training an independent model for every n i , and we have the error bounds unaffected. Moreover, by doing so, we preserve the translation symmetry in the PEAL prediction: If instead of the original initial condition ⃗ Q (0) = ( Q (0) 0 , Q (0) 1 , . . . , Q (0) L -2 , Q (0) L -1 ) we input a shifted initial condition ( Q (0) 1 , Q (0) 2 , . . . , Q (0) L -1 , Q (0) 0 ), then all outputs from the symmetry-preserving PEAL algorithm will be also shifted correspondingly comparing to the original outputs.

## IX. Details of Numerical Experiments

We are interested in the quantum dynamics of the Holstein model in Eq. 8 with a random initial distribution of Q i (0), and we set initial momentum P i (0) = 0. Under the adiabatic approximation, the quantum dynamics can be

further described by three equations:

$$Q _ { i } ^ { ^ { \prime } } ( t ) = \frac { 1 } { M } P _ { i } \quad \quad \.$$

$$P _ { i } ^ { ^ { \prime } } ( t ) & = - k Q _ { i } + g ( n _ { i } - \frac { 1 } { 2 } ) - \gamma P _ { i } \\ & \quad \ \ 1,$$

$$H _ { e } ( \{ Q _ { i } ( t ) \} ) = - t \sum _ { i, j } c _ { i } ^ { \dagger } c _ { j } - g \sum _ { i } ^ { \mu } ( c _ { i } ^ { \dagger } c _ { i } - \frac { 1 } { 2 } ) Q _ { i }$$

where n i = ⟨ c c † i i ⟩ and γ is the damping coefficient.

In the experiment, we observed three regimes of g values. For a small g value, the final stable configuration of the system has no CDW. This is due to the finite size effect. For the system size L = 50 we used, the critical value to generate stable CDW is g crit = 1 152. . This critical value can be analytically solved and we present the calculation already in earlier section. For a large g value, domain walls can be generated in the system. The characteristic length of the domain walls decreases when g increases. Therefore, when the system size is much larger than the domain wall characteristic length, domain walls can be relatively far from each other, which makes their interaction exponentially small and the equilibrium time exponentially large. For the system size L = 50 we used, the value to generate longlasting domain walls is g DW ∼ 1 6. . For the value g crit &lt; g &lt; g DW , we can see a clear CDW developing during the dynamics. Therefore, we collect a data set for g ∈ { 1 3 . , 1 32 . , 1 34 1 36 1 38 1 4 . , . , . , . } , mimicking a uniform distribution in [1 3 . , 1 4]. .

For each g value, we generate 18 ED dynamical evolution paths with random initial conditions, serving as the training data. We also generate 10 more ED dynamical evolution paths with random unseen initial conditions for each g ∈ { 1 3 . , 1 32 1 34 . , . , 1 36 . , 1 38 1 4 . , . } and for each g ∈ { 1 31 . , 1 33 . , 1 35 1 37 1 39 . , . , . } , serving as the standard learning test set and the transfer learning test set, respectively. In Figure 3 in the main text, we show how the test error scales with the number of samples used in model training. We trained models with number of samples from 1, 2, 4, ..., to 1,024. (Each sample contains 50 data pairs based on the system size we used.) We see an error scaling law agrees with the predicted logarithmic scaling. The transfer learning test error is only slightly higher than the standard learning test error. The scatter plots of the target values show the model has learned nicely for both the normal case and the transfer learning case.
Changnan Peng Jin-Peng Liu Gia-Wei Chern Di Luo
Provably Efficient Adiabatic Learning for Quantum-Classical Dynamics
Quantum-classical hybrid dynamics is crucial for accurately simulating
complex systems where both quantum and classical behaviors need to be
considered. However, coupling between classical and quantum degrees of freedom
and the exponential growth of the Hilbert space present significant challenges.
Current machine learning approaches for predicting such dynamics, while
promising, remain unknown in their error bounds, sample complexity, and
generalizability. In this work, we establish a generic theoretical framework
for analyzing quantum-classical adiabatic dynamics with learning algorithms.
Based on quantum information theory, we develop a provably efficient adiabatic
learning (PEAL) algorithm with logarithmic system size sampling complexity and
favorable time scaling properties. We benchmark PEAL on the Holstein model, and
demonstrate its accuracy in predicting single-path dynamics and ensemble
dynamics observables as well as transfer learning over a family of
Hamiltonians. Our framework and algorithm open up new avenues for reliable and
efficient learning of quantum-classical dynamics.