Adaptive Variational Quantum Dynamics Simulations

We propose a general-purpose, self-adaptive approach to construct variational wavefunction ans\"atze for highly accurate quantum dynamics simulations based on McLachlan's variational principle. The key idea is to dynamically expand the variational ansatz along the time-evolution path such that the ``McLachlan distance'', which is a measure of the simulation accuracy, remains below a set threshold. We apply this adaptive variational quantum dynamics simulation (AVQDS) approach to the integrable Lieb-Schultz-Mattis spin chain and the nonintegrable mixed-field Ising model, where it captures both finite-rate and sudden post-quench dynamics with high fidelity. The AVQDS quantum circuits that prepare the time-evolved state are much shallower than those obtained from first-order Trotterization and contain up to two orders of magnitude fewer CNOT gate operations. We envision that a wide range of dynamical simulations of quantum many-body systems on near-term quantum computing devices will be made possible through the AVQDS framework.


I. INTRODUCTION
One of the primary scientific quantum computing focuses has been to simulate ground-state, excited-states, and dynamical properties of spin and fermion systems [1][2][3][4][5][6][7][8][9][10][11][12][13]. The ultimate goal is to accurately model physical systems of classically prohibitive size by efficient encoding and manipulation of many-body states. For ideal fault-tolerant quantum computers where deep circuits can be executed, algorithms built on Trotterized adiabatic state preparation and dynamics simulations, along with quantum phase estimation, can address a wide class of physical and chemical problems in and out of equilibrium [3,4,7,10,11,13].
In the near term with noisy intermediate-scale quantum (NISQ) computers [14], practical calculations on quantum devices are limited to short circuits. To exploit the emerging NISQ technology, the variational quantum eigensolver (VQE) has been developed and demonstrated to prepare the ground state of a time-independent Hamiltonian [15][16][17][18][19][20], or generally to minimize a static cost function with individual terms whose expectation values can be measured efficiently on quantum devices [21]. Preconstructed fixed-form variational ansätze, such as the unitary coupled-cluster ansatz [16,17,22], have been employed in early VQE calculations of molecules on quantum devices. However, the accuracy is often limited by the form of the ansatz [23], and the number of variational parameters and the circuit depth can grow as a high-order polynomial as the number of orbitals and atomic sites increases [16,24]. To alleviate these issues arising from a fixed form of the variational wavefunction, several adaptive VQE methods have been proposed. These methods use a form of the variational wavefunction that is adaptively optimized for the specific problem [23,25], leading to highly accurate results with much simpler variational * ykent@iastate.edu circuits.
Variational approaches to quantum dynamics simulations (VQDS), including fast-forwarding methods, have also been proposed and applied to quantum spin models [26][27][28][29][30], with proof-of-principle applications on real devices [31]. The proposed variational wavefunction forms are described by a set of relatively shallow quantum circuits, which are tailored for execution on NISQ devices. The quality of variational quantum simulations is tied to the ability of the variational ansatz in describing the time-evolved wavefunction. For VQDS, the ansatz should be flexible enough to represent the quantum state along its time-evolution path. This is typically much more challenging than constructing an accurate variational ansatz for the ground state, because the nature of the wavefunction can change significantly during time evolution. Attempts to construct ansätze of fixed variational circuits for VQDS have been reported [31], but the simulation accuracy quickly deteriorates as the system grows from 2 sites to a few sites. This highlights the need for flexible variational circuits that can adapt to the changes of the wavefunction during time-evolution, while still keeping the circuit sufficiently shallow to be run on NISQ quantum processing units (QPUs).
Here, we propose a novel time-dependent adaptive variational method to perform accurate quantum dynamics simulations of fermionic and quantum spin models. Going beyond a variational approach with a fixed cost function, the proposed scheme constructs an efficient time-dependent variational ansatz of the time-evolved wavefunction. This directly addresses the challenge of constructing efficient variational ansätze for dynamics simulations. This is generally a difficult task as the time-evolved wavefunction explores different, and a priori unknown, regions of Hilbert space. The proposed adaptive variational quantum dynamics simulation (AVQDS) approach is built on the McLachlan variational principle for real-time quantum dynamics simulations [26,32], and automatically generates and dynamically expands the variational ansatz by minimizing the McLachlan distance L 2 along the time-evolution path. The form of the variational ansatz is incrementally expanded by choosing optimal operators from a predefined operator pool to construct additional unitary gates, in the same spirit as in the (static) ADAPT-VQE method [23,25]. The crucial difference is that here we optimize a time-dependent cost function, the McLachlan's distance L 2 (t) at time t. This allows the complexity of the ansatz to increase as needed during the time-evolution in order to accurately represent the wavefunction. It is worth noting that our key idea of adaptively generating variational ansätze by minimizing a time-dependent cost function for dynamics simulations is widely applicable beyond the McLachlan approach. For example, it can also be used to generalize the projected-Variational Quantum Dynamics (p-VQD) simulation method discussed in Ref. [33] by replacing the McLachlan distance with the step-infidelity function.
We apply AVQDS to study linear ramp quantum dynamics of the integrable Lieb-Schultz-Mattis (LSM) spin model [34,35], and sudden quench dynamics of the nonintegrable mixed-field Ising model (MFIM). In both cases, we find that dynamical quantities of interest (such as local observables, total energy and the Loschmidt echo) are described accurately with the adaptively generated variational ansatz. Notably, the state preparation circuits for the time-dependent wavefunction require two orders of magnitude fewer two-qubit gates than first-order Trotter dynamics simulations that achieve comparable accuracy. We demonstrate an initial time t-linear growth of the number of CNOT gates N cx ∝ t in the AVQDS circuits, before N cx saturates after a critical time t s . The saturation time t s is found to increase with system size N . The system size scaling of the saturated N cx (t > t s ) changes from quadratic (∝ N 2 ) to higher order (∝ N α with α > 4) as the integrability of the model is broken. In contrast, at fixed simulation times t < t s , N cx scales approximately linearly with large N > N c (t) (where the circuits are in the presaturation regime). The crossover system-size N c grows slowly with t. This implies the practical scalability of general AVQDS simulations over finite time intervals.
The paper is organized as follows. The AVQDS algorithm, calculation procedures and important technical details are presented in Sec. II. For completeness, the section also contains brief reviews of the Trotterized dynamics and the previously introduced VQDS formalism. Sec. III discusses results for the integrable LSM model using a linear-ramp quench protocol and the observed gate count scaling of the AVQDS method. AVQDS simulations of quench dynamics in the nonintegrable MFIM are presented in Sec. IV. We summarize our results and give concluding remarks in Sec. V.

