Zero and Finite Temperature Quantum Simulations Powered by Quantum Magic

Andi Gu¹, Hong-Ye Hu¹, Di Luo^1,2,3, Taylor L. Patti⁴, Nicholas C. Rubin⁵, and Susanne F. Yelin^1,3

We introduce a quantum information theory-inspired method to improve the characterization of many-body Hamiltonians on near-term quantum devices. We design a new class of similarity transformations that, when applied as a preprocessing step, can substantially simplify a Hamiltonian for subsequent analysis on quantum hardware. By design, these transformations can be identified and applied efficiently using purely classical resources. In practice, these transformations allow us to shorten requisite physical circuit-depths, overcoming constraints imposed by imperfect near-term hardware. Importantly, the quality of our transformations is tunable: we define a ‘ladder’ of transformations that yields increasingly simple Hamiltonians at the cost of more classical computation. Using quantum chemistry as a benchmark application, we demonstrate that our protocol leads to significant performance improvements for zero and finite temperature free energy calculations on both digital and analog quantum hardware. Specifically, our energy estimates not only outperform traditional Hartree-Fock solutions, but this performance gap also consistently widens as we tune up the quality of our transformations. In short, our quantum information-based approach opens promising new pathways to realizing useful and feasible quantum chemistry algorithms on near-term hardware.

A central task in quantum chemistry is determining ground state energies and finite temperature free energies for electronic Hamiltonians. While numerous algorithms aim to leverage quantum hardware for solutions [21, 31, 47, 48], the constraints of near-term hardware, particularly limited circuit depth, pose challenges. One approach to address this difficulty is classically assisted quantum algorithms, wherein computation is offloaded to a classical computer where possible, and the quantum hardware only supplies the quantum resources which are absolutely necessary to identify an accurate solution. A natural way in which a classical computer might assist a quantum algorithm is to preprocess the electronic Hamiltonian into a more favorable form before it is analyzed on quantum hardware. For instance, one can use a classical algorithm to partition the terms of a Hamiltonian $H$ into commuting groups which can be measured simultaneously, so that calculating the expectation $\langle H \rangle$ can be done with fewer samples [42].

In this work, we consider another preprocessing method: similarity transformations, wherein our aim is to identify a unitary change of basis $U$ such that the transformed Hamiltonian $H_{\mathrm{eff}} \equiv U^{\dagger}HU$ is more amenable to analysis on near-term quantum hardware. Indeed, similarity transforms have long been used to simplify Hamiltonians [4, 8, 13, 22, 24, 33, 37, 44, 46]. In this work, we take a novel quantum information theoretic approach to similarity transforms: we center our techniques around two fundamental measures of ‘quantumness’ – namely, entanglement and nonstabilizerness (popularly termed ‘quantum magic’). These notions of quantumness formalize the intuition that the degree to which a state is truly quantum can be measured by how difficult it is to simulate classically. For instance, low entanglement states can be simulated by matrix product states [10]. However, entanglement alone is insufficient to completely characterize quantumness – classically simulable (i.e., using classical resources that scales polynomially in system size) circuits [1], known as Clifford circuits can generate an unbounded amount of entanglement. Indeed, it has been shown that random Clifford circuits will almost always produce volume-law entangled states [43]. Yet, these circuits are unable to express a different quantum resource: quantum magic. More precisely, Clifford circuits are generated by a non-universal gateset (composed of the Hadamard, phase, and CNOT gates); this has motivated the formal definition of magic as a measure of the amount of non-Clifford resources required to prepare a particular state [18, 35]. Like entanglement, quantum magic has also been shown to be directly connected with the hardness of classical simulation [5].

Our main insight in this work is that similarity transforms can be used to eliminate spurious quantumness (i.e., classically simulable entanglement and magic) in the Hamiltonian’s eigenbasis before it is analyzed on quantum hardware. Then, the real hardware only needs to supply the fundamentally necessary quantumness for the problem, hence helping shorten requisite circuit depths. We introduce a class of similarity transformations, represented as Clifford circuits augmented with a few (non-Clifford) single-qubit $Z$-rotation gates. For brevity, we call these circuits Clifford + kRz circuits. Conveniently, this design allows good similarity transforms to be both identified and applied efficiently using purely classical resources. The Clifford component of the circuit can be interpreted as eliminating spurious entanglement from $H_{\text{eff}}$, while the non-Clifford gates eliminate spurious magic. Importantly, by varying the number of non-Clifford gates, we show that the power of these transformations is tunable, so that we can identify increasingly simple eigenbases (i.e., exhibiting less entanglement and magic) at the cost of increased classical computational resources. Since this tunability comes in discrete steps, we term this the “quantum magic ladder” (see Fig. 1). By applying our approach to zero and finite temperature free energy calculations, we demonstrate that this reduction in quantumness can have significant practical benefits. Specifically, we show that our similarity transforms find energies which are significantly more accurate than other popular classical solutions, and moreover, the performance gap only improves as we climb the quantum magic ladder. Furthermore, we present a classical algorithm which outputs a new (simplified) effective Hamiltonian $H_{\mathrm{eff}}$, related to the original Hamiltonian via a similarity transform. We show that the resultant simplified Hamiltonians $H_{\mathrm{eff}}$ produced by our method enables significantly reduced circuit depths on quantum hardware. Under this quantum information theoretic lens, these results demonstrate a promising pathway for designing practical quantum chemistry algorithms in the NISQ era.

