Digital quantum simulation, learning of the Floquet Hamiltonian, and quantum chaos of the kicked top

The kicked rotor and the kicked top as periodically driven quantum systems represent paradigmatic models in studying quantum chaos (see chapter 8 in [1]), which have played a central role in the discussion of phenomena like quantum localization and relation to random matrix theory (RMT). In recent collaborative work with Fritz Haake [2] we have pointed out that the well-studied problem of the transition from regular to chaotic dynamics, as observed in the kicked top as a function of the driving frequency, sheds new light on and provides a physical interpretation of Trotter errors in digital quantum simulation (DQS)—a highly relevant problem in the focus of today's effort to 'program' quantum many-body dynamics on quantum computers. In DQS, the unitary evolution U(t) = e^−iHt generated by a Hamiltonian H is simulated by decomposition into a sequence of quantum gates [3]. This can be achieved via a Suzuki–Trotter decomposition [4, 5], which approximately factorizes the time evolution operator in Trotter time steps of size τ (see figure 1(a)). While the 'correct' evolution operator U(t) will emerge in the limit of Trotter stepsize τ → 0, in practice the finite fidelity of quantum gates makes it desireable to take as large Trotter steps as possible [6]. The observation of references [2, 7, 8] was that the error associated with a finite Trotter step size shows a sharp threshold behavior: while for a small time step τ < τ_* Trotterized time evolution provides a faithful representation of the desired dynamics, in the regime τ > τ_* Trotter errors proliferate. This behavior is in correspondence to a transition from regular motion to quantum chaos in Floquet systems when the driving frequency is decreased [9, 10]. While the ultimate goal of DQS is to simulate complex quantum many-body systems with finite range interactions (see recent advances in references [11–23]), the features of a Trotter threshold are already visible in simple models. This leads us back to the kicked top, which—while being intrinsically a single particle problem (for a 'large' spin S)—can also be interpreted as Trotterized time evolution of a many-body spin-model with infinite range interactions. Thus, the kicked top can serve as a testing ground, both theoretically and experimentally, for the phenomenology of the Trotter threshold.

Zoom In Zoom Out Reset image size

In this work, we discuss the Trotter threshold from the perspective of 'Hamiltonian learning (HL)', and we choose the kicked top as a simple model system displaying pertinent features. The discussion is based on the HL framework for the characterization of Trotterized DQS of many-body systems developed by us in reference [24], ⁵ which we adapt and extend here to study the Trotter threshold and the transition to quantum chaos in the kicked top [2]. HL is presently being developed as a new tool in quantum information theory in the context of quantum many-body systems and quantum simulation [25–37]. In the present context, we can phrase HL as follows: we consider quench dynamics of a many-body spin system, |ψ(t)⟩ = e^−iHt |ψ(0)⟩, where the system is initially prepared in the state |ψ(0)⟩ and evolves under a (time-independent) Hamiltonian H to a final state |ψ(t)⟩ at time t. The goal is to learn the operator content of H, i.e., a tomographic reconstruction of H from measurements on |ψ(t)⟩. The key to an efficient learning of H from experimental observations is that a physical many-body Hamiltonian consists of a small (polynomial) number of terms, i.e., the operator content of H will be limited to one-body, quasi-local two-body terms etc, while the many-body wave function lives in a Hilbert space of dimension scaling exponentially with the number of constituents. HL thus becomes efficient by having to learn only a sufficiently small number of coupling coefficients in the Hamiltonians, while testing for presence of additional terms, and thus verifying the learned Hamiltonian structure with more data. In the following we apply these ideas to the kicked top, viewed as DQS of a collective spin system.

The kicked top combines precession of the spin S of the top around the x-axis with τ-periodic non-linear 'kicks' around the z-axis, according to the time-dependent Hamiltonian

$$\begin{equation}{H}_{\text{KT}}(t)={H}_{x}+\tau {H}_{z}\sum\limits _{n\in \mathbb{Z}}\delta (t-n\tau ),\end{equation} \tag{ 1 }$$

where H_x = h_x S_x and ${H}_{z}={J}_{z}{S}_{z}^{2}/(2S+1)$ with quantum angular momentum operators S_μ with μ = x, y, z. The evolution operator generated by H_KT(t) over a single period of duration τ can be equivalently described in terms of a Floquet operator

$$\begin{equation}{U}_{\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{z}\tau }{\mathrm{e}}^{-\mathrm{i}{H}_{x}\tau }\equiv {\mathrm{e}}^{-\mathrm{i}{H}_{\text{F}}(\tau )\tau },\end{equation} \tag{ 2 }$$

as illustrated in figure 1(b). The dynamics of the kicked top is quantum chaotic iff the spectral statistics of U_τ can be described by one of Dysons's ensembles of random matrices [1]. Indeed, while RMT was initially applied in physics by Wigner to understand the distribution of nuclear spectra (see [38] for a review on RMT in nuclear physics), it forms now the basis for the study of quantum chaos. It is indeed a defining feature of quantum chaotic systems that their spectral statistics are universal and obey predictions from RMT [1, 39].

Alternatively equation (2) constitutes the elementary gate sequence of a DQS that aims at approximating the Hamiltonian H = H_z + H_x according to ${\mathrm{e}}^{-\mathrm{i}Ht}\approx {U}_{\tau =t/n}^{n}$, where t denotes the total simulation time which is split into n steps of duration τ = t/n. We emphasize that the accuracy of this approximation does not only depend on the Trotter step τ, but also on the Trotter sequence that can be chosen to compensate Trotter errors up to a given order $\mathcal{O}({\tau }^{k})$ [40, 41]. Equation (2) also defines the Floquet Hamiltonian H_F(τ). For sufficiently small τ, the Floquet Hamiltonian can be written as a FM expansion, i.e., employing Baker–Campbell–Hausdorff formulas,

$$\begin{equation}{H}_{\text{F}}(\tau )={H}_{x}+{H}_{z}+\mathrm{i}\,\frac{\tau }{2}[{H}_{x},{H}_{z}]+\cdots \end{equation} \tag{ 3 }$$

which is a series expansion in the Trotter stepsize τ, written here up to first order, with the higher order terms taking the form of nested commutators of H_x and H_z (see section 2 below). The question of convergence of this series is intrinsically connected with the transition from regular to chaotic dynamics at a specific τ_*. From a DQS point of view, in the limit τ → 0, the Floquet Hamiltonian reduces to H = H_z + H_x , as the desired Hamiltonian 'to be simulated' on a quantum device. Higher order terms in τ represent Trotter errors.

While in traditional discussions of DQS, Trotter errors are quantified in terms of Trotter bounds, ${\Vert}{U}_{\tau }-{\mathrm{e}}^{-\mathrm{i}H\tau }{\Vert}\equiv \mathcal{O}({\tau }^{k})$ [42–44], HL goes beyond in quantifying these errors by learning terms order by order of the FM expansion in an experimentally feasible protocol [24]. Central to our work below is the toolset built into HL which quantifies errors of an Ansatz Hamiltonian, e.g., as a truncated FM series equation (3), to represent the experimental H_F(τ). This will be our key quantifier in studying the Trotter threshold τ_* (see figure 1(c)). Thus our work goes significantly beyond references [2, 7, 8], where the Trotter threshold was studied only for low-order observables. While our discussion below will focus on the kicked top as a simple model system, we emphasize that the main results and conclusions carry over to characterizing Trotter errors in DQS of quantum many-body systems from many-body models in condensed matter physics, to quantum chemistry or high energy physics.

This paper is organized as follows. In section 2 we briefly summarize our previous work [2] on the Trotter threshold of the kicked top. Section 3 provides the main results of the present work. We will start with a description of HL protocols to learn order by order, up to given truncation cutoff K, the Floquet Hamiltonian of the kicked top, followed by a discussion of numerical results illustrating the technique. We conclude with section 4.

In preparation for the discussion in section 3 on HL applied to the Floquet Hamiltonian for the kicked top, we start by reviewing previous work, and, in particular, our collaborative work with Fritz Haake [2].