A. Trotterized state evolution
For the convenience of later comparisons, we summarize the quantum dynamics simulation method using a discretetime propagator based on the Trotter decomposition [36,37]. For a system described by a generic time-dependent HamiltonianĤ the quantum state evolves by a fixed time step δt as Here, e −iδtĥµ[t] is a unitary that can be efficiently implemented on QPUs. Therefore, the Trotter circuit depth grows linearly with the number of time steps, with the incremental depth determined by the Hamiltonian terms {ĥ µ }. For calculations on quantum computers, each individual termĥ µ is a tensor product of Pauli operators. The exponential of a Pauli term coupling p qubits scaled by an imaginary number constitutes a general Pauli rotation gate, which can be compiled to a series of one-qubit rotations and 2(p − 1) two-qubit entangling gates for real quantum device applications [37].
Although it is not feasible to directly compile a general discrete-time Hamiltonian evolution operator in real quantum devices, on classical devices it can be more efficient to calculate the discrete-time evolution according to . This approach is more tolerant to finite step sizes, and becomes exact for simulations of sudden quench dynamics where the system is evolved by a Hamiltonian that is piecewise constant in time. We adopt this approach to get numerically exact data for comparison with the simulation results to be discussed later.

B. Variational Quantum Dynamics Simulation
The formalism of variational quantum dynamics simulations (VQDS) has been systematically presented in Ref. [26]. To facilitate later discussions, we summarize the main points of variational real-time dynamics simulations based on McLachlan's variational principle [32]. We work at the level of the density matrix ρ, which eliminates the necessity of keeping track of the global phase of the state |Ψ . For a system in a pure state evolving under a timedependent HamiltonianĤ, the density matrix ρ = |Ψ Ψ| evolves according to the von-Neumann equation . The circuits to measure matrix M and vector V can be found in reference [26] (see also Fig. 2). Note that in the ansatz adaptive procedure, one only needs to measures the incremental elements in M and V , which are added in a given step.
Here ρ F ≡ Tr[ρ † ρ] is the Frobenius norm of the matrix ρ. The matrix M is real symmetric and defined as which is equivalent to the quantum Fisher information matrix associated with the state fidelity [38]. The vector V is given by where Ĥ θ ≡ Ψ[θ]|Ĥ |Ψ[θ] , and which describes the energy variance ofĤ in the variational state |Ψ[θ] . The second term in the bracket of Eqs. (5) and (5) originates from the global phase contribution [26].
In the geometric picture of quantum evolution [39], the integral of energy variance with respect to time is independent of the specific Hamiltonian and defines a distance along evolution path measured by the Fubini-Study metric [40]. The minimization of the cost function Eq. (4) with respect to {θ µ } leads to the following equation of motion for the variational parameters: The McLachlan distance L 2 of the variational ansatz Ψ[θ] at optimal parameter values can then be calculated as As pointed out in Ref. [26], L 2 provides a natural measure for the accuracy of quantum dynamics simulations.

