Reliability of analog quantum simulation

Analog quantum simulators (AQS) will likely be the first nontrivial application of quantum technology for predictive simulation. However, there remain questions regarding the degree of confidence that can be placed in the results of AQS since they do not naturally incorporate error correction. Specifically, how do we know whether an analog simulation of a quantum model will produce predictions that agree with the ideal model in the presence of inevitable imperfections? At the same time, there is a widely held expectation that certain quantum simulation questions will be robust to errors and perturbations in the underlying hardware. Resolving these two points of view is a critical step in making the most of this promising technology. In this work we formalize the notion of AQS reliability by determining sensitivity of AQS outputs to underlying parameters, and formulate conditions for robust simulation. Our approach naturally reveals the importance of model symmetries in dictating the robust properties. To demonstrate the approach, we characterize the robust features of a variety of quantum many-body models.

Analog quantum simulators (AQS) will likely be the first nontrivial application of quantum technology for predictive simulation. However, there remain questions regarding the degree of confidence that can be placed in the results of AQS since they do not naturally incorporate error correction. Specifically, how do we know whether an analog simulation of a quantum model will produce predictions that agree with the ideal model in the presence of inevitable imperfections? At the same time there is a widely held expectation that certain quantum simulation questions will be robust to errors and perturbations in the underlying hardware. Resolving these two points of view is a critical step in making the most of this promising technology. In this work we formalize the notion of AQS reliability by determining sensitivity of AQS outputs to underlying parameters, and formulate conditions for robust simulation. Our approach naturally reveals the importance of model symmetries in dictating the robust properties. To demonstrate the approach, we characterize the robust features of a variety of quantum many-body models.
Quantum simulation is an idea that has been at the center of quantum information science since its inception, beginning with Feynman's vision of simulating physics using quantum computers [1]. A quantum simulator is a tunable, engineered device that maintains quantum coherence among its degrees of freedom over long enough timescales to extract information that is not efficiently computable using classical computers. The modern view of quantum simulation differentiates between digital and analog quantum simulations. Specifically, the former performs simulation of a quantum model by using discretized evolutions (i.e., gates) [2][3][4] whereas the latter uses a physical mimic of the model to infer its properties [5]. A crucial issue is that while quantum error correction can be naturally incorporated into digital quantum simulation, this does not seem to be possible for AQS, which are essentially special-purpose hardware platforms built to model systems of interest. However, digital quantum simulators are extremely challenging to build, whereas AQS are more feasible in the near future, with several experimental candidates already under study [6][7][8][9][10]. Thus a critical question for the quantum simulation field is: as AQS become more sophisticated and begin to model systems that are not classically simulable, can one verify or certify the accuracy of results from systems that are inevitably affected by noises and experimental imperfections? [11].
In response to this challenge, we develop a technique for analyzing the robustness of an AQS to experimental imperfections. We specialize to AQS that prepare ground or thermal states of quantum many-body models since these are the most common types of AQS currently under experimental development. * Electronic address: mnsarov@sandia.gov † Electronic address: zhangjun12@sjtu.edu.cn

I. DEFINITIONS
Define a quantum simulation model, notated (H, O), as consisting of a Hamiltonian H and an observable of interest O (both Hermitian operators). We write a general Hamiltonian in parameterized form as H(λ) = K k=1 λ k H k , where λ = (λ 1 , . . . , λ K ) T denotes the vector of parameters ( = 1 throughout this paper). H k are the terms in the Hamiltonian that are individually tunable through the parameters λ k . In addition, we decompose the observable into orthogonal projectors representing individual measurement outcomes O = M m=1 θ m P m with P m P n = P m δ mn 1 .
The goal of an AQS is to produce the probability distribution of a measurement of O under a thermal state or ground state of a system governed by H(λ 0 ), where λ 0 denotes the ideal, nominal values of the system parameters. That is, to produce the distribution p m (λ 0 ) = tr(P m (λ 0 )), m = 1, · · · , M , where (λ 0 ) = e −βH(λ 0 ) /tr e −βH(λ 0 ) , for some inverse temperature β = 1/k B T , if the goal is to predict thermal properties of the model; or (λ 0 ) = ψ g (λ 0 ) ψ g (λ 0 ) with ψ g (λ 0 ) being the ground state of H(λ 0 ), if the goal is to predict ground state properties. However, due to inevitable environmental interactions, miscalibration, or control errors, the parameters λ k can deviate from their nominal values, which can potentially corrupt AQS predictions. We quantify the reliability of an AQS by the robustness of this probability distribution with respect to the deviations of λ from its ideal value λ 0 .
In general, there is no reason to expect that the prepared state (λ) will be robust to perturbations of λ. In fact, we know that for Hamiltonians that possess a quantum critical point, thermal and ground states can be extremely sensitive to λ around that point [12][13][14]. However, reliable AQS does not require robustness of (λ) around λ 0 , but only robustness of the probability distribution of observable outcomes, {p m } M m=1 . The fact that this is a less demanding requirement is the fundamental reason to expect that some models may be reliably simulated using AQS.