Quantum Chemistry. Quantum chemistry is a natural domain for assessing the efficacy of our protocol. The Hamiltonians produced by quantum chemistry systems have immediate physical interpretations, allowing us to apply physical intuition to better understand how our protocol modulates the ‘quantumness’ of a problem. Furthermore, there is a wealth of established quantum chemistry methods against which we can benchmark our protocol, such as the Hartree-Fock (HF) approximation. In short, quantum chemistry allows us to demonstrate the advantages of our approach in a well-understood physical context.

Figure 1: A general approach for similarity transformations using the quantum magic ladder. Starting with a second quantized Hamiltonian $H$, we use classical simulable quantum circuits (e.g., Clifford + kRz circuits) to reduce the quantumness of the Hamiltonian by the basis transformation $H_{\mathrm{eff}} = U^{\dagger}HU$. We climb up the quantum magic ladder by increasing $k$, which increases the degree to which $U$ can simplify the Hamiltonian (by reducing the quantum magic of the $H_{\text{eff}}$ eigenbasis). Note that for finite temperature calculations, when working in the grand canonical ensemble, we replace $H$ with $H - \mu N$.

In quantum chemistry, one key goal of the electronic structure problem is to solve the energy of many-body electronic Hamiltonians at zero and finite temperature. These electronic Hamiltonians are often expressed in a molecular orbital basis:

where $p,q,r,s$ label the spin-orbitals and $\hat{c}_{\alpha}, \; \hat{c}_{\alpha}^{\dagger}$ are the annihilation and creation fermionic operators. Anticipating analysis on qubit-based quantum computers, we apply a fermion-qubit mapping (e.g., parity mapping [36], Jordan-Wigner [20], Bravyi-Kitaev [41]) to map the fermionic Hamiltonian to a qubitized Hamiltonian $H = \sum_{i} \alpha_i P_i$, where $P_i$ is an $n$-qubit Pauli operator. For the molecules considered in this work (LiH and $\text{H}_2\text{O}$), the qubitized Hamiltonians contain $O(10^2)$ Pauli strings. The explicit construction of these Hamiltonians was performed using the OpenFermion library [28]. In this work, we focus on two paradigmatic problems involving these electronic Hamiltonians, which are to calculate the ground state energy of $H$, and the free energy of the Gibbs state at finite temperature [9]. In this work, we work in the grand canonical ensemble, so that the Gibbs state is $\rho \propto \exp(-\beta(H - \mu N))$, with inverse temperature $\beta$, chemical potential $\mu$, and particle number $N$.

The Quantum Magic Ladder. Similarity transformations have historically occupied a key role in the development of quantum mechanics. An early example is the Schrieffer-Wolff method, which generates low-energy effective Hamiltonians by constructing similarity transformations that separate high-energy and low-energy subspaces [4, 8, 33, 44]. In keeping with this tradition, we argue that similarity transforms can also play a central role in bringing quantum chemistry algorithms closer to practical feasibility in the NISQ era. One way of understanding the role of a similarity transformation $U$, particularly in the context of near-term hardware, is as a virtual circuit that is appended to the end of a ground state preparation circuit $V$. That is, if we define a transformed Hamiltonian $H_{\text{eff}} = U^{\dagger}HU$, we observe that $\langle V^{\dagger} H_{\text{eff}} V \rangle = \langle (UV)^{\dagger} H(UV) \rangle$, so that, in effect, the applied circuit is $UV$. This is especially important for NISQ devices, which face limited circuit depths (hence limited expressive power of $V$). An appropriate similarity transformation $U$ can act like an appended circuit that significantly eases the burden of the real quantum circuit $V$.