2.1. Quantum many-body models and the kicked top

The motivation for the present study is the quantitative characterization of Trotter errors, and the Trotter threshold in particular, in quantum many-body systems. Therefore, we find it useful to recall some of the basic features of DQS of quantum many-body systems, which we illustrate here for 1D spin models, in relation to the kicked top as the model system studied below.

Typical model systems of interest are one-dimensional chains of N spin-1/2 such as the long-range Ising model with Hamiltonian

$$\begin{equation}H={H}_{x}+{H}_{z},\quad {H}_{x}={h}_{x}\sum\limits _{i=1}^{N}{\sigma }_{i}^{x},\quad {H}_{z}={J}_{z}\sum\limits _{i< j=1}^{N}\frac{{\sigma }_{i}^{z}{\sigma }_{j}^{z}}{{\left\vert i-j\right\vert }^{\alpha }},\end{equation} \tag{ 4 }$$

where ${\sigma }_{i}^{\mu }$ with μ = x, y, z are Pauli operators for spins on lattice sites i = 1, ..., N. Power-law interactions with 0 ⩽ α ⩽ 3 are routinely implemented with trapped-ion quantum simulators [45]. The kicked top emerges in DQS in the limit α → 0.

Before proceeding to a study of the kicked top, we recall some of the basic features of Trotter errors and the Trotter threshold, which have emerged in our previous work. In DQS, time evolution generated by H is represented by a sequence of elementary quantum gates. This is achieved through the approximate factorization of the time evolution operator within each Trotter step, ${\mathrm{e}}^{-\mathrm{i}H\tau }\approx {\mathrm{e}}^{-\mathrm{i}{H}_{z}\tau }{\mathrm{e}}^{-\mathrm{i}{H}_{x}\tau }$: individual terms in the sums in H_x and H_z , respectively, commute with each other, so that ${\mathrm{e}}^{-\mathrm{i}{H}_{x}\tau }$ and ${\mathrm{e}}^{-\mathrm{i}{H}_{z}\tau }$ can directly be decomposed into single-spin and two-spin gates. However, due to the non-commutativity of the components of the target Hamiltonian, [H_x , H_z ] ≠ 0, the factorization of the time evolution operator within a single Trotter step is exact only in the limit τ → 0, and any finite value τ > 0 leads to the occurrence of Trotter errors. On the level of the generator of the dynamics, the FM expansion equation (3) suggests that Trotter errors are perturbatively small in τ. However, a rigorous sufficient condition for the convergence of the FM expansion [46–48] indicates that the radius of convergence scales with system size as $\sim 1/N$, which would imply that the FM expansion is not applicable in the thermodynamic limit N → ∞, and puts into question whether Trotter errors can be controlled in DQS of quantum many-body systems. Addressing this question requires a suitable measure of Trotter errors. In DQS, quantities of physical interest are typically expectation values of few-body observables, i.e., (sums of) products of spin-1/2 operators acting on only a few different spins, such as ${\sigma }_{i}^{\mu }$, ${\sigma }_{i}^{\mu }{\sigma }_{j}^{\nu }$, ${\sigma }_{i}^{\mu }{\sigma }_{j}^{\nu }{\sigma }_{k}^{\rho }$, etc. Therefore, it is natural to quantify Trotter errors in terms of deviations of expectation values of few-body observables from their target values that are obtained by time evolution generated by the target Hamiltonian H. For Ising spin chains, in the limit α → ∞ of short range interactions and for 0 ⩽ α ⩽ 3, such quantitative studies of Trotter errors of few-body observables were carried out in references [2, 7], respectively. As we will illustrate with a concrete example below, these studies found sharp threshold behavior, where Trotter errors remain controlled for small Trotter steps and proliferate for τ larger than a threshold value τ_*.

The Trotter threshold was observed consistently over the entire range of values of power-law interaction exponents α considered. In the limit α → 0, DQS of the spin model in equation (4) is directly related to the dynamics of a kicked top: for α = 0, the components of the Ising Hamiltonian in equation (4) can be cast, in terms of collective spin operators ${S}_{\mu }={\sum }_{i=1}^{N}\;{\sigma }_{i}^{\mu }/2$, as H_x ∼ S_x and ${H}_{z}\sim {S}_{z}^{2}$—just as in the Hamiltonian of the kicked top in equation (1). Then, the collective spin ${\boldsymbol{S}}^{2}={S}_{x}^{2}+{S}_{y}^{2}+{S}_{z}^{2}$ becomes a constant of motion. Consequently, the many-body Hilbert space with dimension 2^N is decomposed into decoupled subspaces of fixed total spin S, and within each subspace, the Trotterization of equation (4) reproduces the dynamics of a kicked top of size S. In this sense, the kicked top becomes a single-particle (i.e., a single collective spin) toy model for the many-body Trotter threshold, where the notion of few-body observables introduced above to quantify Trotter errors translates to low-order products of spin operators S_μ . In the following, we give a detailed account of the Trotter threshold in the context of the kicked top.

2.2. The kicked top and the Trotter threshold

2.2.1. Model system and Magnus expansion

Before we enter a detailed discussion of the Trotter threshold, we introduce the following extension of the Floquet operator given in equation (2):

$$\begin{equation}{U}_{\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{z}\tau }{\mathrm{e}}^{-\mathrm{i}{H}_{y}\tau }{\mathrm{e}}^{-\mathrm{i}{H}_{x}\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{\text{F}}(\tau )\tau }.\end{equation} \tag{ 5 }$$

This extended Floquet operator, on which our discussion of the Trotter threshold will focus, corresponds to DQS of a target Hamiltonian H = H_x + H_y + H_z , where

$$\begin{equation}{H}_{\mu }=\frac{{J}_{\mu }{S}_{\mu }^{2}}{2S+1}+{h}_{\mu }{S}_{\mu }.\end{equation} \tag{ 6 }$$

The rationale behind this choice of model, as opposed to the model for the kicked top in the introduction, is explained in appendix A: the key difference between the Floquet operators in equations (2) and (5) are the absence of time-reversal and geometrical symmetries as well as resonant driving points in the latter case [1]. We choose J_z as the unit of energy, and, for concreteness, we fix J_y = 0, J_x = 0.4J_z , h_z = h_y = 0.1J_z , and h_x = 0.11J_z . Further, we note that while the FM expansion always takes the form of a power series in the Trotter step size τ,

$$\begin{equation}{H}_{\text{F}}(\tau )=\sum\limits _{k=0}^{\infty }{\tau }^{k}\,{C}_{k},\end{equation} \tag{ 7 }$$

for the specific case of a three-step Floquet drive in equation (5), the first few terms are given by

$$\begin{equation}{C}_{0}=\sum\limits _{\alpha }{H}_{\alpha },\end{equation} \tag{ 8 }$$

$$\begin{equation}{C}_{1}=\frac{\mathrm{i}}{2}\sum\limits _{\alpha < \beta }[{H}_{\alpha },{H}_{\beta }],\end{equation} \tag{ 9 }$$

$$\begin{equation}{C}_{2}=-\sum\limits _{\alpha \ne \beta }\frac{[{H}_{\alpha },[{H}_{\alpha },{H}_{\beta }]]}{12}-\frac{[{H}_{x},[{H}_{y},{H}_{z}]]}{6}-\frac{[{H}_{z},[{H}_{y},{H}_{x}]]}{6},\end{equation} \tag{ 10 }$$

where α, β ∈ {x, y, z}.

2.2.2. Quantifying Trotter errors

Among different choices of few-body observables to quantify Trotter errors, a special role is played by the target Hamiltonian itself: for τ → 0, the Floquet operator in equation (5) reduces to time evolution generated by the time-independent target Hamiltonian H, and the energy as measured by the expectation value of H becomes a constant of motion. To quantify the degree to which this conservation law is obeyed in DQS, we define the simulation accuracy [7, 49]:

$$\begin{equation}{Q}_{\text{E}}(n\tau )=\frac{{E}_{\tau }(n\tau )-{E}_{0}}{{E}_{T=\infty }-{E}_{0}}.\end{equation} \tag{ 11 }$$