Algorithm and flow chart
The adaptive variational quantum dynamics simulation method is illustrated in Fig. 1. The technique dynamically constructs a variational ansatz of the following pseudo-Trotter form: where |Ψ 0 is a reference state and {θ µ } (µ = 0, . . . , N θ − 1) are the time-dependent variational parameters.Â µ is a Hermitian operator. The set of N θ operators {Â µ } will be dynamically expanded by including additional operators from an operator pool to maintain the McLachlan distance L 2 (8) below a threshold L 2 cut , if necessary, as the system evolves. In practice, we find L 2 cut = 10 −3 is enough to get highly accurate results.
Without loss of generality, let the dynamics simulation start at t = 0 with the system in a pure state ρ 0 = |Ψ 0 Ψ 0 |, which we choose as the reference state of the variational ansatz. More specifically, at t = 0 the variational ansatz |Ψ[θ] = |Ψ 0 has no variational parameters (i.e., N θ = 0). After an incremental time step t = δt, the Hamiltonian becomesĤ t . The McLachlan distance L 2 is measured with the previously obtained ansatz state |Ψ[θ] through the MMD module, which is specified in Fig. 1(b) according to Eqs. (5)- (6). In the initial case where the variational ansatz has no parameters, the McLachlan distance is determined by the energy variance ofĤ t only. The energy variance is generally larger than zero, because the ansatz state is not an eigenstate ofĤ t as the system evolves in time. The McLachlan distance L 2 is then compared against the threshold L 2 cut , and the ansatz adaptive procedure is triggered if L 2 ≥ L 2 cut . The adaptive procedure starts with evaluating the McLachlan distance L 2 with respect to a series of new variational ansätze. Each new ansatz, e −iθ Â µ | θ =0 |Ψ[θ] , is composed of the existing ansatz, multiplied by the exponential of an generatorÂ µ with a coefficient θ initialized to zero at the current time step. Although the new ansatz with θ = 0 does not change the ansatz state, the McLachlan distance L 2 can change due to nonvanishing derivatives with respect to the additional parameter in Eqs. (5) and (5). Here, the choice ofÂ µ runs through all the operators in a preconstructed (fixed) operator pool of size N p .
For each operator e −iθ Â µ that is added to the ansatz, the dimension of θ increases from N θ to N θ + 1. Accordingly, the dimension of the symmetric matrix M (5) increases from N θ × N θ to (N θ + 1) × (N θ + 1), and the dimension of the vector V (5) increases from N θ to N θ +1. Because the additional parameter in the new ansatz is always initialized to zero, the represented ansatz state remains the same. Therefore, only the additional (N θ + 1) th row of the matrix M and the final element of V need to be evaluated. The obtained {L 2 µ } (µ = 0, . . . , N p − 1) values are compared, and one selects the operatorÂ ν for which L 2 ν is minimal. The ansatz form is updated to |Ψ[θ] → e −iθνÂν |Ψ[θ] , with θ ν initially set to zero and the number of variational parameters N θ increased by one. Note that setting θ N θ +1 = 0 initially ensures that the wavefunction is smooth during time evolution. The McLachlan distance L 2 is then updated and compared with the threshold L 2 cut . For L 2 ≥ L 2 cut the ansatz adaptive procedure is repeated until L 2 < L 2 cut is satisfied. Only then the variational parameters are updated, θ → θ + δθ, at current time step according to In typical AVQDS simulations reported here, the number of unitaries added to the ansatz at some initial time steps can be as much as about 10 per step to gain sufficient expressibility, but quickly decreases to about 1 per step if the McLachlan distance L 2 goes beyond the threshold L 2 cut . Finally, the system evolves to the next time step and the algorithmic procedure continues until the simulation time period ends.

Initial state preparation
The AVQDS calculation starts with some initial state, which may often be the ground state |Ψ 0 of a Hamilto-nianĤ 0 at t = 0. A general efficient state preparation algorithm for near-term quantum devices remains an open challenge and has attracted much research effort, including adiabatic state preparation [6], VQE [16], and quantum imaginary time evolution [41][42][43]. Naturally, one can replace the original proposal of Trotter dynamics-based adiabatic state preparation with the AVQDS approach. This corresponds to first performing an AVQDS simulation starting from a Hamiltonian with a tensor-product ground state and evolving to the HamiltonianĤ 0 with ground state |Ψ 0 at a sufficiently slow quench rate. Alternatively, the ground state of the initial Hamiltonian H 0 can be obtained using (static) ADAPT-VQE [23,25], which can easily be combined with AVQDS. More specifically, we use the recently proposed qubit-ADAPT-VQE technique [25] to prepare the ground state Ψ 0 of the initial HamiltonianĤ 0 . Compared with AVQDS, which dynamically updates the variational ansatz (9) with the goal of minimizing McLachlan's distance along a time-evolution path, qubit-ADAPT-VQE uses an adaptive scheme to optimize the variational ansatz in Eq. (9) to minimize the expectation value of a time-independent Hamiltonian. Note that when encoding the initial state in this way, the initial number of variational parameters in AVQDS is larger than zero. The dynamical ansatz adaptive procedure can of course be carried out in the same way as described above.
In the AVQDS simulations reported below, we focus on the dynamics simulation part, where we dynamically construct the variational ansatz and update the parameters according to Eq. (10). Further details on state preparations using qubit-ADAPT-VQE and its combination with AVQDS will be given in Section III D.