In what follows, we describe a unified framework to programatically construct a new class of similarity transformations for quantum chemistry algorithms, a system we term “the quantum magic ladder.” These transformations are comprised of Clifford circuits augmented by a few Rz gates. There are several reasons for the attractiveness of the Clifford + kRz approach. Firstly, Clifford circuits are classically simulable with complexity $O(n^3)$, where $n$ is the system size [1]. As we discuss below, these Clifford circuits can be dramatically enhanced by doping them with just a few ($k$) non-Clifford gates. While this enhancement comes at the cost of an increase in simulation complexity that is exponential in $k$, Clifford + kRz transforms with $k \leq O(\log n)$ can still be simulated with polynomial classical resources [1, 3, 6]. Secondly our circuits strictly generalize more conventional circuits such as Clifford circuits with a few T gates. Indeed, we find that the additional continuous degrees of freedom introduced by Clifford + kRz circuits (i.e., the rotation angle of each Rz gate) is both easy to optimize and yields significant performance improvements. Furthermore, these improvements are free in the sense that the additional simulation cost for kRz gates is exactly the same as that of $k$ T gates. The simulation of these Clifford + kRz circuits was performed using the PyClifford package [19] (for details, see Appendix A). Finally, as illustrated in Fig. 1, ‘climbing’ the quantum magic ladder provides a systematic way to increase the expressive power of our similarity transformation: one simply needs to increase $k$.

The enhanced expressivity of ansatze higher on this ladder can be understood through the lens of quantum magic. Much like entanglement defines one notion of quantumness, quantum magic, or nonstabilizerness, is another metric which measures the degree to which a state can be considered truly quantum [18, 35]. Our main insight in this work is that the optimal Clifford + kRz transform can model the quantumness of $H$ along both axes of quantumness, capturing its entanglement via the Clifford component of the circuit, and its quantum magic via the non-Clifford components. As we climb the magic ladder, we increase the degree to which our transform can model the quantum magic of the energy eigenbasis (hence reducing the magic of the transformed Hamiltonian’s eigenbasis). In the limit as $k \to \infty$, our ansatz can model any unitary [29], which ultimately leads to a full reduction in the Hamiltonian’s quantumness (i.e., complete diagonalization).

Our Clifford + kRz circuits are parameterized as $U(\Theta)$, where $\Theta$ is a combination of discrete and continuous parameters, corresponding to the Clifford and Rz components respectively. To identify good transformations $U(\Theta)$, we make use of the variational principle for both ground state energies and finite temperature free energies. Specifically, consider the zero temperature case: $U(\Theta)$ defines a parameterized effective Hamiltonian $H_{\mathrm{eff}}(\Theta)$, and our aim is then to minimize $\langle 0|H_{\mathrm{eff}}(\Theta)|0\rangle$ (which is equivalent to minimizing the energy, with respect to $H$, of $U(\Theta)|0\rangle$). Optimizing $\Theta$ is a hybrid discrete-continuous problem with a large search space, so we develop a three-step approach. First, we observe that qubitized electronic Hamiltonians strongly resemble disordered Hamiltonians, whose energy coefficients span a wide range of scales. We leverage this formal similarity by using the spectrum bifurcation renormalization group [13, 37, 46] (which produces Clifford similarity transformations effective for modelling strongly-disordered quantum materials) to find an initialization point $\Theta_0$ for our Clifford + kRz transformation. We then perform a combination of greedy search and gradient descent to further optimize $\Theta$ (see Appendix B for details).

Results. We apply our Clifford+kRz transformation for both ground state and finite temperature quantum chemistry calculations [7, 12, 16, 32, 39] on LiH ($n = 8$) and $\text{H}_2\text{O}$ ($n = 12$). While these system sizes are modest compared to the capabilities of state-of-the-art classical numerical methods, they serve as valuable proof-of-concept demonstrations for our quantum-classical hybrid approach, and are representative of the system scales that are expected to be tractable on near-term quantum devices. We start with the zero temperature case, where we aim to identify a ground state energy $E = \langle 0|H_{\text{eff}}(\Theta)|0\rangle$ that closely approximates the true ground state energy $E_0$. In Fig. 2, the energy gap $E - E_0$ is plotted as a function of different bond lengths between atoms, where $E$ is the energy predicted by an optimized from Clifford + kRz transformation. First, compare the accuracy to the accuracy of the Hartree-Fock (HF) solution: since the HF solution is a strict subset of the Clifford ansatz [38], the energy gap of a pure Clifford transformation $(k = 0)$ is upper bounded by that of the HF solution. However, a pure Clifford circuit only surpasses the HF solution when the bond length exceeds the equilibrium bond length. In contrast, as we climb the quantum magic ladder, the energy gap is significantly reduced for all the bond lengths, including bond lengths shorter than the equilibrium bond length. Moreover, as the number of non-Clifford gates increases, the energy gap consistently decreases. This is clear evidence that tuning up the power of Clifford + kRz ansatz can significantly improve accuracy for ground state calculations. It is important to note that the number of non-Clifford gates $k$ should be considered in relation to the total number of qubits $n$. If we allow $k \, = \, \exp(\Theta(n))$, the doped circuit can approximate an arbitrary unitary gate. In contrast, in our study, we focus on the regime where $k\,\ll\,n$ to maintain the classical simulability of the circuits while still enhancing their expressibility. We observe that the convergence of our method improves as $k$ increases, with significant improvements seen for $k=1$ and $k = 2$. However, the improvement becomes marginal for higher values of $k$. There is a tradeoff between the expressibility gained by increasing $k$ and the computational cost associated with it; it appears that there are sharply diminishing returns beyond $k=2$. Future work could investigate this tradeoff in more detail and develop strategies for optimally selecting $k$ based on the problem at hand. We also note the presence of kinks in some curves of Fig. 2, which are artifacts likely due to numerical precision limitations and minor optimization errors, rather than physical phenomena.