Energy is conserved if Q_E(nτ) = 0, and the energy at time t = nτ, given by ${E}_{\tau }(n\tau )=\langle \psi (0)\vert {U}_{\tau }^{{\dagger}n}H{U}_{\tau }^{n}\vert \psi (0)\rangle$, equals the energy of the initial state, E₀ = ⟨ψ(0)|H|ψ(0)⟩. In contrast, the value Q_E(nτ) = 1 indicates that the system absorbs energy from the time-periodic Floquet drive and heats up to infinite temperature T = ∞ such that ${E}_{\tau }(n\tau )={E}_{T=\infty }=\mathrm{tr}(H)/\mathcal{D}$ where $\mathcal{D}=2S+1$ is the Hilbert space dimension.

2.2.3. Trotter threshold

The temporal average $\bar{{Q}_{\text{E}}}(t)=\frac{1}{{n}_{t}}{\sum }_{n=1}^{{n}_{t}}\;{Q}_{\text{E}}(n\tau )$, where ${n}_{t}=\lfloor t/\tau \rfloor$ is the number of Trotter steps corresponding to the simulation time t, is shown in figure 2(a). Here, the initial state is chosen as a spin coherent state, $\vert \theta ,\phi \rangle ={\mathrm{e}}^{\mathrm{i}\theta \left({S}_{x}\mathrm{sin}(\phi )-{S}_{y}\mathrm{cos}(\phi )\right)}\vert S,{S}_{z}=S\rangle$, with θ = 0.1 and ϕ = 0.2. At times J_z t  ≳  20, Trotter errors, quantified by the time-averaged simulation accuracy $\bar{{Q}_{\text{E}}}(t)$, exhibit threshold behavior: while $\bar{{Q}_{\text{E}}}(t)$ increases smoothly for small values of τ, a sudden jump to the saturation value $\bar{{Q}_{\text{E}}}(t)\approx 1$ occurs at J_z τ_* ≈ 3.5. The Trotter threshold persists for t → ∞, and becomes sharper with increasing spin size S [2]. An analogous Trotter threshold can be observed also for other few-body observables [2], and, as pointed out above, applies also to one-dimensional spin chains with algebraic [2] and nearest-neighbor [7] interactions, as well as bosonic models with infinite Hilbert spaces [8]. It should be noted that generic and sufficiently short-range interacting many-body systems are expected to heat up indefinitely when subjected to Floquet driving [9, 51, 52], so that the Trotter threshold is washed out with $\bar{{Q}_{\text{E}}}(t)\to 1$ for any τ > 0 at late times. However, the heating rate is exponentially small in the driving frequency [53–59]. Thus, for given finite simulation run-time, heating can be suppressed efficiently.

Zoom In Zoom Out Reset image size

Below the Trotter threshold, Trotter errors remain controlled, i.e., perturbatively small in τ. In fact, Trotter errors are well-described by a low-order truncation of the FM expansion: treating corrections to the target Hamiltonian in equation (7) in time-dependent perturbation theory [2, 7], the long-time average of the simulation accuracy, $\bar{{Q}_{\text{E}}}={\mathrm{lim}}_{t\to \infty }\,\bar{{Q}_{\text{E}}}(t)$, can be written as a power series in τ:

$$\begin{equation}\bar{{Q}_{\text{E}}}={q}_{1}\tau +{q}_{2}{\tau }^{2}+\mathcal{O}({\tau }^{3}),\end{equation} \tag{ 12 }$$

with, for the parameters chosen for figure 2(a), q₁ ≈ 0.007 and q₂ ≈ 0.035. The good agreement between this series expansion and the numerical data for τ ≲ τ_* shows that a low-order truncation of the FM expansion provides a quantitatively accurate description of the dynamics of few-body observables well beyond rigorous bounds on the radius of convergence of the FM expansion.

2.2.4. Breakdown of the FM expansion and onset of quantum chaos

The breakdown of the FM expansion, signaled by the proliferation of Trotter errors for τ > τ_*, marks the onset of quantum chaos. As a first indication for the connection between the divergence of the FM expansion and quantum chaos, we note that the saturation value $\bar{{Q}_{\text{E}}}=1$ at τ  ≳  τ_* is consistently reproduced by replacing U_τ in each Trotter step by a random unitary matrix. Then, the temporal average is manifestly equivalent to a Hilbert-space average, which by definition yields the infinite-temperature value E_T=∞.

More systematically, quantum chaos of the Floquet operator U_τ manifests in statistics of the eigenphases θ_n and in localization properties of the eigenvectors |ϕ_n ⟩, which obey the eigenvalue equation ${U}_{\tau }\!\vert {\phi }_{n}\rangle ={\mathrm{e}}^{\mathrm{i}{\theta }_{n}}\!\vert {\phi }_{n}\rangle$. In particular, as a defining signature of quantum chaos, the distribution of spacings of eigenphases θ_n is described by RMT, and the eigenvectors |ϕ_n ⟩ are delocalized in a basis of eigenvectors of an operator that generates integrable dynamics, such as the target Hamiltonian. In the following, we study the onset of chaos in the kicked top via comparison of the eigenphases and eigenstates of the Floquet operator to their respective RMT predictions. We refer the reader to reference [1] for a comprehensive introduction to the signatures of quantum chaos.

Eigenphase statistics are conveniently characterized by the average adjacent phase spacing ratio r [60],

$$\begin{equation}r=\frac{1}{\mathcal{D}}\sum\limits _{n=1}^{\mathcal{D}}{r}_{n},\quad {r}_{n}=\frac{\mathrm{min}({\delta }_{n},{\delta }_{n+1})}{\mathrm{max}({\delta }_{n},{\delta }_{n+1})},\end{equation} \tag{ 13 }$$

where δ_n = θ_n+1 − θ_n . This quantity constitutes a measure of the degree of repulsion between the eigenphases of U_τ , and takes characteristic values for unitaries that are drawn from an ensemble of RMT. In particular, for the circular unitary ensemble (CUE) of random unitary matrices, the average adjacent phase spacing ratio is given by r_CUE = 0.5996(1) [61]. Indeed, this value is reproduced for τ  ≳  τ_* as shown in figure 2(b), giving a clear indication for quantum-chaotic dynamics beyond the Trotter threshold. At very small values J_z τ ≲ 1/S, the phase spacing ratio is determined by the target Hamiltonian. The phase spacing ratio drops to r_POI ≈ 0.39 for 1/S ≲ J_z τ ≲ J_z τ_*, which signals the absence of level repulsion in crossings of eigenphases that wind repeatedly around the unit circle, leading to Poisson statistics for the adjacent phase spacings [2]. In this regime, Trotterization leads to only weak mixing between the eigenvectors |ψ_n ⟩ of the target Hamiltonian, and the Floquet states |ϕ_n ⟩ are in one-to-one correspondence with the states |ψ_n ⟩. This localization of Floquet states in the eigenbasis of the target Hamiltonian is quantified by the participation ratio (PR),

$$\begin{equation}\mathrm{P}\mathrm{R}={\left(\sum\limits _{n,m=1}^{\mathcal{D}}\vert \langle {\psi }_{n}\vert {\phi }_{m}\rangle {\vert }^{4}\right)}^{-1}.\end{equation} \tag{ 14 }$$

The PR is shown in figure 2(c). The small values of PR at τ ≲ τ_* indicate a high degree of similarity between the eigenbases of H and U_τ . For τ  ≳  τ_*, the PR saturates to a large value that indicates equal absolute overlaps between all eigenstates, as expected for the eigenvectors |ϕ_n ⟩ of a random matrix U_τ . In this regime, the precise numerical value of the PR is determined by the corresponding ensemble of RMT. For the CUE, and in the limit $\mathcal{D}\to \infty$, it is given by PR_CUE = 1/2 [50].

Previous work, as summarized in section 2, has focused on the Trotter threshold by monitoring the simulation accuracy Q_E, PR, and the level spacing ratio of the Floquet operator U_τ . Instead, the present section will study the phenomenology of the Trotter threshold, and the transition to quantum chaos, from the view point of learning the Floquet Hamiltonian H_F(τ) in the form of the FM expansion equation (7), as a function of τ.