Important technical details
Let us discuss some important technical details in the practical implementation of the AVQDS algorithm presented above, which can be used for dynamics simulations of generic Hamiltonian systems, including fermionic and spin systems. Since the system Hamiltonian will always be transformed to a qubit representation for calculations on QPUs, the operator pool {Â µ } can be constructed from a set of Pauli terms, i.e., tensor products of Pauli operators at different spin-orbital sites. The pool can be made of a complete list of Pauli terms up to some fixed length, or only involve Pauli terms that appear in the system HamiltonianĤ. Symmetries leading to conservation of particle number and spin quantum numbers, time-reversal symmetry and point group symmetries, can be considered to further reduce the pool size [47,48]. For fermionic systems, the operator pool can also be constructed using fermion operators, i.e., tensor products of fermion creation and annihilation operators, subject to symmetry constraints, before translating into qubit operators. This can potentially reduce the size of the operator pool for simulations of fermionic systems, at the cost that each operator maps, in general, to a sum of several Pauli terms. In the numerical calculations of spin models that we present below, we adopt the Hamiltonian pool which is composed of Pauli terms present in the Hamiltonian.
The AVQDS approach amounts to numerically integrat-ing a system of ordinary differential equations (7). Within the well-known Euler method, the local truncation error at a single time step is of order of (δt) 2 [49], where δt is the step size. This leads to a global truncation error over the total simulation period that scales linearly with δt. Higher-order methods such as the Runge-Kutta technique yield more favorable scaling [49]. In the numerical simulations presented here, we adopt the Euler method for simplicity and a fair comparison with the first order Trotter dynamics simulations. It is useful to contrast the change of the wavefunction that occurs during one timestep within Trotterized and AVQDS simulations. Within Trotter dynamics, the size of the timestep δt controls the change of the wavefunction during a single step (see Eq. (2)). In contrast, in AVQDS the change is controlled by the change of the variational parameters, δθ (see Eq. (10)). This leads to an effective way to stabilize AVQDS simulations by fixing a maximal step size δθ max . In other words, the time step size δt is dynamically adjusted at each time step such that max µ (|δθ µ |) ≤ δθ max . The dynamical time step size δt in AVQDS is therefore set by δθ max and the maximal absolute value |θ| max of elements in the vectoṙ θ in Eq. (10). When comparing AVQDS to Trotterized dynamics simulations, δθ max should be chosen to be the same as the fixed time step size δt in Trotter dynamics to achieve a similar accuracy. An additional consequence is that |θ| max determines the total number of time steps required in AVQDS calculations to reach a given final time.
We note that in the numerical simulations presented below, we find that the evaluation ofθ = M −1 V at intermediate time steps can involve a matrix M with large condition number [50], which are often encountered in numerical statistical analysis and machine learning [51]. The issue can be alleviated by adding a small diagonal element (ξ = 10 −6 ) to the matrix M before inversion following Tikhonov regularization, which effectively penalizes the large magnitude ofθ. As a result, we find that |θ| max falls in the range of (1, 10) in the following calculations. For AVQDS calculations in the presence of noise (both quantum mechanical shot noise and gate noise), we expect that the regularization parameter ξ must be increased in order to deal with the resulting fluctuations in M and V [41,52]. Alternatively, the matrix inversion in Eq. (10) can be avoided altogether by solving the equation of motion (7) using optimization techniques, which may also provide additional channels for gate error mitigation [53].
AVQDS simulations reported here are carried out with a classical implementation based on Quantum Toolbox in Python (QuTiP) [54]. Because AVQDS relies on the same set of measurements as regular VQDS with a fixed wavefunction ansatz, the quantum circuit implementation on NISQ QPUs follows that of VQDS discussed previously [26,44]. Specifically, Fig. 2 lists the unique terms in determining the symmetric matrix M in Eq. (5), vector V in Eq. (5) and scalars Ĥ θ , Ĥ 2 θ listed in the left column. The reduced expressions with the choice of pseudo-Trotter |0 for an N -qubit system. Two types of quantum circuits are adopted: green block for the direct measurement circuit, and blue block for generalized Hadamard test circuit [26,44] |∅ is given by 2P |0 − 1, with P |0 the probability that the ancillary qubit is in state |0 .ÛA[θ] can beÛ0,µ−1 orÛ0,ν−1[θ] as dark red color encoded in the blue box of middle column, with similar color encoding for the other unitaries.σB represents a Pauli string inÂµ orĤ, which is typically defined on one or two qubits for spin models. Therefore, the controlled-σ gate represents a controlled-single qubit or controlled two-qubit gate. Alternatively, the generalized Hadamard test circuit can be replaced with direct measurement circuit [45,46]. form (2) for the variational ansatz are listed in the middle column, which can be measured by direct measurement circuits (green) and generalized Hadamard test circuits (blue), as shown in the right column.
For an estimation of quantum resource of VQDS calculations, let us consider a Hamiltonian composed of N H Pauli strings and N θ variational parameters in the wavefunction ansatz (9). The upper bound of the number of distinct direct measurement circuits and generalized Hadamard test circuits is (N H + 2)N θ + N H + N 2 H and N H (N θ − 1) + (N θ )(N θ − 1)/2, respectively. The ansatz adaptive expansion procedure in AVQDS adds a marginal overhead of quantum resources. The additional terms (a, b, c) as highlighted in Fig. 2 are to be measured only for unitaries to be appended in the variational circuit at the same state of the current time step, as discussed previously. Specifically, with a Hamiltonian operator pool composed of Hamiltonian terms as adopted in the following calculations, no additional measurements are required for terms (b, c) and part of terms (a), because all the contributions have already been measured when evaluating the expectation values ofĤ andĤ 2 . Therefore the additional measurements in AVQDS from the ansatz adaptive procedure are for terms (a), Re ∂ Ψ[θ]| ∂θµ ∂|Ψ[θ] ∂θν , with ν running through the Hamiltonian operator pool and 0 ≤ µ < N θ − 1, which amount to N H (N θ − 1) generalized Hadamard test circuits. Each of the generalized Hadamard test circuits includes at most two controlled two-qubit gates for spin models presented here, which can also be replaced with direct measurement circuits [45,46]. For a system described by N qubits with N H ∝ N p , such as local spin models with N H ∝ N , and N θ ∝ N q , the leading order of distinct measurement circuits scales as N 2 max(p,q) , with the circuit depth tied to N θ . The overall measurement cost of AVQDS approach due to finite sampling shot noise follows the analysis of general metricaware variational quantum algorithms [52]. To keep the shot noise induced error below , the number of samples generally scales as O(1/ 2 ), where the prefactor can be reduced by optimal measurement distribution [52,55].