Figure 2: $E-E_0$ is the gap between $E$, the energy calculated by the optimized Clifford + kRz circuit, and $E_0$, the true ground state energy. We illustrate this gap for LiH and $\text{H}_2\text{O}$ as a function of their bond lengths. Note that the HF solution is consistently higher in energy, as it is strictly less expressive than the Clifford ansatz, even when $k = 0$, underscoring the enhanced representational power of the Clifford circuits in capturing the ground state properties. Furthermore, the accuracy (as measured by $E-E_0$) improves as we climb the quantum magic ladder. The inset shows the stabilizer Rényi entropy $M$, which quantifies the magic of the ground and first excited states at different bond lengths.

Interestingly, the energy gap peaks around a particular bond length, indicating a region of maximal quantum magic. To show this, we employ a popular metric specifically designed to quantify the magic of a state $|\psi\rangle$. This metric $M(|\psi\rangle)$ is known as the stabilizer entropy [23], and is defined as

where $\mathbb{P}_n$ is the $n$-qubit Pauli group. Notably, it obeys $M(|\psi\rangle) = 0$ if and only if $|\psi\rangle$ is a stabilizer state, and it is strictly positive otherwise. In general, a higher value of $M(|\psi\rangle)$ indicates a greater degree of nonstabilizerness [23]. We see in the inset of Fig. 2 that quantum magic (as measured by $M$) and the energy gap peak near the same bond length, indicating that regions where the ground state has high quantum magic cannot accurately be reduced by basis transformations with insufficient non-Clifford gates. Interestingly, we observe that the magic of the ground state is typically very small at extreme bond lengths, and peaks at intermediate bond lengths. This peak in magic and energy gap likely corresponds to the Coulson-Fisher point [11], where the electronic structure transitions from single-reference to multi-reference character. This transition point is known to be challenging for computational methods [25], aligning with our observation of increased quantum resources being necessary in this region. Our quantum information theoretic approach thus provides a new perspective on this well-established concept in quantum chemistry. This suggests that in practice, hyperparameters like $k$ ought to be set accordingly — for extreme bond lengths, a small $k$ should suffice, while intermediate bond lengths will require a larger $k$ to accurately capture the increased magic. Near these intermediate bond length regions, accurate results demand either more powerful classical similarity transformations (i.e., larger $k$) or further analysis on real quantum devices.

Figure 3: Pulse time reduction with the transformed Hamiltonian, depicted as a violin plot. This displays the probability density of a solution having a particular energy gap, with the width representing the solution density at each energy level. We model an analog quantum machine with a $^{87}$Rb atomic chain for LiH ground state preparation using fixed lattice spacing and total evolution time. The Rabi frequency and the local detuning of each atom are optimized. For each evolution time, twenty randomly initialized variational ground states are optimized.

We demonstrate this by applying our method on analog quantum computers. Specifically, we classically simulate the dynamics of a one-dimensional Rydberg array, which serves as a programmable quantum simulator (PQS) for preparing a ground state of LiH at a bond length 3.4Å. Note that this bond length exhibits high magic (Fig. 2), indicating that additional quantumness, supplied by the PQS, is necessary to prepare the ground state. To prepare this state, we fix the total PQS evolution time $T$ and evolve the initial state $|0\rangle$ under a time-dependent Rydberg Hamiltonian [2] with variationally optimized site-dependent Rabi frequencies $\Omega_i(t)$ and local detunings $\Delta_i(t)$. For each fixed total evolution time $T$, we variationally optimize 20 randomly initialized ground states for both the bare Hamiltonian $H$, and for the Clifford-transformed Hamiltonian $H_{\rm eff}$, using automatic differentiation with Jax [15] and Tensorcircuit [50] to directly optimize the pulse shapes $\Omega_i(t), \Delta_i(t)$. The results are shown in Fig. 3. Consider just the case of $T = 0.1 \mu s$, for which dramatic differences are already clear. The optimized energy with the original Hamiltonian is concentrated at levels significantly above the HF energy, because the entanglement generated by the PQS within this small timeframe is limited. However, the Clifford + kRz similarity transform already significantly reduces the entanglement of the ground state before evolution on the PQS; hence, even under the limited evolution time $T = 0.1 \mu s$, the PQS can quickly identify a solution that vastly surpasses the accuracy of the Hartree Fock solution. This demonstrates that Clifford + kRz transformations effectively reduce requisite entanglement and pulse times, a significant advantage for near-term quantum devices.