Below we first describe the technique of HL (section 3.1), which we then adapt to 'learn' the FM expansion of the kicked top. Central to our discussion is the ability of the HL protocol to provide a quantitative assessment of errors in HL, in particular in learning H_F(τ) with τ approaching the Trotter threshold τ_*. Corresponding results for the kicked top will be presented in section 3.2.

3.1. Hamiltonian learning

3.1.1. HL protocols for quantum many-body systems

Recently, HL has been developed as a technique to efficiently recover an unknown Hamiltonian of an isolated quantum many-body system via measurements on quantum states prepared in the laboratory. Motivation for HL is provided by the ongoing development of controlled quantum many-body systems, e.g., as analog quantum simulators, where HL serves to characterize and thus verify the functioning of quantum devices. Various methods for recovering the Hamiltonian have been described, from learning the Hamiltonian from a single eigenstate [25–28], or stationary states [29, 32], or from observation of short time dynamics [31, 33].

To be specific, we give a brief description of HL in quench dynamics. We focus on HL as proposed in reference [31] as our discussion below, on learning the Floquet Hamiltonian H_F(τ), directly builds on this method. Reference [31] considers quench dynamics of an isolated many-body system, |ψ(t)⟩ = e^−iHt |ψ(0)⟩, where an initial state |ψ(0)⟩ evolves in time t according to a time-independent Hamiltonian H to a final state |ψ(t)⟩. The task is to recover the operator structure of H by measuring observables at various times t. The key to an efficient reconstruction is the locality of the physically implemented many-body Hamiltonian. This means that H can typically be expanded as a sum of few-body operators, and can thus be specified by a number of coupling parameters that scales at most polynomially in system size N. In the specific example of the long-range interacting spin-1/2 model discussed in section 2 (see equation (4)), of which the kicked top can be seen as the limit α → 0, the Hamiltonian is expressed as a sum of Pauli operators ${\sigma }_{i}^{z}{\sigma }_{j}^{z}$ with coefficients ${J}_{i,j}^{z}\propto \frac{{J}_{z}}{\vert i-j{\vert }^{\alpha }}$, and depends thus on $\mathcal{O}({N}^{2})$ coupling coefficients, far less than the $\mathcal{O}({4}^{N})$ parameters needed to express a generic operator.

We thus seek to reconstruct the physically implemented Hamiltonian from an Ansatz

$$\begin{equation}H(\boldsymbol{c})=\sum\limits _{{h}_{j}\in \mathcal{A}}{c}_{j}{h}_{j},\end{equation} \tag{ 15 }$$

specified by a chosen set $\mathcal{A}={\left\{{h}_{j}\right\}}_{j=1}^{{N}_{\mathcal{A}}}$ of ${N}_{\mathcal{A}}=\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}(N)$ few-body Ansatz operators h_j , and depending on coefficients c_j , with c denoting the vector of such coefficients, which we want to determine from experimental measurements. In the context of quantum simulation, the choice of the Ansatz is based on the target Hamiltonian one seeks to implement on the quantum device, possibly complemented with additional terms representing deviations from this target Hamiltonian, whose presence we want to test. Reference [31] proposes to determine the coefficients c_j that yield the best approximation to the implemented Hamiltonian from the condition of energy conservation:

$$\begin{equation}\langle \psi (0)\vert H\vert \psi (0)\rangle =\langle \psi (t)\vert H\vert \psi (t)\rangle \quad \forall \,\,\vert \psi (0)\rangle ,t.\end{equation} \tag{ 16 }$$

Imposing energy conservation at the level of the Ansatz H( c ) amounts to choosing ${N}_{\text{con}} > {N}_{\mathcal{A}}$ different initial states ${\left\{\vert {\psi }_{i}(0)\rangle \right\}}_{i=1}^{{N}_{\text{con}}}$ as 'constraints', and to minimize the energy differences |⟨ψ_i (0)|H( c )|ψ_i (0)⟩ − ⟨ψ_i (t)|H( c )|ψ_i (t)⟩| w.r.t. the coefficients. The optimal c , which we denote with c ^rec in the following, is thus determined as [31]

$$\begin{equation}{\boldsymbol{c}}^{\text{rec}}=\mathrm{arg}\,\underset{\boldsymbol{c}}{\mathrm{min}}\frac{{\Vert}M\boldsymbol{c}{\Vert}}{{\Vert}\boldsymbol{c}{\Vert}},\end{equation} \tag{ 17 }$$

where the elements of the constraint matrix M are given by

$$\begin{equation}{M}_{i,j}=\langle {\psi }_{i}(t)\vert {h}_{j}\vert {\psi }_{i}(t)\rangle -\langle {\psi }_{i}(0)\vert {h}_{j}\vert {\psi }_{i}(0)\rangle ,\end{equation} \tag{ 18 }$$

which can be inferred from experimental measurements collected from a series of quantum quenches starting from the chosen states ${\left\{\vert {\psi }_{i}(0)\rangle \right\}}_{i=1}^{{N}_{\text{con}}}$. The reconstructed Hamiltonian, denoted with H_rec ≡ H( c ^rec), is the one that best approximates the implemented H within the Ansatz space spanned by $\mathcal{A}$, in the sense of energy conservation. Comparison of H_rec with the target Hamiltonian allows one to assess the quality of the quantum device in simulating the model of interest. Note that in equation (17) the overall scale of c ^rec remains undetermined, but there exist efficient ways of determining it (see references [29, 32]).

It is essential that the technique above also provides a way of assessing errors, or confidence in the reconstructed parameters c ^rec for a given set of experimental data, which are provided here as correlation functions entering the constraint matrix M. The quantity that one has access to from HL is the optimal value of the cost function ||M c || in equation (17),

$$\begin{equation}{\lambda }_{1}=\frac{{\Vert}M{\boldsymbol{c}}^{\text{rec}}{\Vert}}{{\Vert}{\boldsymbol{c}}^{\text{rec}}{\Vert}}=\frac{\sqrt{{\sum }_{i}{\left\vert \langle {H}_{\text{rec}}{\rangle }_{i,t}-\langle {H}_{\text{rec}}{\rangle }_{i,0}\right\vert }^{2}}}{{\Vert}{\boldsymbol{c}}^{\text{rec}}{\Vert}},\end{equation} \tag{ 19 }$$

with ${\langle {H}_{\text{rec}}\rangle }_{i,t}\equiv \langle {\psi }_{i}(t)\vert {H}_{\text{rec}}\vert {\psi }_{i}(t)\rangle$, which measures how well the reconstructed Hamiltonian is conserved during the dynamics governed by H. In practice, λ₁ is calculated as the smallest singular value of M, which is shown to be equivalent to the minimum of the 'energy cost-function' ||M c || in equation (17) (see reference [31] and appendix B for details). In section 3.2 below, λ₁ will play a key role in our study of the Trotter threshold.

3.1.2. Learning the generator of a Trotter step of the kicked top 'order by order' in the FM expansion

We now extend the HL protocol described above to the learning of the Floquet Hamiltonian of the kicked top, defined via ${U}_{\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{\text{F}}(\tau )\tau }$ with U_τ given in equation (5). As outlined in the introduction, the kicked top constitutes a particular example of more general DQS of genuinely many-body systems (see section 2.1), where the HL method proposed here would acquire its most (experimentally) relevant application.

In this context, we choose to rephrase the kicked top as an 'experiment' implementing a Trotter evolution cycle which is presented to us as a black box. We denote the implemented Trotter cycle as ${U}_{\tau }^{\mathrm{exp}}={\mathrm{e}}^{-\mathrm{i}{H}_{\text{F}}^{\mathrm{exp}}(\tau )\tau }$, parametrized in terms of the generator ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$. In an ideal experiment ${H}_{\text{F}}^{\mathrm{exp}}(\tau )={H}_{\text{F}}(\tau )$, but here we wish to infer the operator content of ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ from experimental measurements, for two main reasons. First, in DQS we are interested in learning the Trotter errors, contained in ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ as higher order terms in τ (see FM expansion equation (7)). Second, in an experimental context there might also be control errors, which will be reflected in the operator structure of ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ deviating from the ideal H_F(τ): the reconstruction of ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ enables the detection of such experimental control errors, and thus a characterization of the Trotter block which we can view as a process tomography of the corresponding quantum circuit ⁶ (see figure 1).