II. QUANTIFYING AQS ROBUSTNESS
To quantify the reliability or robustness of an AQS, we begin by utilizing the Kullback-Leibler (KL) divergence to measure the difference between the measurement probability distributions p(λ) and p(λ 0 ) [15]: . Assuming that the deviation in parameters from the ideal, ∆λ = λ − λ 0 , is small, we expand the KL divergence to second order to obtain The positive semidefinite matrix F is the Fisher information matrix (FIM) for the model, whose elements are given by [15]: In Appendices A & B we describe how to compute the FIM for a quantum simulation model in closed-form, without using numerical approximations to derivatives. Note that even though we adopt the KL divergence to motivate the FIM,Cencov's theorem states that the FIM is the unique Riemannian metric for the space of probability distributions under some mild conditions [16], and is therefore a general measure of the sensitivity of the parameterized outcome distribution around λ 0 . We first note that if the parameter deviations, ∆λ, are Gaussian distributed with zero mean then the expected KL-divergence can be approximated to second-order by the trace of the FIM. This follows from Eq. (1), and the fact that 1 i Az i is an estimate of the trace of A when the elements of z i are independent, standard normal variables [17]. However, we are interested in not only obtaining such an average measure of AQS robustness, but also in understanding the factors that determine robustness, or lack thereof, of a particular model. For this purpose we turn to a spectral analysis of the FIM associated with a quantum simulation model. Consider the set of eigenvalues ζ k and eigenvectors v k of F , with k indexing the eigenvalues in descending order. Since F is a symmetric matrix, we have F = K k=1 ζ k v k v † k . Then the simulation error caused by the deviated parameter λ can be approximated to the second order by This error is influenced by two quantities: the magnitude of the eigenvalues, and the overlap of the eigenvectors with the parameter deviation. We can use this structure to quantify the robustness of AQS outputs to the system parameter deviations around the ideal λ 0 . A quantum simulation model is trivially robust to parameter deviations if all ζ k ≈ 0; i.e., F ≈ 0. In the high temperature limit, β → 0, we can show that F (λ 0 ) → 0 at the rate of β 2 generically and so all models become trivially robust Appendix E. This is expected since the equilibrium state becomes dominated by thermal fluctuations at high temperatures, and observables become insensitive to underlying Hamiltonian parameters.
A more interesting way a model can be robust is if the FIM possesses only a small number of dominant eigenvalues that are separated by orders of magnitude from other eigenvalues. In this case, only parameter deviations in the directions given by the eigenvectors of dominant eigenvalues affect the simulation results. For instance, if ζ 1 is the dominant eigenvalue, then the composite parameter deviation (CPD) v † 1 ∆λ has the major influence on simulation errors. We refer to AQS models that have FIMs with a few dominant eigenvalues separated by orders of magnitude from the rest as sloppy models. This terminology is adopted from statistical physics, where it has been recently established that a wide variety of physical models possess properties that are extremely insensitive to a majority of underlying model parameters , a phenomenon termed parameter space compression (PSC) [18,19].
Model sloppiness is a prerequisite for non-trivial AQS robustness, since without this property an AQS can only be robust if most or all Hamiltonian parameters can be precisely controlled, an impractical task as quantum simulation models scale in size. In contrast, given a sloppy quantum simulation model, one only has to control and stabilize a few ( K) influential CPDs. However, model sloppiness alone is not sufficient for AQS robustness since the practicality of controlling these influential CPDs has to be evaluated within the context of the particular AQS experiment at hand, including its control limitations and error model. In this work we aim for a general analysis and do not focus on any particular AQS implementation. Instead, we demonstrate that many quantum simulation models exhibit model sloppiness, the prerequisite for robustness, and how this can help to identify the parameters that must be controlled in order to produce reliable AQS predictions.