Figure 4: Free energy gap for LiH at a bond length 3.4Å. The inset shows the entanglement negativity $\mathcal{N}$ averaged over each site of the Gibbs state for the transformed Hamiltonian.

For finite temperature calculations, we apply our approach to the LiH molecule for grand free energy estimation. Our ansatz for the density matrix is

where $U(\Theta)$ is a Clifford + kRz circuit. Again applying the variational principle, we solve for $\Theta$ by minimizing the free energy $F = \langle H \rangle - \mu \langle N \rangle - TS \text{ of } \rho(\Theta)$ (see Appendix B for details). The results in Fig. 4 show that the Clifford + kRz ansatz consistently reduces the free energy gap $F - F_0$ (where $F$ is the free energy of the optimized density matrix $\rho(\Theta)$, and $F_0$ is the free energy of the true Gibbs state) as we climb the quantum magic ladder. Similar to the zero temperature case, all Clifford transformations outperform the finite temperature HF method, including the pure Clifford transformation ($k = 0$), which reduces the free energy gap by $\sim 50\%$.

To quantify the extent to which our transforms reduce the quantumness of the transformed Hamiltonian, we examine the entanglement features and off-diagonal terms of the exact Gibbs state for the transformed Hamiltonian $H_{\mathrm{eff}}$. We quantify the entanglement of a mixed state with its entanglement negativity $\mathcal{N}(\rho) = (\|\rho^{\Gamma_A}\|_1 - 1)/2$, where $\rho^{\Gamma_A}$ is the partial transpose of $\rho$ with respect to subsystem $A$ [14, 17, 30]. Using this metric, in the inset of Fig. 4, we show that our similarity transformations indeed consistently reduce the entanglement of the

Figure 5: Density distribution of the off-diagonal term for the Gibbs state $\rho$ of LiH across different temperatures for bond length 3.4Å. The width of features is proportional to solution density.

Gibbs state. Furthermore, Fig. 5 shows that, across a wide temperature range, our transformation’s ability to suppress off-diagonal elements in the exact Gibbs state $\rho$ grows with $k$. This provides further evidence that Clifford + kRz transformations provide an effective method to reduce the quantumness of a Hamiltonian, and moreover, that this reduction is tunable by varying $k$.

Conclusions. We have developed a general approach for zero and finite temperature quantum calculations from a quantum information theoretic perspective, facilitated by the Clifford + kRz transformation. One of the primary benefits of our approach is that it is fully classically simulable for small $k$, providing an convenient method to conserve costly quantum resources. Furthermore, by tuning $k$, we provide a means of customizing the similarity transforms produced by our technique: larger $k$ enables greatly increased expressive power and accuracy, at the cost of increased classical computation. This tunability forms the basis of what we term the quantum magic ladder. We apply these methods in the context of ground state and free energy calculations for quantum chemistry, demonstrating superior performance to HF techniques and consistent improvement as we ascend the quantum magic ladder, for both digital and analog quantum computers. We examine how the Clifford + kRz transformation reduces quantumness (as measured by entanglement and magic) before any algorithms are run on quantum hardware, and demonstrated how this reduction directly translates to reduced requirements for quantum resources. While our approach shows promise for small-scale systems like LiH and $\mathrm{H}_{2}\mathrm{O}$, it is necessary to investigate its scalability to larger systems. The classical simulation cost of the Clifford+kRz circuits grows exponentially with $k$, which could limit the applicability of our method to larger systems. However, it is encouraging to note that we are already able to find significant improvements for $k \leq 2$, which suggests that even a small number of strategically placed non-Clifford gates can provide significant benefits. As quantum hardware continues to improve, we expect that the insights and techniques developed in this work, such as the use of classically-optimized similarity transformations to reduce the “quantumness" of the problem, will be valuable for scaling up to larger system sizes.