To achieve these goals, we seek for a reconstruction of ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ based on the operator Ansatz

$$\begin{equation}{H}_{\text{F}}^{\mathrm{exp}(K)}(\tau )=\sum\limits _{k=1}^{K}{\tau }^{k}\,{C}_{k}^{\mathrm{exp}},\end{equation} \tag{ 20 }$$

formulated as series in τ truncated at order K, i.e., we learn the generator 'order by order' by increasing stepwise K. The operator content of the operators ${C}_{k}^{\mathrm{exp}}$ is specified based on our expectation H_F(τ) (see the C_k in equations (8)–(10)), but might also contain additional terms whose presence we want to test. Using the notation of the previous subsection, we may rewrite the Ansatz as

$$\begin{equation}{H}_{\text{F}}^{\mathrm{exp}(K)}(\tau )=\sum\limits _{{h}_{j}\in {\mathcal{A}}_{K}}{c}_{j}\,{h}_{j},\end{equation} \tag{ 21 }$$

with ${\mathcal{A}}_{K}$ denoting the Ansatz set comprising the operator content of the ${C}_{k}^{\mathrm{exp}}$ up to order K in τ. In this notation, ${\mathcal{A}}_{0}$ is the Ansatz for the operator content of the zeroth order terms, corresponding to the generators of the individual gates in ${U}_{\tau }^{\mathrm{exp}}$ (e.g., the H_μ in equation (5)), while ${\mathcal{A}}_{k > 0}$ contains all operators generated by the k-nested commutators of the terms in ${\mathcal{A}}_{0}$, with k = 1, ..., K.

The key quantity representing the quality of our reconstruction is λ₁(τ) introduced in equation (19), which quantifies the error of our—typically low order in τ—Ansatz in representing the 'experimental' Trotter block. The structure of the Ansatz (21) and the ability of measuring λ₁(τ) for different values of τ allows us to (i) discriminate between control errors (appearing as 'unwanted' terms in ${\mathcal{A}}_{0}$) and Trotter errors (the higher order terms), and (ii) detect the value of τ_* at which these Trotter errors proliferate and quantum chaotic dynamics emerges, corresponding to the stepsize at which our Ansatz (20) fails in approximately capturing the stroboscopic DQS dynamics. This is in analogy to the schematic figure 1(c) where λ₁(τ) will become the proxy for the Hamiltonian distance.

Below, we will apply these ideas to the study of the Trotter threshold in the kicked top, using the behavior of λ₁(τ) as the quantifier signaling the transition from regular to chaotic dynamics when approaching τ → τ_*.

3.2. Trotter threshold from HL: results for the kicked top

We now present our results for learning the generator ${H}_{\text{F}}^{\mathrm{exp}}(\tau )$ of the Trotter step as a function of τ via HL, as a method to detect and characterize the Trotter threshold in an experimentally feasible protocol. Specifically, we simulate the above protocol for learning the generator of the kicked top dynamics, imagined here as a DQS implementing the Trotter cycle U_τ defined in equation (5). In order to illustrate the main ideas, and to compare to the results of section 2, we consider U_τ without the addition of 'unwanted' terms (experimental control errors), i.e., ${H}_{\text{F}}^{\mathrm{exp}}(\tau )={H}_{\text{F}}(\tau )$, and simulate the HL protocol in absence of measurement noise in the matrix elements M_i,j of equation (18). To reconstruct H_F(τ), the Ansatz operators h_j for the HL protocol are chosen such that they capture the operator content of the first few orders of the FM expansion in equation (7). We denote with ${\mathcal{A}}_{k}$ the set of operators ${\left\{{h}_{j}\right\}}_{j}$ corresponding to a kth-order truncation of the FM series. For the kicked top of the present example, the Ansatz sets for the first few orders read as ⁷

$$\begin{align*}\hfill & {\mathcal{A}}_{0}=\left\{{S}_{x}^{2},{S}_{x},{S}_{y},{S}_{z}^{2},{S}_{z}\right\},\hfill \\ \hfill & {\mathcal{A}}_{1}={\mathcal{A}}_{0}\cup \left\{{S}_{x}{S}_{y},{S}_{y}{S}_{z},{S}_{x}{S}_{z},{S}_{x}{S}_{y}{S}_{z}\right\},\hfill \\ \hfill & {\mathcal{A}}_{2}={\mathcal{A}}_{1}\cup \left\{{S}_{x}^{2}{S}_{y},{S}_{y}{S}_{z}^{2},{S}_{x}^{2}{S}_{z}^{2},{S}_{z}^{4},{S}_{x}^{2}{S}_{z},{S}_{z}^{3},{S}_{x}{S}_{z}^{2},{S}_{x}^{4},{S}_{x}^{3}\right\},\hfill \end{align*}$$

containing all the linearly independent products of spin operators coming from the commutators in equation (8) ⁸ . The elements M_i,j of the constraint matrix M defined in equation (18) are determined from expectation values of these operators over the states $\vert {\psi }_{i}(n\tau )\rangle ={U}_{\tau }^{n}\vert {\psi }_{i}(0)\rangle$. To compare to the results of section 2, we choose $n\tau =100{J}_{z}^{-1}$ and random coherent states as initial states: our results are however independent of these specific choices. The Trotter threshold discussed in section 2 will be revealed by measuring λ₁ for several values of the Trotter step τ, keeping the Ansatz ${\mathcal{A}}_{k}$ fixed as τ is changed. Specifically, the behavior of λ₁(τ) serves as an indicator of the quality of HL, which is bound to fail for τ > τ_* where the dynamics of the kicked top cannot be approximated by a truncated FM expansion.

The behavior of λ₁ as a function of τ is exemplified in figure 3(a) for several ${\mathcal{A}}_{k}$. Two different behaviors are clearly visible, separated by a grey vertical line denoting the estimated value of τ_*, corresponding to the Trotter threshold discussed in section 2. Let us now describe these regimes and relate them to the observations presented in section 2.

Zoom In Zoom Out Reset image size

3.2.1. Pre-threshold regime

For τ < τ_* we see that, depending on the Ansatz, λ₁ obeys different scaling with τ, i.e., ${\lambda }_{1}\in \mathcal{O}({\tau }^{k+1})$ for a kth-order Ansatz ${\mathcal{A}}_{k}$. This behavior of λ₁(τ) confirms that for τ < τ_* a low-order truncation of the FM expansion is sufficient to accurately describe the stroboscopic dynamics of the kicked top. Indeed, it is intuitively clear that an Ansatz ${\mathcal{A}}_{k}$ capturing only the lowest k orders of H_F(τ) results in a reconstructed Hamiltonian ${H}_{\text{rec}}(\tau )={\sum }_{j}\;{c}_{j}^{\text{rec}}(\tau ){h}_{j}$ violating energy conservation by terms of order $\mathcal{O}({\tau }^{k+1})$ during the stroboscopic dynamics, i.e., ${\Vert}{H}_{\text{rec}}(\tau )-{H}_{\text{F}}(\tau ){\Vert}\in \mathcal{O}({\tau }^{k+1})$. Since λ₁ is precisely what captures how well H_rec(τ) is conserved (see equation (19)), the scaling of ||H_rec(τ) − H_F(τ)|| must be reflected in that of λ₁(τ), hence ${\lambda }_{1}\in \mathcal{O}({\tau }^{k+1})$. This is also confirmed by the behavior of the parameter distance || c ^FM − c ^rec||, shown by the dark red line in figure 3(b), where c ^FM denotes the vector of coefficients that are calculated from the analytical FM expansion (truncated to order k for an Ansatz ${\mathcal{A}}_{k}$). As long as τ < τ_*, || c ^FM − c ^rec|| remains small and its behavior as a function of τ reflects that of λ₁(τ).