A. Model and phase diagram
To demonstrate and benchmark the performance of the AVQDS method, we perform a series of finite-rate quantum quench dynamics simulations of the integrable N -site spin-1 2 Lieb-Schultz-Mattis chain with open boundary conditions, which describes an anisotropic XY Hamiltonian in a transverse magnetic field [34] whereX,Ŷ andẐ are single-site Pauli operators. The coupling constant J is set to one as the energy unit. h z is the magnetic field strength. The anisotropy in the xy plane is controlled by the parameter γ, and rotational symmetry around the z-axis is obtained with γ = 0. The ground state phase diagram of the model in the thermodynamic limit (N → ∞) is well known [35], and shown in Fig 3. The phase diagram is composed of two ferromagnetic phases with magnetic moment in x-direction (FM x ) and y-direction (FM y ) and a paramagnetic phase (PM), along with multiple phase boundaries and tricritical points.

B. Quench protocol and operator pool
We consider the linear ramp protocol: γ(t) = 1 − 2t T with 0 ≤ t ≤ T , as shown in Fig. 3. The longitudinal mag-netic field is set to h z = −0.7 and 1.6 to avoid degenerate ground states that occur along the vertical line h z = 0. In the adiabatic limit, T → ∞, and in the thermodynamic limit, N → ∞, the system evolves from the FM x phase, crosses a phase boundary and enters the FM y phase. In the following, we choose a finite quench speed of T = 3, such that nontrivial spin dynamics is developed in the linear ramp process. System sizes N ∈ [2, 11] have been considered. We restrict the operator pool to Pauli terms that appear in the Hamiltonian (11). The quantum gates to be applied on the reference state |Ψ 0 , which are the exponentials of the scaled Pauli terms appearing in Eq. (9), include single-site and two-site Pauli rotation gates. Since the Hamiltonian (11) only contains nearest-neighbor coupling terms, an expanded operator pool, which includes a complete set of two-site Pauli terms, has also been investigated. Interestingly, we find that AVQDS with the expanded pool generally produces longer variational circuits, although the simulation results are of the same accuracy. The reason can be attributed to the fact that the list of new McLachlan distances {L 2 ν }, see Eq. (8) and Fig. 1, often contains almost degenerate minimal values among several operators. The "biased" choice of operators in the physical pool shows some advantage due to direct connections with the Hamiltonian, which governs the quantum dynamics of the system.