III. ANALYZING THE FIM
A low rank FIM immediately indicates a sloppy model, and since the rank is an analytically accessible quantity, we can use the FIM rank to study model sloppiness beyond numerical simulations. In particular, in this section we discuss two useful methods for bounding the rank of the FIM for a quantum simulation model.
We begin by rewriting the FIM in a compact form. Define a matrix V ∈ R K×M , whose km-th entry is ∂pm(λ) ∂λ k , and Λ = diag{p 1 (λ), p 2 (λ), · · · , p M (λ)}. Then the FIM can be written as F = V Λ −1 V † . Here we assume that all p m (λ) are non-zero. In the case when some p m (λ) equal 0, these elements and the corresponding rows in V should be removed.
This factorized form of the FIM immediately provides a useful bound on its rank. Notice that the row sum of V is zero, therefore the rank of V is at most M − 1, which is an upper bound on the rank of F . In many physical situations, it is common that the number of distinct measurement outcomes is much less than the number of model parameters, i.e., M K. In this case, the rank bound of M − 1 can immediately signal a sloppy model. An example of this that we shall encounter later is a spin-spin correlation function observable, whence M = 2 and K typically scales with n, the number of spins in the model.
Next we will show that fundamental symmetries of the quantum simluation model can reduce the rank of the FIM, and further, that symmetries can be used to deduce the structure of the FIM eigenvectors and characterize the influential CPDs. To do this, we define the symmetry group of a quantum simulation model, G, as the largest set of symmetries shared by the Hamiltonian and the observable in the model -i.e., the maximal group of space transformations that leave the Hamiltonian and the observable invariant. Let {U g } g∈G be a faithful unitary representation of this symmetry group for the quantum simulation model 2 , and suppose U g H k U † g = H j for some k, j, g. Then in Appendix C we show that ∂pm(λ) ∂λ k = ∂pm(λ) ∂λj for all m, under ground or thermal states. Therefore, the spatial symmetry of the model leads to identical rows in V , and we see an immediate connection between model symmetry and model sloppiness: a high degree of symmetry yields a significant redundancy in the FIM and only a few non-zero eigenvalues.
This observation suggests a constructive procedure to formulate an upper bound on the rank of FIM based on model symmetries. Specifically, compute the orbit of H k under the symmetry group for the quantum simulation model; i.e., {U g H k U † g |g ∈ G}, for all 1 ≤ k ≤ K. The number of orbits will be the maximum number of distinct rows in the matrix V , and therefore provides an upper bound to the rank of the FIM.
The repeated rows in V resulting from model symmetries also informs us about the structure of the eigen-  vectors of the FIM, and as a result, the structure of the influential CPDs. Explicitly, the CPD takes the form (see where s indexes the unique orbits, and µ k s is a scalar dependent on the orbit, nominal parameter values and temperature. Although the forms of the CPDs are always determined by the eigenvectors of F and therefore by the symmetries of the model, i.e., Eq. (3), the coefficients µ k s (λ 0 , β) are temperature-dependent and the structure of the CPD can simplify further if these coefficients become alike or approach zero as temperature changes. We will encounter instances of this in the next section.

IV. APPLICATIONS
In this section we use the rank bounds derived above and numerical simulations to understand the sloppiness and robustness of several quantum simulation models. In addition to the applications presented here, we analyze several other quantum simulation models in Appendix G.

A. 1D transverse-field Ising model
The well-known transverse field Ising model in one dimension (1D-TFIM) is described by the Hamiltonian: where σ i α is a Pauli operators acting on spin i with α = x, y, or z, and is normalized such that {σ α , σ β } = δ αβ I 2 .
We are interested in the uniform version of this model with B 0 i = B 0 and J 0 i = J 0 for all i; however, when this model is simulated by an AQS, the actual values of B i and J i may fluctuate around these nominal values. The boundary conditions for this model can be either periodic, i.e., σ n+1 x ≡ σ 1 x , in which case the Hamiltonian will be denoted as H per 1 ; or open, i.e., J n = 0, in which case the Hamiltonian will be denoted as H open 1 . Although this model is efficiently solvable [20][21][22], its role as a paradigmatic quantum many-body model with a nontrivial phase diagram makes it a useful benchmark for quantum simulation. Moreover, it exhibits many generic phenomena related to robust AQS, as we will show below.
Two observables of interest in this model are the net transverse magnetization S z = n i=1 σ i z and two-point correlation functions C z (i, j) = σ i z σ j z . It is feasible to measure these observables experimentally, and importantly, they probe the magnetic order in the system. For example, both of these observables can be used to characterize a quantum phase transition that occurs in the ground state of the uniform 1D-TFIM when swept past its quantum critical point at J 0 /2B 0 = 1 [23].
First we consider the quantum model {H per 1 , S z } with fixed J 0 , and sweep the parameter B 0 to explore the behavior of the model across its phase diagram. This quantum simulation model has full translational invariance. The orbit of any σ i z under the (lattice) translation group contains all σ j z , 1 ≤ j ≤ n, and the orbit of any for all m and 1 ≤ i, j ≤ n; that is, all the rows in V corresponding to B and J are identical, respectively. Hence, an upper bound on the rank for the FIM of this model is 2, for all possible J 0 , B 0 , β, and n. This is a very sloppy model, especially for large n.
To illustrate this general result, in Fig. 1 we show the eigenvalues of the FIM for a 10-spin 1D-TFIM with J 0 = 1, as B 0 is swept. The rank bound derived above is evident in this figure -there are two dominant eigenvalues -and the negligible eigenvalues shown in Fig. 1 (gray lines) are actually numerical artifacts. In fact, the largest eigenvalue is also orders of magnitude above the second largest, except in the region of the quantum critical point, where the second eigenvalue approaches it (although still many orders of magnitude smaller).
The eigenvectors associated to the two dominant eigenvalues prescribe the parameter deviations that the model is most sensitive to, and due to the full translational invariance of the model we find that they exhibit particularly simple structure (regardless of β). Namely, the two dominant eigenvectors take the form [µ, · · · µ, η, · · · , η] T and [−η, · · · , −η, µ, · · · , µ] T , where µ and η are two scalars depending on the value of B 0 . This implies that across all phases, the model is sensitive only to the CPDs i ∆B i and i ∆J i . Hence, this quantum simulation model will be robust to parameters deviations as long as these two sums are maintained at zero; i.e., local fluctuations of the microscopic parameters that (spatially) average to zero are inconsequential.
Next we examine the AQS model {H per 1 , C z (i, j)}i.e., the 1D-TFIM with periodic boundary and a correlation function observable. Noticing that the observable has only two outcomes immediately indicates that the rank of F is at most one, and hence this model is also very sloppy, especially for large n. To illustrate this in Fig. 2(a) we show eigenvalues of the FIM for a 10-spin example, with the observable being the correlation function C z (2, 6), for zero and intermediate temperature. As expected, only one eigenvalue is significant and all the others are zero up to numerical precision across the whole phase diagram (values of J 0 /2B 0 ).
The structure of the dominant eigenvector is more complex in this case, since although the Hamiltonian is translationally invariant, the observable is not. The eigenvector structure can be extracted from symmetry considerations, but for simplicity we plot its components for the n = 10 case in Fig. 2(b), (c), for β = ∞, β = 1, respectively. Focusing on the zero temperature case first ( Fig. 2(b)), we see that the CPD takes the form Unlike the previous quantum simulation model {H per 1 , S z }, the form of the linear combination of underlying model parameters that the AQS is sensitive to not only depends on B 0 , but this dependence is not the same for all 20 parameters. Another interesting aspect of Fig. 2(b) is that away from the quantum critical point, the composite parameter is mostly composed of model parameter variations near the spins whose correlation is being evaluated. More specifically, the AQS model is most sensitive to (∆B 2 + ∆B 6 ) + (∆B 1 + ∆B 3 + ∆B 5 + ∆B 7 )/2 and (∆J 1 + ∆J 2 + ∆J 5 + ∆J 6 ) (i.e., the parameters local to spins involved in the correlation function C z (2, 6)). However, near the quantum critical point, all underlying parameter changes enter into the definition of the influential CPD. This is a novel manifestation of collective phenomena in quantum many-body systems: whereas local correlations are typically influenced by local parameters, near a critical point, local correlations are influenced by all the parameters in the system.
The complexity of the influential CPD for this model is most evident when the system is in its ground state 3 , but these features persist for small finite temperatures also. However, as shown in Fig. 2(c), the structure of the CPD simplifies with increased simulation temperature. The sensitivity to all parameter variations in the model around the region near the quantum critical point dis- (b) Influential CPD when system is in ground state appears at intermediate temperature, as expected, since thermal fluctuations overwhelm signatures of quantum criticality as the temperature increases [24]. Moreover, the influential CPD becomes composed of only the parameter changes at the spins involved in the correlation function (∆B 2 +∆B 6 and ∆J 1 +∆J 2 +∆J 5 +∆J 6 ) across the whole phase diagram.
We pause to reflect on the differences between the two models examined so far. Whereas {H per 1 , S z } and {H per 1 , C z (i, j)} are both sloppy quantum simulation models, the influential CPD for the former is much simpler in form -its form remains invariant across the phase diagram and with varying temperature. An immediate consequence is that if the goal of a quantum simulation of the 1D-TFIM is to characterize the phase diagram and the phase transition, one should utilize the transverse magnetization as an experimental observable as opposed to correlation functions since the former is more robust to independent local parameter fluctuations. Another option is to probe the site averaged correlation function (C z (j) = 1 n i σ i z σ i+j z ) in which case the translational invariance, and consequently robustness to independent local parameter fluctuations of the quantum simulation model is restored.
To study a model with a lower degree of symmetry, we now turn to the 1D-TFIM with open boundary conditions, with the observable of interest being transverse magnetization again; i.e., the quantum simulation model {H open 1 , S z }. This model is no longer translationally invariant, but has reflection symmetry about the center spin (for odd n) or center coupling (for even n). Under this symmetry, each orbit contains at most two elements -e.g., the orbit of σ j z contains itself and σ n+1−j z -and hence an upper bound on the rank of the (2n−1)×(2n−1) matrix F is n. In this case symmetry considerations do not completely reveal the sloppiness of the model, that is, the FIM rank bound is weak, as n is not a lot less than 2n − 1. We explicitly calculate the FIM for this model with n = 10 at low temperature, and Fig. 3(a) shows its eigenvalues as a function of B 0 . As expected from the symmetry rank bound, the model has at most n = 10 eigenvalues that are nonzero (within numerical precision). Furthermore, the first eigenvalue is several orders of magnitude larger than the others at all phases, although there is a pronounced aggregation of eigenvalues around the quantum critical point. Hence the model is sloppy although not to the same degree as the previous two models examined. The influential CPDs for this model takes the form: where µ i and η i are B 0 -dependent real numbers. Therefore this model is robust to parameter fluctuations that are negatively correlated across its center spin (or coupling for even n). As a result of the complexity of these CPDs and the overall lower degree of sloppiness, we conclude that an AQS implementation of this model will be less robust to parameter fluctuations than the previous two 1D-TFIM models considered.

B. 2D transverse field Ising model
Now we study the uniform 2D-TFIM on an n×n square lattice: with net magnetization S z as the observable of interest.  and J ij = J 0 . In this case the model has two types of planar symmetries: rotational symmetry about the center of the lattice and mirror reflection symmetry about four reflection lines. The net magnetization observable is invariant under the above symmetries. This is not an exactly solvable model as in the 1D-FTIM case and is therefore of more fundamental interest for AQS.
Several local terms (σ i z ) and coupling terms (σ i x σ j x ) in the Hamiltonian are mapped to the same orbit under the action of the symmetry transformations for {H 2 , S z }. For example, Fig. 4 shows the lattice sites and couplings that lie in the same orbit for a 3 × 3 lattice. There are a total of five distinct orbits in this case and thus the rank the 19 × 19 FIM is upper bounded by five. Also, according to Eq. (3) fluctuations of the local magnetic fields or spinspin couplings that act on identically colored site or edges in Fig. 4 will be grouped together in the influential CPD. Explicit computations of eigenvalues and CPDs for this model are included in Appendix G.1.

C. Fermi-Hubbard model
The Fermi-Hubbard Hamiltonian, a minimal model of interacting electrons in materials, is of significant interest to the AQS community since it is thought that understanding emergent properties of this model could explain some high-T c superconducting materials [25]. The Hamiltonian takes the form: where c † iσ (c iσ ) creates (annihilates) an electron with spin σ ∈ {↑, ↓} on site i, n iσ = c † iσ c iσ is the electron number operator for site i. We consider this Hamiltonian defined over a two-dimensional lattice, and the i, j indicates that the first sum runs over nearest neighbor sites. Moreover, t ij represents the coupling energy between sites that induces hopping of electrons, and U i > 0 represents the repulsive energy between two electrons on the same site. We are interested in the uniform version of this Hamiltonian with nominal parameters U i = U 0 , for all i and t ij = t 0 , for all i, j. The observable of interest is the double occupancy fraction, D = 2 n i n i↑ n i↓ , where n is the total number of sites, which for example can be used to probe metal to insulator transitions in this model.
In Fig. 5 we show FIM properties for this AQS on a 2 × 3 lattice with periodic boundary conditions. We show results from simulations of the Hubbard model at half-filling ( i n i↑ = i n i↓ = 3), but the results are qualitatively the same for the slightly doped cases as well. Fig. 5(a) shows sites and coupling energies that lie within the same orbit under symmetry transformations for this model, which are lattice translations in the x and y directions. All Hamiltonian terms that act locally are mapped between each other and all couplings are mapped between each other, and thus there are three distinct orbits for this model implying an upper bound on the rank of the FIM of 3. Fig. 5(b) shows eigenvalues of the model with t 0 = 1, as a function of U 0 . As expected, there are always at most three non-zero eigenvalues (to numerical precision) and the model is extremely sloppy. In contrast to the models examined so far, the low temperature version of this model is sloppier than the intermediate temperature version. Finally, Fig. 5(c) confirms that the influential composite parameter deviations take the form expected from the symmetry analysis, with the model only showing sensitivity to the sum of local fluctuations i ∆U i , and sum of vertical coupling terms or sum of horizontal coupling terms.

V. SCALING TO LARGE SYSTEMS
Quantum simulation is most compelling for large-scale quantum models since difficulty of classical simulation typically increases exponentially with the model scale 4 . Obviously, evaluation of model robustness through classical computation of the FIM is not possible for large-scale models. However, we will show how analysis of smallscale systems can be bootstrapped by various techniques to draw useful conclusions about their large-scale versions.
First, we note that the bounds on the rank of the FIM that we derived earlier can be useful for models of any scale. For example, the rank bound derived from symmetry considerations allows us to determine the sloppiness of the quantum simulation model {H per 1 , S z } at any scale (i.e., for any number of spins); and further, symmetry considerations yield the form of the CPD that the model is sensitive to. More generally, we observe that the FIM for any quantum simulation model is greatly simplified by translational invariance, and this can be used to determine sloppiness of the model at any scale.

Consider a general (finite-dimensional) translationally invariant Hamiltonian
is an operator acting on degrees of freedom in the spatial neighborhood N , and of type α. As an example, consider the following general spin-1/2 Hamiltonian on a 3D lattice with nearest-neighbor interactions and periodic boundary conditions in all directions: where i, j indicates the sum runs over nearest neighbors in all three directions. Here α ∈ {x, y, z, xx, yy, zz} and the neighborhoods are local sites or edges of the 3D lattice. Translational invariance implies that under the action of the translation symmetry group for these models, all Hamiltonian terms of a given type α lie in the same orbit. Therefore, the number of orbits is the same as the number of types of interaction, and assuming that the observable of interest is also translationally invariant, A is an upper bound on the rank of the FIM for such models at any scale. Thus such models are guaranteed to be sloppy, except at very small scales (where the number of parameters is comparable to A). Furthermore, the AQS will be most susceptible to the CPDs ∆λ α N for each α. For example, for the spin-1/2 Hamiltonian H 4 above, if the observable is also translationally invariant, e.g., S x , S y or S z , then the FIM for this quantum simulation model will have rank at most 6, for any number FIG. 6: Influential CPD for the model {H per 1 , Cz(2, 10)} evaluated with n = 70 spins, when the system is in ground state. This model has 140 microscopic parameters, only the ones that significantly contribute to the influential CPD are labeled for clarity.
of spins. Note that this example covers a wide range of models including tilted and transverse field Ising models and a variety of Heisenberg models.
The rank bound obtained by counting the number of observable outcomes is also useful in determining sloppiness at any scale. For example, the spin-1/2 correlation C α (i, j) = σ i α σ j α has only two possible outcomes ±1, thus the FIM rank is always one, regardless of the Hamiltonian and number of spins. Unfortunately, this bound does not also inform us about the structure of the CPD that the model is sensitive to.
Second, even in cases where a complete symmetry analysis is not possible, an analysis of the small-scale model can be informative about the robustness of the corresponding large-scale model. In particular, since the form of the CPDs is determined by symmetries of the model, one can extrapolate from the form of the CPDs from small-scale models to large versions. For example, for the model {H per 1 , C z (i, j)} studied above, we can examine large-scale behavior by using the well-known exact solution to the 1D-TFIM [20,21] (see Appendix F for details), and confirm that the form of the influential CPD remains the same at large n as for the small-scale version. In Fig. 6 we plot entries of the dominant eigenvector for the model {H per 1 , C z (2, 10)} for n = 70 spins in the ground state. The influential CPD is mostly composed of parameters around the spins whose correlation function is being evaluated, except near the quantum critical point when other parameters also contribute. These trends agree with results for the small-scale version of the model shown in Fig. 2(b).
Third, we note that in some cases we can approximate a quantum simulation model with one of higher symme- try in order to gain more information from the FIM. An example of such an approximation is the common practice of imposing periodic boundary conditions on finite lattices in order to make calculations tractable. This approximation can also be useful for assessing robustness of large-scale models using our approach. To illustrate this, we turn to the exact solution of the 1D-TFIM again, and confirm that the model {H open 1 , S z } can be approximated by {H per 1 , S z } as the number of spins increases. Our numerical investigations show that when n is large, e.g., n > 50, the largest eigenvalue of the FIMs for these two models become almost identical, and the forms of the influential CPDs for the two models approach each other. Hence for some large-scale models one can infer sloppiness and robustness from analysis of approximations with higher degree of symmetry. Of course such approximations are not always possible and one should be aware of their accuracy across parameter regimes.
Finally, we pose a conjecture regarding the behavior of sloppiness with scale: if a small-scale AQS model with a lattice quantum many-body Hamiltonian is sloppy, then its large-scale version will also be sloppy. Although we currently lack a proof of this statement, it is well supported by numerical evidence. For example, consider the model {H open 1 , S z } that was shown to be sloppy at small scales earlier. By utilizing the exact solution to the 1D-TFIM, we can analytically calculate the FIM for a large number of spins. We choose B 0 = 0.45, J 0 = 1, and β = 1, and in Fig. 7 plot the largest 10 eigenvalues of the FIM for this model as a function of the number of spins, n. The model remains sloppy across all scales that were simulated.

VI. DISCUSSION
We have developed and applied a formalism for analyzing the robustness of analog quantum simulators. Many quantum many-body models are potentially robust for AQS, especially if they possess a high degree of symmetry, which we have shown leads to model sloppiness, a necessary condition for robustness. In addition, our techniques allow one to determine which underlying parameter(s) impact simulation results the most, which could help to focus experimental effort when designing AQS platforms. In a sense, our work can be thought of providing a formal justification of the commonly encountered intuition that bulk properties should be immune to microscopic fluctuations, and elucidating the connection between this intuition and system symmetries.
For brevity we have only presented results from applying our approach to uniform models above. However, we have analyzed a large variety of more general models, including ones with random parameters and long-range couplings, and some of the results from these studies are presented in Appendix G. Application of our approach to these more complex cases with less symmetry illustrates how any symmetries in the underlying ideal model can be exploited to understand sloppiness and robustness. While nearly all the quantum simulation models we studied were sloppy (the exception being models with complete disorder, i.e., random parameters), in some cases the influential CPD is complex, and engineering robust AQS for these models could be challenging. This finding is mirrored by the ubiquity of sloppiness in the classical models studied by Sethna et al. [18,19].
The intent of this work is to introduce the notion of sloppy models to AQS, demonstrate its relation to robust simulation and illustrate that certain quantum simulation models can be robust to uncertainties in parameters. There are many promising directions to extend this work. For example, while we have focused on AQS that prepare ground or thermal states of quantum many-body models, the approach can be extended to analyze quantum simulations that predict dynamic properties of quantum models by considering probability distributions for the dynamical variables of interest. Finally, we have restricted ourselves in this work to investigating the robustness of analog simulation of Hamiltonian models with calibration uncertainties because these uncertainties can in fact dominate the behavior of existing cold-atom analog quantum simulation platforms, e.g., [7][8][9][10], where decoherence due to environmental coupling is very small. However, for a complete picture of robustness, it is desirable to extend this analysis to diagnose robustness of quantum simulation models with decoherence.

Appendix A: Calculation of FIM for thermal states
We can analytically simplify the partial derivatives required to compute the FIM when the system is in a thermal state (λ) = e −βH(λ) /Z, where Z = tr e −βH(λ) . Now we have (A1) In order to calculate ∂e −βH(λ) ∂λ k , we utilize Eq. (78) in Ref. [26] to obtain: Note that we drop the λ-dependence when it is clear from the context. Now we diagonalize the Hamiltonian as where T is a unitary matrix of eigenvectors and Γ = diag{γ 1 , γ 2 , · · · } is a diagonal matrix of eigenvalues. Substituting this decomposition into Eq. (A2), we get where denotes the Hadamard product, i.e., , element-wise product, and Θ pq (τ ) = e (γq−γp)τ is the pq-th element of Θ. The τ dependence is entirely in this matrix, and therefore we can evaluate this integral to yield: where Φ is a matrix with elements: Consequently, Inserting these expressions into Eq. (E1) allows us to evaluate the derivatives required to calculate the FIM for thermal states in a manner that is numerically stable.
Appendix B: Calculation of FIM for ground states.
The FIM when the system is in its ground state, |ψ g , can also be obtained in an analytical manner. We must calculate where gs ≡ |ψ g ψ g |. For a Hamiltonian with a simple (non-degenerate) minimum eigenvalue, the minimum eigenvalue and the associated eigenvector are infinitely differentiable in a neighborhood of H, and their differentials at H(λ) are [27] dE = ψ g | (dH) |ψ g (B2) and d |ψ g = (E 0 I n − H(λ)) + (dH) |ψ g , where + denotes the Moore-Penrose (MP) pseudoinverse. We then obtain and therefore, V , the matrix of partial derivatives can then be written in a compact matrix form as: These analytical expressions for the derivatives for thermal and ground states are faster and more numerically stable to evaluate than approximations using difference equations.
Appendix C: FIM and model symmetries.
In the main text, we stated that if a quantum simulation model has a symmetry transformation that relates H k and H j , then This has consequences for the rank of the FIM for the model. To prove the above, we start with the explicit expressions for the partial derivatives under thermal states, given in Eq. (E1). The two k dependent quantities in this expression can be written, using Eq. (A2) as: Then suppose the quantum simulation possesses a symmetry with unitary representation (we assume the symmetry group is compact) = 0 for all g. Furthermore, given the decomposition of the observable, [U g , P m ] = 0, ∀g, m. Now, suppose the symmetry maps H j to H k , meaning H k = U g H j U † g , then using the commutation properties stated above, Also, tr P m ∂e −βH ∂λ k = − tr P m e −βH/2 Therefore, all k-dependent terms in Eq. (E1) are the same if we exchange k with j, and hence we arrive at Eq. (C1) for thermal states. To prove the same property when the system is in its ground state, we turn to the expression for the partial derivatives given in Eq. (B5): Since [U g , H(λ)] = 0, and both of these operators are normal, they share an eigenbasis, implying [U g , gs ] = 0. Therefore, Using [U g , H(λ)] = 0, it is easy to verify that U g (E 0 − H(λ)) + U † g is also the MP pseudoinverse of E 0 I − H(λ), and from the uniqueness of MP pseudoinverse, we have that From this equality and Eq. (C3), Eq. (C1) follows for ground states as well.
Appendix D: Structure of the eigenvectors of F As discussed in the main text, spatial symmetries of a quantum simulation model render some rows of the matrix V equal. Here we show that this induces a certain structure on the Fisher information matrix (FIM), namely that the corresponding entries of each eigenvector of F are equal.
Without loss of generality, we assume that V can be written as where 1 k is a column vector with dimension n k and all entries being 1, and v T k are pairwise distinct row vectors. As a result, and p T = p 1 · · · p s is an eigenvector of M D with eigenvalue α. Then Therefore, p 1 1 T 1 p 2 1 T 2 · · · p s 1 T s T is an eigenvector of F . From Eq. (D1), we know that the rank of V is s, and thus the ranks of M and F are both s. Hence, all the eigenvectors of F can be written in the form p 1 1 T 1 p 2 1 T 2 · · · p s 1 T s T , that is, they have the same structure of repeated entries as V in Eq. (D1).
We will show that in the limit of high temperature, the FIM approaches 0 at the rate of β 2 . For simplicity, we consider an n-qubit system. From Appendix A, and therefore we know that when the system is in a thermal state (λ) = e −βH(λ) /Z, we have where Z = tr e −βH(λ) . In the high temperature limit, β → 0, we expand to the first order Further, using this approximation and ignoring higher order terms in β, we get and tr(P m e −βH ) tr ∂e −βH ∂λ k Z 2 ≈ −β tr P m (I − βH) tr H k (2 −2n + 2 −3n+1 β tr H) Combining these two equations, we have where Define a matrix U whose km-th element is u km . Then F = β 2 U Λ −1 U † . Hence, as β → 0, the FIM approaches the zero matrix as β 2 and thus the quantum simulation is robust. Furthermore, U Λ −1 U † is a constant matrix that is independent of the system parameters, which indicates that at high temperature the quantum simulation is completely insensitive to the nominal values of the underlying parameters.

Net magnetization distribution for the 1D-TFIM
Recall that the Hamiltonian for the 1D-TFIM is given by Consider the observable S z = n j=1 σ j z = m θ m P m , where in the second equality we have decomposed the observable as a sum of projectors. We wish to compute p m = tr(P m ), and we use a two-step procedure to calculate this quantity. First, we express each P m as a linear combination of {S 1 , · · · , S n }: where . . .
Second, we calculate the expection values of S j , i.e., S j = tr(S j ). Combining these two steps, we have We now elaborate on the details of these two steps. First, we express P m in terms of S j . The observable P m can be written as where |κ j is a state with m − 1 spins in the ground state |0 and n − m + 1 spins in the excited state |1 , and N m = n m−1 . For simplicity, we use the case m = 2 to illustrate the approach. In this case, we have P 2 =|01 · · · 1 01 · · · 1| + |101 · · · 1 101 · · · 1| + · · · + |1 · · · 10 1 · · · 10| Since |0 0| = I/2 + σ z and |1 1| = I/2 − σ z , we have To find the coefficients ξ mj , we replace I ⊗n /2 by 1 2 and σ j z by a scalar variable x j in Eq. (F8) and obtain the following polynomial: The polynomial p 2 is symmetric and thus can be represented by elementary symmetric polynomials s j : x k1 , x k1 x k2 , The coefficients to represent P 2 in terms of S j are identical to those that represent p 2 in terms of s j , that is, In fact, to obtain ξ mj , we can choose all the variables x j to be the same x. Then, we have n m Equating the coefficients in both sides of Eq. (F12), we can obtain ξ mj . Next we show how to compute S j . From Refs. [20,21], we define two matrices P and Q as Let φ T k be a normalized row eigenvector of (P −Q)(P +Q), i.e., φ T k (P −Q)(P +Q) = Λ 2 . Juxtapose φ T k and ψ T k into two matrices Φ and Ψ. For the calculation of ground state, we define and for the thermal state, we let From Wick's theorem and Ref. [21], we know that S j is the sum of all the j-by-j principle minor of G. Moreover, from Ref. [28], we have det(tI − G) = t n − S 1 t n−1 + S 2 t n−2 − · · · ± S n .
Hence we can determine S j by calculating the characteristic polynomial of G. With these two steps, we can now obtain p m .

Correlation function distribution for the 1D-TFIM
When the observable is the correlation function C z (i, j) = σ i z σ j z , we know from Eq. (2.33c) in Ref. [21] that under the ground state, and under the thermal state, where G g and G t are defined in Eqs. (F14) and (F15), respectively. We then consider to analytically calculate the FIM for ground state. Since σ i z σ j z has two eigenvalues ± 1 4 , we obtain that for ground state, Then ij dλ l and dp 2 dλ l = − dp 1 dλ l . (F20) We now derive dG g /dλ l . Since G g = Ψ T Φ, we have The matrix (P − Q)(P + Q) is simple, meaning that it has pairwise distinct eigenvalues. Then its eigenvalue and the associated eigenvector are infinitely differentiable in a neighborhood of H(λ) and their differentials are where + denotes the Moore-Penrose pseudoinverse. From the definition of P and Q in Eq. (F13), it is straightforward to derive dP/dλ l and dQ/dλ l and thus Moreover, we have that where Combining these equations, we can calculate dp 1 /dλ l and dp 2 /dλ l for ground state analytically. For thermal states, we just need to calculate an additional derivative of tanh( β 2 Λ) in G t and can obtain the results similarly. When the observables are σ i x σ j x and σ i y σ j y , their mean values can be obtained from Eq. (2.33a) and (2.33b) in Ref. [21]. And following similar procedures as above, we can derive analytical expressions for derivatives of the measurement probabilities.   In this section we report the behavior of the FIM for some quantum simulation models that were not included in the main text for conciseness.

2D transverse field Ising model
In the main text we demonstrate how symmetry analysis of the 2D-TFIM with open boundary conditions and net magnetization as the observable enables one to determine the rank of FIM for this model, and show that it is sloppy. For more details on the symmetry analysis for this model, see section H in this Appendix. Here in Fig. 8, we explicitly present the eigenvalues and eigenvectors of the FIM for a 3 × 3 square lattice version of this model. It is evident from Fig. 8(b) that the FIM eigenvalues agree with the rank bound (rank ≤ 5) derived from symmetry. Furthermore, Figs. 8(c) and (d) show that the forms of the influential CPDs respect the symmetry of the model. A bound derived from considering the number of measurement outcomes tells us that the rank of the FIM is at most 10, and therefore we show the ten largest eigenvalues in color and the others (numerical artifacts) in gray. (b) The form of the first influential CPD.

1D random Ising model
To examine a model with disorder, consider the 1D transverse field Ising model with random local fields and coupling energies, i.e., with periodic boundary conditions (σ n+1 x ≡ σ 1 x ), and B 0 i = B 0 +δB i , J 0 i = J 0 +δJ i , where δB i and δJ i are independent zero-mean Gaussian random variables with standard deviation σ. As for the observable of interest, consider the net magnetization S z again. This quantum simulation model has no symmetries due to the random parameters and so the FIM rank bounds based on symmetry are not informative. The number of measurement outcomes for this observable is M = n + 1, and therefore the rank of the FIM is at most n. In Fig. 9(a) we show the eigenvalues of the FIM for a 10-spin example of this quantum simulation model, with J 0 = 1, disorder variance σ = 0.2 and β = 10. This figure shows the FIM eigenvalues for one representative sample of δB i and δJ i . As evident from this figure, while the dominant eigenvalue is roughly two orders of magnitude above all others, this model cannot be considered sloppy except for small or large values of B 0 . In Fig. 9(b) we also show the form of the first influential CPD (we do not label the points on this plot since we only wish to illustrate the complexity of the behavior of this quantity for this model).

J1-J2 antiferromagnetic Heisenberg model
Now we turn to a quantum simulation model based on a Hamiltonian that contains non-nearest-neighbor interactions and geometric frustration. The J 1 -J 2 antiferromagnetic Heisenberg model is defined by the following Hamiltonian governing spin-1/2 systems on a two-dimensional lattice: where the first sum is over nearest-neighbor spins and the second is over next-nearest-neighbor spins. We are interested in the uniform nominal operating point for this model where J 0 ij = J 0 and K 0 ij = K 0 with J 0 , K 0 > 0 5 . Fig. 10 shows a single plaquette in the square lattice in the nominal model.  The magnetic order in this system is complex with different phases of magnetic ordering being driven by competition between the two different kinds of interactions. The magnetic order parameter is different in different K 0 /J 0 regimes. For small values of this ratio (∼ 0) the magnetization is Néel ordered (the model resembles a conventional Heisenberg antiferromagnet on a square lattice in this regime), and as this ratio approached unity one has so-called "striped magnetization" [29]. Our observables of interest is the staggered magnetization, which probes the Néel order in the system: where n is the total number of spins in the system. The quantum simulation model {H 5 , M s } with open boundary conditions on the lattice has several symmetries despite the complicated form of the observable of interest. For square lattices, this model has rotational symmetry about the center of the lattice and reflection symmetry about four reflection lines. In Fig. 11(a) we explicitly show the symmetries in this model for a 3 × 3 square lattice. Note that since n is odd, all these symmetry transformations take odd (even) labeled spins to odd (even) labelled spins, and hence leave the observable of interest invariant. From this symmetry analysis, we obtain a rank bound on the FIM of rank(F ) ≤ 4. Fig. 11(b) shows the eigenvalues of the FIM for this 3 × 3 example for β = 10 and β = 1, and it is clear that the rank bound is respected. Finally, Fig.  11(c) shows the primary influential CPD for this model when β = 10. The first four eigenvectors of the FIM all define influential CPDs since the first four eigenvalues are non-negligible. We only plot the primary influential CPD here for simplicity, but all the others have the same symmetry properties.