3.2.2. Trotter threshold

At τ ≈ τ_*, λ₁(τ) transitions from a $\mathcal{O}({\tau }^{k})$ behavior to a τ-independent value. The reason of this transition is the fact that a low-order truncation of H_F(τ) becomes insufficient to accurately describe the dynamics. More precisely, since our HL protocol is based on monitoring the dynamics of few-body observables (the ${h}_{j}\in {\mathcal{A}}_{k}$) which are measured in order to construct the matrix M, the onset of the plateau in λ₁(τ) signals the proliferation of Trotter errors in the dynamics of such observables. The behavior of λ₁(τ) is therefore analogous to that of the simulation accuracy $\bar{{Q}_{\text{E}}}(\tau )$, defined in equation (11), showing the threshold behavior at the same value τ_*. In fact, $\bar{{Q}_{\text{E}}}(\tau )$ and λ₁(τ) encode very similar physical information, in that they both are related to how well energy, as defined by the time-averaged Hamiltonian in one case and by a low-order truncation of FM expansion in the other, is conserved. Specifically, it is clear from its definition (equation (19)) that λ₁(τ) can be seen as a simulation accuracy Q_E(nτ) (equation (11)) averaged over several initial states, upon identifying the energy E_τ with the reconstructed Hamiltonian H_rec(τ) (up to an overall normalization). We notice that threshold behavior can be observed also in the parameter distance defined before, and shown in figure 3(b), as well as in the reconstructed Hamiltonian parameters which, while showing almost perfect agreement with the analytically determined FM expansion for τ < τ_*, become essentially random for τ > τ_*. We show the behavior with τ of some of the reconstructed parameters in figure 3(c). Furthermore, we observe that the value of τ_* does not strongly depend on the spin size S, as we show in the inset of figure 3(a). This is consistent with the results of our previous work [2], where we observed that τ_* remained finite even in the classical limit S → ∞.

3.2.3. Post-threshold regime

For τ > τ_* we observe a constant value of λ₁(τ), independently of the Ansatz chosen. As shown in section 2, the regime τ > τ_* corresponds to the chaotic regime where the spectral statistics of U_τ is faithful to RMT. As a consequence, the value of λ₁(τ > τ_*) shows perfect agreement with the value of λ₁ that can be estimated using RMT. In this sense, λ₁(τ) can be seen as an indicator of the onset of chaos, showing a behavior similar to the level spacing ratio r and the PR introduced in section 2. Here, the RMT estimate of λ₁ has to be understood as the value that one would obtain when carrying out the HL protocol above with the same Ansatz ${\mathcal{A}}_{k}$, but replacing U_τ with a random U drawn from the relevant matrix ensemble (CUE for the example presented here). We refer the reader to appendix C for the details of this calculation. In the plots, only the CUE estimate for the zeroth order Ansatz is shown, corresponding to the horizontal grey dashed line in panels (a) and (b): changing the Ansatz to reconstruct higher orders of the FM expansion leads to different values for the CUE estimate, which match with the observed plateaus in the plots.

In the present paper, we have shown how HL can be employed as a method for analyzing the transition to chaotic dynamics in the kicked top, interpreted here as a DQS of a collective spin system where this transition manifests itself in the proliferation of Trotter errors for a certain Trotter step size τ_* [2]. In general, our method gives a recipe for interrogating (via measurements) a quantum device implementing a Trotterized evolution, asking the question: can the implemented dynamics be approximately described by a time-independent few-body Hamiltonian? A low-order truncation of the FM expansion constitutes an example of such a Hamiltonian, whose structure motivates our choice of the Ansatz for HL, and thus of the measurements performed to interrogate the system. The answer to the question above is provided by λ₁, assessing the quality of the Ansatz: by measuring λ₁ as a function of τ one can detect the value τ_* at which a description in terms of a few-body Hamiltonian breaks down, corresponding to the proliferation of Trotter errors in DQS.

When applied to the kicked top, our protocol requires only the measurement of products of few collective spin operators, and may be readily implemented in experimental setups where the kicked top dynamics can be realized. These include existing realizations of the kicked top as the spin of single atoms [62], nuclear spins [63], composite spins in nuclear magnetic resonance [64], and collective spins of ensembles of superconducting qubits [65], together with potential realizations with magnetic atoms [66, 67] and trapped ions [68, 69].

We emphasize that the method presented here applies also to the case of DQS of genuinely many-body systems [24], of which the kicked top can be seen as a limit where the dynamics, for suitable initial conditions, admits a description in terms of collective spin operators (see section 2.1). In view of the experimental applicability of our HL protocol in this more general context, we conclude by commenting on its relation to the previously proposed signatures of Trotter threshold and chaos transition in Trotterized dynamics.

First, while we have pointed out a similarity between the physical information encoded in λ₁ and the simulation accuracy Q_E (see section 3.2.2), we note that there exist a fundamental difference between the two. The definition of the simulation accuracy is reliant on the exact knowledge of the implemented Hamiltonian which, in an experimental setting, may be unknown or only approximately known. Conversely, the value of λ₁ is determined from the reconstruction of the implemented Hamiltonian, and hence does not rely on previous knowledge of Hamiltonian parameters. Thus, measuring λ₁(τ) simultaneously achieves the goal of (i) characterizing a DQS in terms of the implemented stroboscopic Hamiltonian, and (ii) detecting the regime in which Trotter errors on few-body observables are controlled in τ.

Second, we point out that while the PR and level spacing ratio studied in section 2 are standard ways of probing the notion of quantum chaos, their measurement in the context of general quantum many-body systems is non-trivial. In this context, we would like to mention recent advances in this direction [70] based on the use of randomized measurements [71–73]. Conversely, although the behavior of λ₁ is only an indicator of the onset of chaos, it is measurable with costs that in general scale polynomially in the system size.

Finally, we mention that even if the results presented here do not contain any simulated experimental imperfections nor measurement shot noise, the presented protocol retains its validity as long as the dynamics of the system can be treated as approximately unitary for the relevant timescales. Furthermore, via the value of λ₁, our method provides means of assessing the quality of this approximation. A detailed analysis of the influence of measurement noise and imperfections is beyond the scope of this work, and we refer the reader to reference [24] for a detailed study of realistic applications of this protocol in many-body DQS.

While the emphasis of the present work has been on connecting the framework of HL with the regular-chaotic dynamics of the kicked top in the more general context of DQS of many-body systems, we believe that these concepts and techniques can also be used to shed new light on quantum chaos in the kicked top per se, a field that has seen Fritz Haake as one of its most active and brilliant contributors.

We thank Lata K Joshi, Markus Heyl, Philipp Hauke, Nathan Langford and Cahit Kargi for valuable discussions. This work was supported by Simons Collaboration on Ultra-Quantum Matter, which is a Grant from the Simons Foundation (651440, PZ), the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 817482 (Pasquans), and by LASCEM by AFOSR No. 64896-PH-QC. LMS acknowledges support from the Austrian Science Fund (FWF): P 33741-N. The computational results presented here have been obtained (in part) using the LEO HPC infrastructure of the University of Innsbruck.

The data sets generated and analyzed during the current study are available from the corresponding author on reasonable request.