C. Simulation results
To characterize the time-evolution of the quantum state |Ψ(t) under a quench with inverse speed T = 3, we calculate the instantaneous total energy and spin correlation functions. Results are presented in Fig. 4(a-d) for the LSM model with N=8, where the system is further evolved for an additional time period of T after the linear ramp under the final HamiltonianĤ(T ), in order to assess the flexibility of the variational ansatz in describing the postquench dynamics. The upper (lower) panels show results for system size h z = −0.7(1.6). In Fig. 4(e,f), we show a comparison of the number of two-qubit gates in the quantum circuits that describe the Trotterized dynamics and AVQDS. The AVQDS results, shown as solid lines, are in excellent agreement with the exact data indicated by symbols over the full time range. At the end of linear ramp t = T , the ansatz fidelity f ≡ | Ψ[θ(T )]| Ψ exact (T ) | 2 is beyond 99.9%. Numerically exact results for reference are obtained by propagating the state by matrix exponentiation, |Ψ[t + δt] = e −iĤ[t]δt |Ψ[t] , for a fine time mesh with step size δt exact = 5 × 10 −4 . In the adiabatic limit, the results for the energy and spin correlation functions are symmetric under the transformation t → T − t for 0 ≤ t ≤ T , but the finite quench speed breaks that symmetry. In the initial state, the dominant spin correlations are S x S x , as exemplified by the correlations between the first and second site, S x 0 S x 1 , and between the first and the last site, S x 0 , S x N −1 . During time-evolution, these correlations decrease as the parameter γ decreases, which reduces the Compared with linear circuit growth of Trotter dynamics simulations, AVQDS circuits are much shallower, with a sublinear circuit growth rate that quickly decreases, as plotted in panels (e) and (f). The ansatz at t = T is composed of 50 two-qubit Pauli rotation gates (100 CNOTs) for both cases, which slightly increases to 53 two-qubit gates at t = 2T . A multiplying factor of two for CNOT gates is included, since each two-qubit Pauli rotation gate can be compiled to two CNOTs along with single-qubit gates assuming full connectivity [37]. Inset in (f): variation of Hamiltonian parameter γ as a function of t.
strength of the S x i S x i+1 interactions. Instead, the spin correlations among the y-components increases, as exemplified by S y 0 S y 1 . The system energy increases and reaches a maximum close to γ = 0, where the phase transition occurs in the thermodynamic limit. For γ < 0, the energy begins to decrease as the S y i S y j -spin correlations continue to grow. Duo to the nonadiabatic finite speed quench, the long-range correlation S y 0 , S y N −1 , which requires longer time to establish, remains far from equilibrium value at t = T . Although the dynamical spin correlations are remarkably distinct between the two h z points of parameter space, the adaptive variational circuits reach about the same complexity of 50 two-qubit Pauli rotation gates at the end of linear ramp. The post-ramp dynamics in the following time period of T slightly expands the ansatz by 3 two-qubit gates.
To compare the circuit complexity of AVQDS with the Trotter approach, we perform Trotterized simulations with a uniform time step δt = 5 × 10 −3 , applying a series of one and two-qubit Pauli rotation gates to the state. To characterize the accuracy of the dynamics simulations, we evaluate the standard deviation of an observableÔ along the linear ramp dynamical path for 0 ≤ t ≤ T according to where the summation is over the entire time mesh of dimension N t for Trotter simulations. Note that the expectation value Ô t at a specific time t in the simulation is calculated rigorously without any noise, which corresponds to infinite repeated measurements of the observableÔ on ideal fault-tolerant quantum computers. The standard deviation is 0.003 for both energy and the shown spin correlation functions for a system size of N = 8. As discussed previously in section II C 3, the maximum allowed parameter step size δθ max is the proper controlling factor in AVQDS, replacing δt that is the relevant quantity in Trotter dynamics simulations.
To properly benchmark AVQDS and obtain a reasonable comparison of the resulting circuit complexities, we thus set δθ max = 5 × 10 −3 equal to the Trotter timestep δt in the AVQDS simulations. This results in a standard deviation of 0.003 for spin correlation observables, and 0.010 for the energy, which is nevertheless about three times bigger than the Trotter energy error. The number of two-qubit gates is the defining factor for practical quantum computing in the NISQ era. As shown in Fig. 4(e) and (f), we compare the number of two-qubit gates contained in the variational ansatz and in the Trotterized circuits as a function of time t. Importantly, the AVQDS simulations require up to two orders of magnitude fewer two-qubit gates than Trotter simulations. In contrast to the linear growth of Trotter circuits, the AVQDS circuits mainly grow in the initial quench stage, and approach a plateau as the system evolves further. For N = 8, we find a plateau value of about 100 CNOTs is sufficient to follow the spin dynamics to the final simulation time, whereas Trotterized circuits require execution of about 10 4 CNOTs. This implies that the variational circuit has gained sufficient expressibility to represent the relevant manifold of the Hilbert space for the linear-ramp dynamics studied here.