Having established this framework, we also identify a number of remaining open questions and future directions. Algorithmically, in addition to the greedy search algorithms used in this work, powerful methods such as simulated annealing and differentiable quantum architecture search [49] could be considered. Our similarity transformation based on doped Clifford circuits could also potentially augment well-known numerical approaches such as density matrix renormalization group (DMRG) [34]. Since our method identifies effective Hamiltonians whose ground states have significantly lower entanglement compared to the original problem, it could potentially aid in the construction of more efficient tensor network representations. Moreover, recent work [27] has shown that the stabilizer formalism can be integrated with tensor networks. By transforming the problem into a basis where the ground state is closer to a stabilizer state, our approach could enable the use of these compact representations. Investigating the interplay between our similarity transformations and tensor network methods could lead to powerful new tools for simulating quantum many-body systems. Finally, to scale finite temperature calculations to larger systems, high-dimensional classical probability distributions of the Gibbs state can be simulated using generative models in the spirit of $\beta\text{-VQE}$ [26]. Our approach can also be applied to the finite temperature calculation based on the thermal double field framework [40, 45]. This offers an opportunity to tailor our transformations to specific applications, opening new pathways for near-term quantum computers to address a wide array of quantum chemistry problems.

Acknowledgments

The authors would like to thank Chong Sun, Yi Tan, Linqing Peng, Yuan Liu, Shi-Xin Zhang, Joonho Lee, James Whitfield, and Ryan Babbush for the helpful discussions. DL acknowledges support from the NSF AI Institute for Artificial Intelligence and Fundamental Interactions (IAIFI). This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (C2QA) under contract number DE-SC0012704. HYH is grateful for the support from Harvard Quantum Initiative Fellowship. SFY would like to thank the NSF through a Cornell HDR and the CUA PFC grant.

References

A Classical simulation of Clifford-dominated quantum circuits

In this section, we review the basic idea of the generalized stabilizer representation and how to simulate non-Clifford gates and general quantum states with stabilizer frames. First, let us define Pauli group and stabilizer group:

Definition A.1 (Pauli group). The Pauli group on $n$-qubits is the set of all tensor products of Pauli matrices, i.e.

Definition A.2 (Stabilizer group). Stabilizer group $S$ is a subgroup of Pauli group satisfies the following properties:

1. $S$ is an abelian group. For $\forall s_i, s_j \in S, \ [s_i, s_j] = 0$.

2. It doesn’t contain negative identity, i.e. $-I \notin S$.

The stabilizer group can always be represented by its generators. For example, a stabilizer group containing $2^r$ elements can be represented by $r$ generators of the group: $S = \langle s_1, s_2, \dots, s_r \rangle$. One should notice the generators are not unique. For example, $S = \langle s_1, s_2, \dots, s_r \rangle = \langle s_1 s_2, s_2, \dots, s_r \rangle$. With the stabilizer group, one can define stabilizer states:

Theorem A.1 (Pure stabilizer state). A full stabilizer ($n$ generators for a $n$-qubit system) stabilizes a unique pure state $|\psi_S\rangle$

One can check the density matrix can be represented as

Another concept that is useful for describing general quantum state is the destabilizers:

Definition A.3 (Destabilizer group). A destabilizer group $\mathcal{D}$ is a subgroup of Pauli group associated with the stabilizer group $S = \langle s_1, s_2, \dots, s_r \rangle$. It has the following properties:

1. $\mathcal{D}$ has the same number of generators as $S$: $\mathcal{D} = \langle d_1, d_2, \dots, d_r \rangle$.

2. For each generator $d_i$, it anti-commute with $s_i$, but commute with all the other generators $s_j, j \neq i$.

Definition A.4 (Full tableau). A full stabilizer tableau $\mathcal{T} = (S, \mathcal{D})$ is a pair of full stabilizer and destabilizer.

Definition A.5 (Stabilizer basis). The stabilizer basis with full tableau $\mathcal{T}$ is $\mathcal{B}(\mathcal{T}) = \{d \left| \psi_S \right\rangle \left| d \in \mathcal{D} \right\}$

The stabilizer basis is an orthonormal and complete basis for $n$-qubit Hilbert space. Since the destabilizers are note unique, if $\mathcal{D}$ and $\mathcal{D}'$ are both destabilizers for $S$, then $\mathcal{B}(S,D)$ and $\mathcal{B}(S,D')$ are the same set of basis modulo phase difference. Since $\mathcal{B}(\mathcal{T})$ is a complete basis for Hilbert space, then any quantum state can be represented by the stabilizer basis. Therefore, we introduce the generalized stabilizer representation.

Definition A.6 (Generalized stabilizer representation). Any pure quantum state can be represented with stabilizer basis. More particularly, the density matrix $\rho$ of any pure quantum state can be written as