The time-dependent Hamiltonian in equation (1) describes the prototypical and—with H_x being purely linear and H_z purely quadratic in spin operators—simplest version of the kicked top capable of global chaos [1]. Simplicity brings about a high degree of symmetry. This is best discussed in terms of the Floquet operator equation (2): time-reversal covariance of U_τ , which is expressed through the relation $T{U}_{\tau }{T}^{{\dagger}}={U}_{\tau }^{{\dagger}}$ where $T={\mathrm{e}}^{\mathrm{i}{H}_{x}\tau }K$ and K is complex conjugation, implies that U_τ belongs to the COE. Further, U_τ is symmetric under rotations by π around the x axis, ${R}_{x}{U}_{\tau }{R}_{x}^{{\dagger}}={U}_{\tau }$ with ${R}_{x}={\mathrm{e}}^{\mathrm{i}\pi {S}_{x}}$. Hence, an RMT analysis, e.g., of level spacings, has to be carried out independently within the subspaces of eigenstates of U_τ that are even and odd under R_x , respectively. Finally, there are isolated points of even higher symmetry: if h_x τ is a multiple of 2π, then the entire spectrum of H_x τ = h_x τS_x contains only multiplies of 2π, such that ${\mathrm{e}}^{-\mathrm{i}{H}_{x}\tau }=1$ and the Floquet operator reduces to ${U}_{\tau }={\mathrm{e}}^{-\mathrm{i}{H}_{z}\tau }$. Consequently, S_z is conserved, and the dynamics becomes trivially integrable. These points of resonant driving are necessarily encountered when the driving period τ is varied over a wide range, as we do when we study the Trotter threshold, and, while being well understood, they somewhat spoil the phenomenology of the threshold behavior. Further, in view of the fact that the Trotter threshold occurs in a variety of models [20], the existence of such resonant driving points may be regarded as an artefact of the choice of H_x having an equidistant spectrum. Therefore, in reference [2] and in the main text, we focus on the more generic situation, in which H_x is augmented by a term that is quadratic in S_x . Further, by adding a term that is linear in S_z to H_z , we rid the Floquet operator of geometric symmetries, and thus simplify the RMT analysis. Finally, the RMT analysis of HL presented in appendix C, while proceeding along the same lines, simplifies considerably for the CUE as compared to the COE. Hence, to not overburden the presentation of our results with unnecessary technicalities, we choose to focus on a version of the kicked top without time reversal symmetry, which is given in equation (5) [1].

In this appendix, we provide additional details on the HL protocol explained in section 3.1, showing how the reconstruction of the Hamiltonian parameters is achieved via singular value decomposition (SVD). We write the Ansatz H( c ) for the Hamiltonian to be reconstructed as $H(\boldsymbol{c})={\sum }_{j=1}^{{N}_{\mathcal{A}}}\;{c}_{j}\,{h}_{j}$, specified by a chosen set $\mathcal{A}$ of Ansatz operators $\mathcal{A}={\left\{{h}_{j}\right\}}_{j=1}^{{N}_{\mathcal{A}}}$, and we want to find the optimal coefficients c ^rec by solving

$$\begin{equation}{\boldsymbol{c}}^{\text{rec}}=\mathrm{arg}\,\underset{\boldsymbol{c},{\Vert}\boldsymbol{c}{\Vert}=1}{\mathrm{min}}\sum\limits _{i=1}^{{N}_{\text{con}}}{\left\vert \langle H{(\boldsymbol{c})\rangle }_{i,t}-\langle H{(\boldsymbol{c})\rangle }_{i,0}\right\vert }^{2},\end{equation} \tag{ B.1 }$$

with ⟨H( c )⟩_i,t = ⟨ψ_i (t)|H( c )|ψ_i (t)⟩, for a set of N_con chosen initial states |ψ_i (0)⟩. This is equivalent to

$$\begin{equation}{\boldsymbol{c}}^{\text{rec}}=\mathrm{arg}\,\underset{\boldsymbol{c},{\Vert}\boldsymbol{c}{\Vert}=1}{\mathrm{min}}{\Vert}M\boldsymbol{c}{{\Vert}}^{2},\end{equation} \tag{ B.2 }$$

with the constraint matrix M having elements defined in equation (18), and amounts to determining c ^rec as the right singular vector of M corresponding to the smallest singular value λ₁, that is

$$\begin{equation}M{\boldsymbol{c}}^{\text{rec}}={\lambda }_{1}{\boldsymbol{w}}_{1},\end{equation} \tag{ B.3 }$$

where w ₁ is the corresponding left singular vector with || w ₁|| = || c ^rec||. This can be seen as follows: from the SVD of M one writes M = WSV^†, with W having ${N}_{\mathcal{A}}$ orthonormal columns w _n (each being an N_con-component vector), V^† having ${N}_{\mathcal{A}}$ orthonormal rows ${\boldsymbol{v}}_{n}^{{\dagger}}$ (each being an ${N}_{\mathcal{A}}$-component vector), and $S=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\lambda }_{1},\dots ,{\lambda }_{{\text{N}}_{\mathcal{A}}})$, with $0\leqslant {\lambda }_{1}\leqslant \cdots \leqslant {\lambda }_{{N}_{\mathcal{A}}}$. Thus M c = WSV^† c = ∑_n λ_n ( v _n ⋅ c ) w _n , with ⋅ denoting the scalar product, which yields

$$\begin{equation}{\Vert}M\boldsymbol{c}{{\Vert}}^{2}=\sum\limits _{n}{\lambda }_{n}^{2}\,{({\boldsymbol{v}}_{n}\cdot \boldsymbol{c})}^{2}\geqslant {\lambda }_{1}^{2}{({\boldsymbol{v}}_{1}\cdot {\boldsymbol{v}}_{1})}^{2}={\lambda }_{1}^{2},\end{equation} \tag{ B.4 }$$

where the absolute minimum is obtained for c = v ₁. Thus, c ^rec = v ₁.

In this appendix, we show that for U_τ = U ∈ CUE one can make analytical predictions for the form of the matrix $Q=\frac{1}{{N}_{\text{con}}}{M}^{\top }M$, thereby efficiently estimating the value of the λ₁. The result is in good agreement with the value obtained from our HL protocol for large τ and large number of driving cycles. We start by calculating the expectation value of the elements ${Q}_{j,k}=\frac{1}{{N}_{\text{con}}}{\sum }_{i}\;{M}_{i,j}{M}_{i,k}$ over the CUE as

$$\begin{equation}{Q}_{j,k}^{\text{CUE}}=\frac{1}{{N}_{\text{con}}}\sum\limits _{i}{\mathbb{E}}_{U\in \mathrm{C}\mathrm{U}\mathrm{E}}\left[{M}_{i,j}{M}_{i,k}\right]\equiv \frac{1}{{N}_{\text{con}}}\sum\limits _{i}{q}_{j,k}^{\text{CUE}}[{\rho }_{i}],\end{equation} \tag{ C.1 }$$

with ρ_i denoting the chosen initial states for the protocol (which in general can be also mixed). The result depends on the chosen set {ρ_i }, and can be calculated analytically in the case considered in the main text, namely the ρ_i being coherent states. Expressing the matrix elements M_i,j as ${M}_{i,j}=\mathrm{t}\mathrm{r}({\rho }_{i}{h}_{j})-\mathrm{t}\mathrm{r}({\rho }_{i}{U}^{{\dagger}}{h}_{j}U)\left.\right)$, and using the methods developed in [74], we obtain

$$\begin{align}\hfill {q}_{j,k}^{\text{CUE}}[{\rho }_{i}]=\;& \mathrm{t}\mathrm{r}({\rho }_{i}{h}_{j})\mathrm{t}\mathrm{r}({\rho }_{i}{h}_{k})-\frac{\mathrm{t}\mathrm{r}({h}_{j})}{\mathcal{D}}\mathrm{t}\mathrm{r}({\rho }_{i}{h}_{k})\hfill \\ \hfill & -\frac{\mathrm{t}\mathrm{r}({h}_{k})}{\mathcal{D}}\mathrm{t}\mathrm{r}({\rho }_{i}{h}_{j})+\frac{\mathrm{t}\mathrm{r}({h}_{j})\mathrm{t}\mathrm{r}({h}_{k})}{{\mathcal{D}}^{2}-1}+\frac{\mathrm{t}\mathrm{r}({\rho }_{i}^{2})\mathrm{t}\mathrm{r}({h}_{j}{h}_{k})}{{\mathcal{D}}^{2}-1}\hfill \\ \hfill & -\frac{\mathrm{t}\mathrm{r}({h}_{j}{h}_{k})}{\mathcal{D}({\mathcal{D}}^{2}-1)}-\frac{\mathrm{t}\mathrm{r}({\rho }_{i}^{2})\mathrm{t}\mathrm{r}({h}_{j})\mathrm{t}\mathrm{r}({h}_{k})}{\mathcal{D}({\mathcal{D}}^{2}-1)},\hfill \end{align} \tag{ C.2 }$$