D. System-size scaling of circuit complexity
To reach the goal of performing scalable quantum dynamics simulations on quantum devices, it is crucial to estimate how the required quantum resources scale with the system size N . Although this scaling depends on the complexity of the dynamical problem studied and is therefore model dependent, it is instructive to investigate the scaling of the required resources with system size N for the LSM model. As shown in Fig. 5, the final number of parameters of the variational ansatz in AVQDS simulations scale with N 2 . Specifically, we find that the total number of variational parameters, which are associated with single and two-qubit Pauli rotation gates, scales as 1.16 N 2 . The number of variational parameters that are associated with two-qubit Pauli rotation gates scales as 0.78 N 2 .
The AVQDS approach utilizes an inherent Trotter-type structure for the variational ansätze with unitaries constructed based on Hamiltonian Pauli terms. It is conceptually different from the random parameterized quantum circuit optimization approach, where the cost function gradients can become exponentially small with increasing number of qubits [56][57][58]. In the context of adiabatic state preparation, AVQDS can lead the system to the ground state without resorting to explicit high-dimensional optimization or cost function gradients. Therefore, barren plateaus of cost functions associated with random variational circuits are unlikely to constitute a problem for AVQDS simulations. As further numerical evidence, the element-wise maximal absolute value of the vector V defined in Eq. (5), which is closely related to the cost function gradient, remains close to 0.28 in AVQDS simulations of the LSM model as the number of sites is increased from 4 to 8.
The above AVQDS calculations start with the ground state Ψ 0 of the HamiltonianĤ 0 at t = 0 in the FM x phase. Due to the finite transverse field h z , this state is not a tensor-product state and preparation of the initial state is therefore non-trivial. While it is convenient to initialize any state vector in classical simulations, it is not straightforward to prepare an entangled state on quantum computers. As discussed previously in Section II C 2, we here adopt the qubit-ADAPT-VQE method [25] for the initial state preparation. We use a reference product state with all spins aligned in σ z = 1 and an expanded operator pool that includes all one-and two-qubit Pauli terms. For system size N = 6, we find that an ansatz with 20 variational parameters, which are associated with two-qubit Pauli rotation gates, is able to variationally represents |Ψ 0 with unit overlap up to the 7 th decimal place. At the end of the AVQDS simulation, the number of variational parameters associated with two-qubit Pauli rotation gates has increased to 45, which is slightly smaller than a rough estimation of 20 + 0.78N 2 ≈ 20 + 28 for N = 6. Here 20 is the number of variational parameters in the initial qubit-ADAPT-VQE ansatz as mentioned above, and 28 is the number of variational parameters added during time-evolution.

IV. SUDDEN QUENCH DYNAMICS OF THE MIXED-FIELD ISING MODEL
To further benchmark the AVQDS approach for nonintegrable dynamics, we perform sudden quench dynamics withẐ N =Ẑ 0 for periodic boundary conditions. In the following, we measure energy in units of J = 1. This model is integrable for h z = 0, where it becomes the transverse-field Ising model (TFIM), but is nonintegrable when both h z and h x are finite. Initially, the system is prepared in the ordered state |Ψ 0 = |↑ . . . ↑ , which is a ground state of the MFIM in the absence of magnetic fields h x = h z = 0. We consider two sudden quench protocols: (A) quenching to the TFIM with h x = −2.0 at t = 0, and (B) quenching to the MFIM with h x = −2.0, h z = 0.5 at t = 0. Protocol A has been used in the study of dynamical quantum phase transitions [59,60]. Protocol B allows us to compare the performance of AVQDS for integrable and nonintegrable models. The Loschmidt echo, defined as the probability of a time-evolving system to return to its initial state, is a central concept in the study of dynamical quantum phase transitions [59,60]. The simulations presented here starts with the initial state |Ψ 0 = |↑ . . . ↑ , and the Loschmidt echo can be written as L(t) = Ψ 0 |e −iĤ f t |Ψ 0 2 . The time-dependent Loschmidt echo calculated from AVQDS is plotted in Fig. 6 for (a) the integrable TFIM and (d) the nonintegrable MFIM. Both cases are in excellent agreement with exact results. To better quantify the simulation accuracy, the infidelity 1 − f ≡ 1 − | Ψ[θ(t)]| Ψ exact (t) | 2 is shown in panels (b) and (e), indicating that the fidelity f is generally higher than 99.5%. To make a comparison with naive Trotter dynamics, we performed Trotter simulations with a time step size chosen to provide fidelities comparable to those obtained with AVQDS [see green dashed lines in (b) and (e)]. In panels (c) and (f), we plot the number of CNOTs N cx used in the Trotter and variational AVQDS circuits. Consistent with the results of Sec. III, we find that Trotter circuits have about two orders of magnitude more two-qubit gates than the AVQDS circuits. A tendency towards circuit depth saturation becomes noticeable for the integrable TFIM (h z = 0) for simulation times t 2. The parameterized circuit for |Ψ[θ(t)] reaches 134 CNOTs at final simulation time t = 3, slightly larger the N = 8 LSM simulation result. Moving from the integrable TFIM to the nonintegrable MFIM by introducing finite field h z , the AVQDS circuit increases N cx modestly from 134 to 210 at t = 3.
In Fig. 7, we consider the scaling of the number of CNOT gates N cx in the AVQDS circuit with time t (panel a) and system size N (panel b). This is the number of CNOT gates required to build the post-quench state |Ψ[θ(t)] in the MFIM. Importantly, we observe an initially linear growth, N cx ∝ t, that crosses over into saturation at a system size dependent timescale t s (N ). The initial growth resembles the behavior of the entanglement entropy in the post-quench regime [61]. Since the number of CNOTs in the circuit is proportional to the number of variational parameters, this (polynomial) growth is notably slower than the exponential growth of the number of parameters needed in simulations using matrix-product states (MPSs). This is a direct manifestation of the complexity window, which describes the phenomenon that quantum circuits can generate states with high entanglement more efficiently, i.e. with less parameters, than MPSs [62,63]. More precisely, states in the complexity window are highly entangled yet can be represented by a quantum circuit that grows only polynomially in time [64]. Figure 7 thus shows that such efficient circuits can be automatically generated within AVQDS.
Specifically, Fig. 7(a) shows N cx of the AVQDS circuits for the MFIM at system sizes 2 ≤ N ≤ 10 with simulation times up to t = 12, except for N = 9 and 10, which use shorter final times. We observe that N cx for N ≤ 8 saturates within t = 12. The critical time t s for saturation increases with N . Panel (b) depicts N cx as a function of system size N for different fixed times t (vertical cuts of the data in panel (a)). While the saturated circuit depth, as measured by N cx at t = 12, appears to scale beyond quadratically (with N ) for 2 ≤ N ≤ 8 due to the increased complexity of the nonintegrable model, the exact order of the scaling cannot be obtained due to the small size of the data set. Similar quality of fitting is obtained using a power-law function aN α with α ≈ 4.8 and an exponential function b(e N/β − 1) with β ≈ 1.4. Let us now consider the practically relevant question of how the circuit depth grows with N for a sequence of fixed simulation times 1 ≤ t ≤ 7. Because the saturation time t s increases rapidly with N , the equal-time cut generally exhibits a crossover behavior of N cx from the initial super-quadratic growth with small N to approximately linear growth with large N > N c (t) as the circuits cross over into the presaturation regime (where t < t s (N )). The crossover size N c (t) grows slowly with t. This suggests that practical dynamics simulations of generic quantum models out to fixed final times within the AVQDS approach remains scalable with increasing N . Finally, in the inset of Fig. 7(b), we show that N cx after the saturation time t s in the integrable TFIM grows approximately as N 2 , which is slower than the N 4.8 observed for the MFIM case. The quadratic scaling is similar to the results for the LSM chain shown in Fig. 5, suggesting that it is related to the integrability of the TFIM and the LSM models.