Therefore, we call the two-tuple $(\chi,\mathcal{B}(S,\mathcal{D}))$ the generalized stabilizer representation.

For a general state, in the worst case, we need $\mathcal{O}(4^n)$ classical memory to store and simulate its dynamics. However, this formulation is particularly useful when this decomposition is sparse. For simplicity, we can use the zero norm of $\chi$ matrix to denote its complexity, i.e. $\Lambda(\chi) = ||\chi||_0$. And the simulation complexity scales as $\mathcal{O}(\Lambda(\chi))$. Now we define the Pauli channel and see how to simulate general quantum states with the generalized stabilizer representation.

Definition A.7 (Pauli channel). The Pauli channel of a $n$-qubit system can be written as

where $P_i, P_j \in G_n$, and $\phi_{ij}$ are complex numbers. We use $\Lambda(\mathcal{E}) = ||\phi||_0$ to denote the complexity of the Pauli channel.

A.1 General state evolution under Clifford gates

The evolution of the general state under a Clifford circuit $C$ evolution can be written as

Therefore, we only need to update the full tableau $\mathcal{T} = (S, \mathcal{D})$, which can be done with $\mathcal{O}(n^2)$ time. Under Clifford gates, the general stabilizer frame will update as $(\chi, S, \mathcal{D}) \to (\chi, S', \mathcal{D}')$.

A.2 General state evolution under Pauli channels

First, we use the fact that any Pauli string can be decomposed with stabilizer and destabilizers, $P_m = \alpha_m d_{b_m} s_{c_m}$. An efficient implementation of this decomposition is in Algorithm 1. Suppose we want to simulate a general Pauli channel, $\mathcal{E}(\rho) = \sum_{m,n} \phi_{m,n} P_m \rho P_n^{\dagger}$. For each term in the Pauli channel:

Therefore, for each term in Pauli channel, we need to update $(\chi, S, \mathcal{D}) \to (\chi', S, \mathcal{D})$. And we sum over all the terms in the Pauli channel, which contains $\Lambda(\phi)$ terms. Especially, T-gate can be written as the following channel:

where $\Lambda(T)=4$.

B Optimization of Clifford+kRz circuits

We use a brickwork circuit ansatz (see Fig. 6) with $L$ layers. One layer is composed of $n$ Clifford gates, each of which acts on neighboring qubits. Although a common gate set for the Clifford group is taken to be the Hadamard, phase, and CNOT gates, we instead use the gate set composed of $e^{i\pi P/4}$, where $P \in \{I, X, Y, Z\}^{\otimes 2}$. We dope the circuit with $k$ Rz gates at a fixed layer, but a variable site index (which is optimized over). To summarize, our free parameters belong to the domain $\mathbb{Z}_{16}^{n \times L} \times \mathbb{Z}_n^k \times \mathbb{R}^k$. The first term corresponds to 16 possible choices for the Clifford gates (of which there are $n \times L$). The second corresponds to the $n$ possible choices for the site index of the $k$ Rz gates, and the last is simply the argument for the $k$ Rz gates.

Figure 6: The brickwork circuit ansatz for $L=2$ and $k=2$.

To optimize our circuit, we treat the discrete and continuous parameters separately. We label $X \in \mathbb{Z}_{16}^{n \times L} \! \times \! \mathbb{Z}_{n}^{k}$ the discrete parameters and $\Theta\,\in\,\mathbb{R}^k$ the continuous ones. For any cost function $f(X,\Theta)$, we can define a ‘marginalized’ cost function by minimizing over $\Theta\colon$

This is a useful construction when it is easy to do this minimization; it turns out that for small $k$, this process is very easy. Since the cost is periodic in $\Theta$, we only need to search the small volume $\Theta \in [0, 2\pi]^k$. This search is further simplified by the fact that $f(X,\Theta)$ is continuous in $\Theta$, and does not have many local minima (see Fig. 7). By applying a simple gradient-based optimizer (e.g., gradient descent, conjugate gradient, etc.) to just a few ($\sim 10$) randomly initialized choices of $\Theta$, we can often find the global minimum within very few iterations. This is to say, the marginalized cost function $f(X)$ is easy to evaluate. With this, we optimize over the remaining parameters $X$ with a simple greedy search. Each iteration, we pick one of the parameters in $X$ to optimize over. For instance, if we choose to optimize $P_1$ (refer to Fig. 6), we would evaluate the cost with all other parameters held constant, varying $P_1$ over all 16 possibilities $\{I, X, Y, Z\}^{\otimes 2}$. Then, we simply replace $P_1$ with the choice that yielded the lowest cost $f(X)$. Finally, we do this discrete optimization in parallel for a large number ($n_{init} \sim 1000$) of randomly initialized parameters $X$, and simply take the best performing solution after a fixed number ($n_{iter} \sim 100$) iterations. Although we considered more advanced schemes for this discrete optimization, such as a simulated annealing process, we found that these schemes did not yield better results.

Figure 7: The landscape of $f(\Theta) \equiv f(X_{best}, \Theta)$ for LiH at a bond length 3.4Å, where $X_{best}$ was the best performing solution for $k=2$. On the left, we show the landscape for the ground state calculation, and on the right, we show the landscape at a finite temperature $\beta = 200$ ($T \approx 1500$K). Note that there are very few suboptimal local minima – applying gradient descent to a randomly initialized $\Theta$ will find the global minimum with high probability.

C Symmetry of the Recovered Solutions

A noteworthy drawback of the Clifford + kRz ansatz (Fig. 6) is that it does not necessarily conserve any symmetries of the Hamiltonian. Despite this, we observe that the optimization process prefers to recover solutions that are both eigenstates of the number and spin operators, and furthermore, typically states that conserve symmetries in the Hamiltonian (e.g., total spin zero). Although observe that the latter does not hold for every solution (see Fig. 9), it holds in the vast majority of cases. Remarkably, we observe a greater tendency to preserve symmetries as more Rz gates are injected. However, we did not explicitly investigate the relationship between these preserved symmetries and the magic of the initial state. Understanding this connection could provide valuable insights into the role of symmetries and magic in our approach. Future work could explore this relationship in more detail, potentially leading to better ansatz designs and improved performance.

D Pulse optimization of programmable quantum simulators

The Rydberg Hamiltonian can be modeled as

where $\Omega_j(t), \ \phi_j(t), \ \Delta_j(t)$ are Rabi frequency, laser phase, and local detuning of the driving field for each atom, and $|g\rangle$ is the ground state, $|r\rangle$ is the Rydberg state, $\hat{n}_j \, = \, |r\rangle_j \, \langle r|_j$ is the number operator. We use

Figure 8: We analyze optimized Clifford + kRz solution in the eigenbasis for the qubitized number operator $N$ for both LiH (top) and $\text{H}_2\text{O}$ (bottom), demonstrating that the recovered solution approximately conserves of $N$.

Figure 9: We analyze optimized Clifford + kRz solution in the eigenbasis for the qubitized spin operator $S$ for both LiH (top) and $\text{H}_2\text{O}$(bottom), demonstrating that the recovered solution, for most cases, approximately conserves of $S$.

$V_{ij} = C_6/|x_i - x_j|^6$ with $C_6 \; = \; 2\pi \, * \, 8.627 * 10^5$ MHz $\mu m^6$ to describe the Rydberg interactions. In our optimization, we translate the Hamiltonian into Pauli basis, which gives

where $H_{\text{const}} = \sum_{i,j} \frac{V_{ij}}{4} I - \sum_{j} \frac{\Delta_{j}}{2} I$. Then we optimize the pulse in the Pauli basis, such as $\frac{\Omega_j}{2} \cos \phi_j$. We parameterize each pulse as a piece-wise constant function with ten pieces, and initialize them randomly. The atom positions are fixed to be separated with $9.37 \mu m$. Since the interaction strength of the next nearest neighbors are 64 times smaller than the interaction strength of the nearest neighbors, we only consider the leading nearest neighbor interactions for our model.

E Finite temperature calculation

The variational principle extends beyond ground state calculations; indeed, it can be shown that the free energy satisfies a variational principle. That is, defining

where $S(\rho)$ is the von Neumann entropy, it can be shown that for all $\beta$,

for all $\rho$. We take advantage of this variational principle to estimate the free energy at finite temperatures. The procedure is simple: we fix an orthogonal basis of stabilizer states $\{|\psi_k\rangle\mid k=1,\dots,2^n\}$ (typically derived from SBRG), so that the initial state is

where the $p_k \in [0,1]$ form a probability distribution (i.e., $\sum_k p_k = 1$). We construct a Clifford+kRz unitary $U$ from Clifford + kT gates, and evolve each of the under this unitary. The final state is therefore

Defining $E_k \equiv \text{Tr} \big( (H - \mu N) U \, |\psi_k\rangle\!\langle\psi_k| \, U^\dagger \big)$, we observe that

The optimal $p_k$ can be solved for in closed form (i.e., they do not need to be optimized numerically), yielding a free energy estimate

Since calculating $E_k$ as a function of the unitary $U$ is efficient, we can then apply our optimization techniques over the unitary $U$ to minimize $F_{\beta}(\rho)$.