where $\mathcal{D}$ is the Hilbert space dimension, i.e., $\mathcal{D}=2S+1$. Since we use pure initial states, $\mathrm{t}\mathrm{r}({\rho }_{i}^{2})=1$. In the case of the kicked top, the Ansatz operators h_j are products of spin operators, and we can use the results of reference [75] to obtain the analytical expressions for the traces tr(h_j ) and tr(h_j h_k ). Where μ, ν = x, y, z. The expectation values tr(ρ_i h_j ) are easily calculated for the chosen initial states. If we are interested in the zeroth order terms in the Floquet Hamiltonian (7), we can use the following identities: ⟨θ, ϕ| S |θ, ϕ⟩ = S(sin θ  cos ϕ, sin θ  sin ϕ, cos θ), $\langle \theta ,\phi \vert {S}_{x}^{2}\vert \theta ,\phi \rangle =S\left(S-\frac{1}{2}\right){(\mathrm{sin}\,\theta \,\mathrm{cos}\,\phi )}^{2}+\frac{S}{2}$, and $\langle \theta ,\phi \vert {S}_{z}^{2}\vert \theta ,\phi \rangle =\frac{S}{4}(1+2S+(2S-1)\mathrm{cos}\,2\theta )$ with S = (S_x , S_y , S_z ). For a large number of initial states, the sum ${N}_{\text{con}}^{-1}{\sum }_{i}$ can be approximated by an expectation value over the initial states ensemble (ISE), i.e.,

$$\begin{equation}\frac{1}{{N}_{\text{con}}}\sum\limits _{i}{q}_{j,k}^{\text{CUE}}[{\rho }_{i}]\approx {\mathbb{E}}_{\rho \in \mathrm{I}\mathrm{S}\mathrm{E}}\,{q}_{j,k}^{\text{CUE}}[\rho ]\equiv {\mathcal{Q}}_{j,k}.\end{equation} \tag{ C.3 }$$

which, for coherent states (ISE = CS), is calculated using ${\mathbb{E}}_{\rho \in \mathrm{C}\mathrm{S}}\left[\mathrm{t}\mathrm{r}(\rho {h}_{j})\right]=\frac{1}{4\pi }\int \mathrm{d}{{\Omega}}_{\theta ,\phi }{\langle {h}_{j}\rangle }_{\theta ,\phi }$ and ${\mathbb{E}}_{\rho \in \mathrm{C}\mathrm{S}}\left[\mathrm{t}\mathrm{r}(\rho {h}_{j})\mathrm{t}\mathrm{r}(\rho {h}_{k})\right]=\frac{1}{4\pi }\int \mathrm{d}{{\Omega}}_{\theta ,\phi }{\langle {h}_{j}\rangle }_{\theta ,\phi }{\langle {h}_{k}\rangle }_{\theta ,\phi }$ with ${\langle {h}_{j}\rangle }_{\theta ,\phi }=\langle \theta ,\phi \vert {h}_{j}\vert \theta ,\phi \rangle$ and $\int \mathrm{d}{{\Omega}}_{\theta ,\phi }={\int }_{0}^{2\pi }\mathrm{d}\phi {\int }_{0}^{\pi }\mathrm{d}\theta \,\mathrm{sin}\,\theta$. Using equations (C.1), (C.2) and (C.3), we can obtain an analytical expression for the matrix $\mathcal{Q}$ which, in the case of Ansatz set ${\mathcal{A}}_{0}=\left\{{S}_{z}^{2},{S}_{z},{S}_{y},{S}_{x}^{2},{S}_{x}\right\}$, reads as (we show only the non-zero matrix elements)

$$\begin{align*}\hfill & {\mathcal{Q}}_{{S}_{z}^{2},{S}_{z}^{2}}={\mathcal{Q}}_{{S}_{x}^{2},{S}_{x}^{2}}=\frac{S(8{S}^{3}-4{S}^{2}+6S-3)}{90},\hfill \\ \hfill & {\mathcal{Q}}_{{S}_{z}^{2},{S}_{x}^{2}}=-\frac{S(8{S}^{3}-4{S}^{2}+6S-3)}{180},\hfill \\ \hfill & {\mathcal{Q}}_{{S}_{x},{S}_{x}}={\mathcal{Q}}_{{S}_{y},{S}_{y}}={\mathcal{Q}}_{{S}_{z},{S}_{z}}=\frac{S(2S+1)}{6},\hfill \end{align*}$$

from which we estimate the RMT value of λ₁, denoted with λ_RMT, as

$$\begin{equation}{\lambda }_{\text{RMT}}=\sqrt{{\varepsilon }_{1}(\mathcal{Q})},\end{equation} \tag{ C.4 }$$

where ${\varepsilon }_{1}(\mathcal{Q})$ denotes the lowest eigenvalue of $\mathcal{Q}$. This yields the estimates used in the main text.

As a final remark we emphasize that, while here we are estimating λ₁ as ${\lambda }_{\text{RMT}}=\sqrt{{\varepsilon }_{1}({\mathbb{E}}_{\text{ISE},\mathrm{C}\mathrm{U}\mathrm{E}}\,Q)}$, the correct prediction would be ${\lambda }_{\text{RMT}}={\mathbb{E}}_{\text{ISE},\mathrm{C}\mathrm{U}\mathrm{E}}\sqrt{{\varepsilon }_{1}(Q)}$. However (i) it is difficult to make analytical progress in the second case and (ii) the two calculations agree very well for large Hilbert space dimension $\mathcal{D}$ and large number N_ini of initial states used. Indeed, the variance of the Q_j,k over the CUE, ${\mathrm{Var}}_{\text{CUE}}{Q}_{j,k}\equiv {N}_{\text{con}}^{-1}{\sum }_{i}\;{\mathrm{Var}}_{\text{CUE}}\left[{M}_{i,j}{M}_{i,k}\right]$, which can also be analytically calculated using the methods in [74], is suppressed in $\mathcal{D}$ as $\mathcal{O}({\mathcal{D}}^{-3})$. Furthermore, the variance of Q_j,k over the ISE, estimated as ${\mathrm{Var}}_{\text{ISE}}{Q}_{j,k}\equiv {\mathrm{Var}}_{\rho \in \mathrm{I}\mathrm{S}\mathrm{E}}\,{q}_{j,k}^{\text{CUE}}[\rho ]$, results in random fluctuations of each realization of Q around the mean value $\mathcal{Q}$, which are suppressed as ${N}_{\text{ini}}^{-1/2}$. That is to say, each realization of Q for a given random unitary and a given set of N_con random initial states can be written as $Q=\mathcal{Q}+{R}_{\text{V}}$ where R_V is a random matrix whose elements can be approximated as

$$\begin{equation*}{({R}_{\text{V}})}_{j,k}\approx {({R}_{\text{V}}(\xi ,\zeta ))}_{j,k}=\xi \sqrt{\frac{{\mathrm{Var}}_{\text{ISE}}{Q}_{j,k}}{{N}_{\text{con}}}}+\zeta \sqrt{{\mathrm{Var}}_{\text{CUE}}{Q}_{j,k}},\end{equation*}$$

with $\xi ,\zeta \sim \mathcal{N}(0,1)$, if we treat the corrections coming from the statistical fluctuations over the CUE and ISE as independent and normal distributed random variables. Thus, R_V can be seen as a small perturbation to Q which, in general, has only perturbative effects on its eigenvalues. In principle, one could achieve a more accurate estimate of λ₁ by calculating ${\lambda }_{\text{RMT}}={\mathbb{E}}_{\xi ,\zeta \sim \mathcal{N}(0,1)}\,\sqrt{{\varepsilon }_{1}\left[\mathcal{Q}+{R}_{\text{V}}(\xi ,\zeta )\right]}$, which can also be done efficiently.