V. CONCLUSION
We present a novel adaptive variational approach, AVQDS, to perform quantum dynamics simulations in many-body fermionic and spin models. It builds upon the theory of variational quantum dynamics according to McLachlan's variational principle. The key novelty of the presented AVQDS method is to adaptively and dynamically expand the variational ansatz during time evolution. This allows to account for the changing nature of the time-evolved wavefunction and results in highly accurate variational quantum dynamics simulations. Expansion of the ansatz is controlled by setting a threshold of the maximum allowed McLachlan distance L 2 . This distance describes the difference after one time step between the exact time-evolution of the variational state, as described by the von-Neumann equation, and the time evolution obtained from classically propagating the variational parameters in time. The AVQDS approach does not involve complex optimization in a high-dimensional parameter space.
To benchmark and study the performance of the approach, we use AVQDS to simulate spin dynamics in the integrable LSM model under a finite-rate quantum quench, and in the nonintegrable MFIM under sudden parameter quenches. We consider system sizes up to N = 11 sites. The AVQDS simulations are shown to be in excellent agreement with numerically exact results for local observables, total energy, and wavefunction overlaps. The depth of the resulting AVQDS variational circuits, as characterized by the number of two-qubit CNOT gates N cx , is shown to saturate to values smaller than Trotterized circuits by up to two orders of magnitude for N = 8. We find that N cx after saturation at the end of the linear ramp scales as N cx ∝ N 2 for both the LSM and the TFIM model, suggesting that this polynomial scaling is linked to the integrability of models. For quench dynamics simulations of the nonintegrable MFIM, the saturated number of two-qubit gates in the AVQDS circuits scales as a higher-order polynomial N cx ∝ N 5 approximately, as expected from the increased complexity of the model. For fixed times t < t s (N ), however, we observe that the number of two-qubit gates in the AVQDS circuits reduces to approximately linear growth N cx ∝ N as system size increases, implying the scalability of AVQDS simulations in practice. Finally, we find that the complexity of the AVQDS circuits at fixed N scales initially linearly with time N cx ∝ t, showing that these circuits can efficiently capture the rapid growth of entanglement under nonintegrable dynamics in the system. We envision that the AVQDS approach will have wide applications in the growing field of quantum dynamics and far-from-equilibrium physics. In addition to directly simulating dynamics in other spin and fermionic models, AVQDS can be used as an impurity dynamics solver for quantum embedding approaches for dynamics simulations of large and infinite lattice models [65][66][67][68]. An open question to further explore the finite-size and finitetime scaling of the AVQDS circuit depth and relating it to the entanglement content of the time-evolved state. The prospect of adaptively and automatically generating polynomial depth circuits that generate highly entangled states is intriguing and warrants further investigation. Another important future research direction is to study the noise resilience of the algorithm and noise mitigation strategies, in particular when implementing AVQDS on NISQ QPUs. For the preparation of the initial state of the dynamics simulation, we have explicitly shown that AVQDS can be easily combined with the known qubit-ADAPT-VQE method. Finally, AVQDS can be generalized from real time to imaginary time axis [69], which offers a novel efficient approach to finding ground states of Hamiltonian systems, or to a wider range of optimization problems of static cost functions in the field of machine learning.