Estimating Molecular Thermal Averages with the Quantum Equation of Motion and Informationally Complete Measurements

Abstract

By leveraging the Variational Quantum Eigensolver (VQE), the “quantum equation of motion” (qEOM) method established itself as a promising tool for quantum chemistry on near-term quantum computers and has been used extensively to estimate molecular excited states. Here, we explore a novel application of this method, employing it to compute thermal averages of quantum systems, specifically molecules like ethylene and butadiene. A drawback of qEOM is that it requires measuring the expectation values of a large number of observables on the ground state of the system, and the number of necessary measurements can become a bottleneck of the method. In this work, we focus on measurements through informationally complete positive operator-valued measures (IC-POVMs) to achieve a reduction in the measurement overheads by estimating different observables of interest through the measurement of a single set of POVMs. We show with numerical simulations that the qEOM combined with IC-POVM measurements ensures satisfactory accuracy in the reconstruction of the thermal state with a reasonable number of shots.

1. Introduction

Computing properties of quantum systems at finite temperatures poses a significant challenge due to its inherent resource-intensive nature. The computational cost of such calculations generally scales exponentially with the system size, unless symmetries or other characteristics of specific systems can be leveraged. The efficient classical computation of thermal averages is, therefore, strongly subjected to the type of system one is interested in. The construction and estimation of properties of thermal states is a task required in different fields, ranging from quantum chemistry [1] and many-body physics [2] to high energy physics [3]. Classical methods typically used to evaluate such quantities include algorithms from the quantum Monte Carlo family [4,5], truncated eigenstate sampling [6], Density Matrix Renormalization group (DMRG) [7], quantum belief propagation [8] and minimally entangled typical thermal states (METTS) [9,10].

In recent years, new types of quantum algorithms including adaptations of classical ones have been developed to leverage the capacity of quantum computers [11], both for the near-term and fault-tolerant eras [12]. Algorithms developed for fault-tolerant quantum computers typically require the use of Quantum Phase Estimation (QPE) [13,14,15] as a sub-routine, while those made for near-term devices are typically based on variational routines such as the Variational Quantum Eigensolver (VQE) [16,17] or Quantum Approximate Optimization Algorithms (QAOA) [18,19].

In this work, we focus on molecular systems such as linear polyenes [20,21,22,23,24,25] as testbeds over which to reconstruct an approximate low-temperature thermal state by means of variational quantum algorithms. Among the methods based on the use of the VQE, the Quantum Equation Of Motion (qEOM) [26,27] is a hybrid method that allows for the evaluation of excited states of a system. More generally, the qEOM belongs to a class of methods referred to as Quantum Subspace Methods (QSM) [28,29,30,31,32], which are becoming increasingly popular for both near-term and fault-tolerant quantum computing for the exploration of the low- and the high-lying excited state, energies, as well as molecular properties [33,34,35] (Figure 1).

The “Equation Of Motion” is a classical chemistry algorithm [36] used to infer a finite number of excited states for a molecular system. To do so, it constructs a subspace of states with a limited number of excitations, which can be conducted by applying excitation operators onto the ground state (GS), from which the eigenstates of the Hamiltonian included in the subspace can be found by minimizing the energy. In this way, a Generalized Eigenvalue Problem (GEP) is formed, which can be solved to recover the eigenstates. In their work [26] Ollitrault et al. introduced the hybrid version of the methods, where the GS is prepared on a quantum computer and the terms of the GEP are measured directly from it. Doing so improves the efficiency in preparing the GS [1,37].

In this work, we demonstrate a potential application for the qEOM method, showing that it can be used to evaluate the thermal state of molecular systems. Compared to other hybrid algorithms capable of performing this task, the qEOM method requires a minimal number of qubits, just the bare minimum necessary to prepare the ground state. In contrast, some methods may require additional resources, such as those that try to emulate open system dynamics [38], aim to reconstruct a purification of the thermal state known as the thermofield double state [18], or simply need additional memory registers like the quantum Metropolis sampling [13]. Additionally, the approach used in the qEOM method to reconstruct the thermal state involves deterministic sampling of the energetically lowest-lying excited states, which is highly efficient at low temperatures [6]. Therefore, our approach is particularly advantageous for the point when near-term quantum computers can efficiently run the VQE routine for ground state preparation, especially if high-temperature thermal state preparation is not the primary focus.

A crucial aspect of the optimal usage of near-term devices is the choice of measurement scheme. Indeed, while in the fault-tolerant era, an efficient measurement scheme based on the QPE is known [39], in the near-term era it is not trivial to find an optimal measurement strategy, as different methods may produce different results when it comes to number of required measurements, circuits depths and accuracy. Promising methods involve many ideas to optimize the measurement process, such as grouping of commuting observables for simultaneous measuring [40,41,42,43,44,45], classical shadows [46,47,48,49], Bayesian inference [50] and machine learning assisted tomography [51]. In combination with the qEOM methods, we will perform measurements using informationally complete positive operator-valued measures (IC-POVMs), which have been recently demonstrated to be a general framework to efficiently estimate many-body operators on near-term devices [52]. The qEOM methods require measuring a number of observables that grow rapidly with the system size. In such cases, IC-POVMs allow for the efficient estimation of different observables using a scalable number of shots. Furthermore, IC-POVMs may be adapted variationally to further lower the variance of the estimation of the observables of interest, i.e., to reduce the number of measurement shots [52]. We point out that while in [52] IC-POVMs are implemented through the use of ancillary qubits, the technique put forward in [53] (“dilation-free IC-POVMs“) does not require any additional qubit, thus preserving one of the advantages of the qEOM method. In this work, we employ the IC-POVM measurement scheme without adaptations, leaving such improvement for future studies. Consequently, the performance of our unadapted method is comparable to that of classical shadow tomography, as both methods allow for the estimation of different observables by measuring a single set of informationally complete POVMs.

This work is organized as follows: in Section 2 we review the (q)EOM methods. In Section 3 we review the IC-POVM measurement scheme. In Section 4 we present our result and in Section 5 we give our outlook and conclusion.

2. Quantum Equation of Motion

In this section, we provide a brief overview of both classical and quantum versions of the Equation of Motion (EOM) method for computing electronic excitation energies [26,27]. The primary goal of the EOM approach is to derive the excited states of a molecule from its GS, represented as $|0\rangle$. This is achieved by constructing and solving an eigenvalue problem within a reduced subspace, limiting the size of the problem to be manageable by conventional classical algorithms. The construction of this reduced subspace involves spanning a set of states derived from applying various excitation operators, including both single and multiple excitations, onto the GS.

In the EOM method, we assume an approximate understanding of the GS of the system, denoted as |0〉. The difference between classical and quantum EOM methods lies in how we determine |0〉. Classically, the GS can be inferred using various methods developed for such a purpose. However, attaining the desired accuracy swiftly becomes computationally prohibitive as the size of the molecule increases. This is primarily due to the limitations in scaling that classical computational methods encounter. To enhance efficiency and achieve superior scaling, it becomes imperative to transcend classical computation and harness quantum effects [37,54,55,56].

A quantum computer may be exploited to obtain a more efficient and faithful GS of the molecule. To achieve this, the fermionic problem (i.e., the second-quantized Hamiltonian of the molecule) must first be transformed into a qubit problem using established mappings [57,58,59,60]. Then, the Hamiltonian can be represented using gates on a quantum computer, and the GS can be prepared, for instance by employing the Variational Quantum Eigensolver (VQE) algorithm [37,61].

Let us introduce the generic excitation operator as ${\widehat{O}}_n^{†}=|n\rangle \langle 0|$, where |n〉 is the nth excited state of the molecule. Analogously, $O_n=|0\rangle \langle n|$. We can readily observe that the expectation value for the excitation energies denoted as $E_{0n}=E_n-E_0$, can be expressed as [26]:

$$\begin{array}{c}\hfill E_{0n}={\displaystyle \frac{\langle 0\left|\left[{\widehat{O}}_n,\left[\widehat{H},{\widehat{O}}_n^{†}\right]\right]\right|0\rangle }{\langle 0\left|\left[{\widehat{O}}_n,{\widehat{O}}_n^{†}\right]\right|0\rangle }}.\end{array}$$

(1)

The key idea of the EOM method lies in the representation of the generic excitation operator $O_n$. We choose to write it as a linear combination of a finite number of single and double excitations denoted by ${\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}$

$$\begin{array}{c}\hfill {\widehat{O}}_n^{†}=\sum _{\alpha }\sum _{{\mu }_{\alpha }}\left[X_{{\mu }_{\alpha }}^{\left(\alpha \right)}\left(n\right){\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}-Y_{{\mu }_{\alpha }}^{\left(\alpha \right)}\left(n\right){\left({\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}\right)}^{†}\right]\end{array}$$

(2)

where $\alpha$ refers to the number of excitations considered ($\alpha =1$ for single excitations and $\alpha =2$ for double excitations), ${\mu }_{\alpha }$ is a collective index for all one-electron orbitals involved in the excitation, and $X_{{\mu }_{\alpha }}^{\left(\alpha \right)}\left(n\right)$, $Y_{{\mu }_{\alpha }}^{\left(\alpha \right)}\left(n\right)$ are the coefficients of the linear combination. For instance, a single excitation is represented by ${\widehat{E}}_{{\mu }_1}^{\left(1\right)}=a_m^{†}{a}_i$ for some orbitals m and i. The span of all the excitations ${\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}$ applied onto the GS forms the subspace over which we solve the EOM problem.

We can finally insert the expansion in Equation (1) to obtain a parametric equation for the coefficients of Equation (2). To find the approximate eigenvalue of the Hamiltonian in the generated subspace, and identify the coefficients of the linear expansion, we need to look for minima of the energy in the coefficients space. We use the variational principle for local minima and evaluate $\delta \left(E_{0n}\right)=0$. By doing so, we obtain a Generalized Eigenvalue Problem (GEP) whose solutions lead to the coefficients we are interested in. The GEP reads:

$$\begin{array}{c}\hfill \left(\begin{array}{cc}M& Q\\ Q\ast & M\ast \end{array}\right)\left(\begin{array}{c}{X}_n\\ Y_n\end{array}\right)=E_{0n}\left(\begin{array}{cc}V& W\\ -W\ast & -V\ast \end{array}\right)\left(\begin{array}{c}{X}_n\\ Y_n\end{array}\right),\end{array}$$

(3)

where the matrix elements can be expressed as expectation values over the GS $|0\rangle$:

$$\begin{array}{cc}\hfill M_{{\mu }_{\alpha }{\nu }_{\beta }}=& \langle 0\left|\left[{\left({\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}\right)}^{†},\widehat{H},{\widehat{E}}_{{\nu }_{\beta }}^{\left(\beta \right)}\right]\right|0\rangle ,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill Q_{{\mu }_{\alpha }{\nu }_{\beta }}=& -\langle 0\left|\left[{\left({\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}\right)}^{†},\widehat{H},{\left({\widehat{E}}_{{\nu }_{\beta }}^{\left(\beta \right)}\right)}^{†}\right]\right|0\rangle ,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill V_{{\mu }_{\alpha }{\nu }_{\beta }}=& \langle 0\left|\left[{\left({\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}\right)}^{†},{\widehat{E}}_{{\nu }_{\beta }}^{\left(\beta \right)}\right]\right|0\rangle ,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill W_{{\mu }_{\alpha }{\nu }_{\beta }}=& -\langle 0\left|\left[{\left({\widehat{E}}_{{\mu }_{\alpha }}^{\left(\alpha \right)}\right)}^{†},{\left({\widehat{E}}_{{\nu }_{\beta }}^{\left(\beta \right)}\right)}^{†}\right]\right|0\rangle ,\hfill \end{array}$$

(7)

with the double commutator defined as $[\widehat{A},\widehat{B},\widehat{C}]=\{[[\widehat{A},\widehat{B}],\widehat{C}]+[\widehat{A},[\widehat{B},\widehat{C}]]\}/2$.

In the classical EOM method, the elements of the GEP are numerically computed from the approximate GS. Similarly, in the qEOM variant, the elements are obtained by performing measurements on the GS prepared on a quantum computer. Note, that the number of observables that need to be estimated grows as $\mathcal{O}\left(N^4\right)$ with the dimension of the Hamiltonian [26]. Then, measuring with IC-POVMs a ground state prepared on a quantum computer is much more efficient than computing classically all the observables of interest, which quickly becomes unfeasible for large molecules. Once the GEP has been reconstructed for both methods, it can be solved through classical algorithms to find a set of excited states of the molecule.

3. Informationally-Complete POVMs

In the qEOM method, once the state of the quantum computer has evolved to the GS of the system of interest, it is necessary to recover the information needed to reconstruct the EOM matrices. Full Quantum State Tomography (QST) is the method used to characterize unknown quantum states, but it is particularly resource-intensive. In most cases, the full knowledge of the state is not necessary, so further reduction in the required resource can be achieved by using techniques that allow for partial state tomography, such as classical shadows [46,62] and, more generally, informationally complete POVMs. In our work, we use local single-qubit IC-POVMs to reconstruct the EOM matrices. In this section, we review how such measurements are performed according to this method, as described in [52,63].

A POVM is said informationally complete if its effects are linearly independent and span the space of linear operators in the Hilbert space of the system. The number of effects of the POVM needs to be at least $d^2$ where d is the dimension of the Hilbert space. For a system of N qubits, this means $4^N$ effects.

An operator $\widehat{O}$ can be decomposed on the basis of the effects of a given POVM $\left\{{\mathsf{\Pi }}_m\right\}$:

$$\begin{array}{c}\hfill \widehat{O}=\sum _m{\omega }_m{\mathsf{\Pi }}_m,\end{array}$$

(8)

which makes its expectation value over the state $\rho$

$$\begin{array}{c}\hfill \langle \widehat{O}\rangle =\left[\rho \widehat{O}\right]=\sum _m{\omega }_m{p}_m,\end{array}$$

(9)

where $p_m=\mathrm{T}\mathrm{r}\left[\rho {\mathsf{\Pi }}_m\right]$ is the probability of the outcome m to be measured.

The first step of the method is to construct a stochastic estimator for the expectation value by performing a Monte Carlo sampling. Specifically, by sequentially measuring the system of a quantum computer S times it is possible to obtain a sequence of outcomes $m_s$ from the distribution $\left\{p_m\right\}$. The weight ${\omega }_m$ can be classically computed in polynomial time. This sequence of samples can then be used to construct the estimator

$$\begin{array}{c}\hfill \overline{O}=\left[\rho \widehat{O}\right]={\displaystyle \frac{1}{S}}\sum _{s=1}^S{\omega }_{m_s}.\end{array}$$

(10)

The statistical error for the procedure is [52]

$$\begin{array}{c}\hfill \sqrt{{\displaystyle \frac{{\langle {\omega }_m^2\rangle }_{p_m}-{\langle {\omega }_m\rangle }_{p_m}^2}{S}}}\end{array}$$

(11)

Ref. [52] proposes to update the choice of the IC-POVM used in the measurement procedure so that the above-defined estimator has a minimal statistical error. In our work, we do not perform this optimization process, while we employ a symmetric IC-POVM (SIC-POVM), which is a reasonable choice in the absence of prior knowledge about the state of the system [52]. Symmetric POVM means that the single qubit effects, when rescaled as ${\tilde{\mathsf{\Pi }}}_i=2{\mathsf{\Pi }}_i$, fulfill the relation $\mathrm{T}\mathrm{r}\left[{\tilde{\mathsf{\Pi }}}_i{\tilde{\mathsf{\Pi }}}_j\right]=(2{\delta }_{ij}+1)/3,\forall i,j$. Specifically, we consider the canonical SIC POVM $\{{\tilde{\mathsf{\Pi }}}_i=|{\tilde{\pi }}_i\rangle {\langle {\tilde{\pi }}_i|\}}_{i=0}^3$ given by the projectors onto $|{\tilde{\pi }}_0\rangle =|0\rangle$, $|{\tilde{\pi }}_k\rangle =|0\rangle +\sqrt{2}{e}^{i2\pi (k-1)/3}|1\rangle$ for $k\in \{1,2,3\}$. We leave the assessment of the performance of adaptive measurements to future work.

A well-known issue of near-term quantum computers is measurement noise, i.e., the effects we are implementing in our measurement apparatus on the quantum device may not be exactly the same as the ideal effects ${\mathsf{\Pi }}_m$, and therefore, the estimation in Equation (10) may be biased. To cope with this issue, one may apply Quantum Detector Tomography (QDT) to estimate the noisy effects on the quantum computer. For instance, we refer the readers to Ref. [64] for the analysis of measurement noise and its correction through QDT for SIC-POVMs.

4. Results

In this section, we present results from simulating the qEOM algorithm on a classical computer. We evaluate the algorithm’s performance on linear polyenes with increasing chain lengths, specifically ethylene and butadiene (${\mathrm{C}}_{2n}{\mathrm{H}}_{2n+2} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}$ n = 1,2). A correlation-consistent (cc-pVDZ) basis set [65] is used for both systems. For each molecule, we obtain molecular orbitals and one- and two-electron integrals using the Hartree–Fock method with the quantum chemistry package PySCF [66,67]. Furthermore, we use the same library for subsequent multiconfigurational self-consistent field (MCSCF) calculations [68]. This allows us to select a subset of active space orbitals to further reduce the size of the systems considered, specifically, we choose CAS(2,2) for ethylene and CAS(4,4) for butadiene, with a complete active space CAS(n_e, n_o) of n_e electrons in n_o orbitals. We then utilize Qiskit [69] to transform the second-quantized fermionic molecular Hamiltonian to a bosonic problem, via the Jordan–Wigner mapping, and perform a VQE routine to find the approximate GS needed for the qEOM algorithm. On top of that, we perform the SIC-POVM as described at the end of Section 3, through the software Aurora v.0.1 [70] to simulate the measurement outcomes of the quantum experiment and to estimate the matrix elements required for the construction of the GEP. Once the eigenstates are obtained, further manipulations are performed with the QuTiP library [71,72].

For comparison purposes, we also diagonalize the Hamiltonian, filtering the eigenstates with a defined particle number [26], and obtain the exact thermal states for the ideal molecular systems, written as:

$$\begin{array}{c}\hfill {\rho }_{th}=\sum _n{e}^{-\beta E_n}|n\rangle \langle n|.\end{array}$$

(12)

For the qEOM algorithm, simulations are performed with both the ideal GS obtained from the diagonalization and an approximate GS obtained through the VQE. The VQE has been run with the UCCSD ansatz [73,74] assuming we are using a noiseless quantum computer. The parameters optimized in the VQE are 3 for the ethylene and 26 for the butadiene. Then, we use IC-POVMs to estimate the observables needed to reconstruct the EOM matrices for both molecules. The number of observables we need to estimate is 100 for the ethylene and 10,816 for the butadiene. Simulations are repeated multiple times to account for statistical errors due to shot noise. Additionally, for the sake of comparison, we also show the unrealistic case of simulations of qEOM with infinite sampling. We demonstrate the agreement between the ideal thermal state and the thermal state reconstructed with the qEOM method using the trace distance between them as a figure of merit [11]:

$$\begin{array}{c}\hfill D(\rho ,\sigma )={\displaystyle \frac{1}{2}}\mathrm{T}\mathrm{r}\left[\sqrt{{(\rho -\sigma )}^{†}(\rho -\sigma )}\right].\end{array}$$

Moreover, we also calculate the absolute value of the energy difference between ideal and reconstructed thermal states, $\Delta E$.

For ethylene, as shown in Figure 2, the method achieves a good approximation regardless of the temperature. The trace distance between the exact thermal state and the thermal state found through the qEOM method, when starting from the exact GS and using infinite sampling, is within numerical error. Noise from limited sampling is more dominant at higher temperatures, where excited states have a heavier weight.

While not reported, statistics comparing trace distance spread for the thermal state in the case where the GS is exact and for the VQE-approximated GS do not display any noticeable difference. Although excited states carry errors propagated from the approximated GS, the scale of these errors is below that of the sampling noise. According to the authors of the original qEOM article [26], for the energy at least, the ratio between the absolute error for the excited state and the GS becomes smaller as the imprecision on the GS increases. In all cases, the excited states are estimated with higher accuracy compared to the GS.

In our case, we find a similar result, as the excited states measured with infinite sampling are much more precise than the VQE-approximated GS. This can be seen from the red curve in Figure 2b, which tends to zero as the temperature increases, where the weight of the approximated GS decreases and so does the error it carries. Conversely, in the other cases displayed in the same plot, the trace distance increases for $\beta$ approaching zero, with the main source of error driving the increment being the finite number of measurement samples.

The different behavior in the cases considered demonstrates that the error due to limited sampling is orders of magnitude higher than the error due to the approximated GS. The error in the excitation energy is always below chemical precision [1], regardless of temperature. While not shown in the plot, our data indicate that the same statement holds true even for a lower number of sampling shots, down to a value of ${10}^4$ shots.

For butadiene, similar plots are reported in Figure 3. For the trace distance in the case of an exact GS, see Figure 3a, the curve with infinite sampling has a value of zero within numerical accuracy for values of $\beta$ over approximately 40, under which the limited number of excited states limits the accuracy of the approximation. The other curves with limited sampling follow the trend, with evident lower accuracy. In the case of the VQE-approximated GS, the curve with infinite sampling is lower bounded by the accuracy of the GS at higher values of $\beta$.

We report that the instability of the qEOM method, subjected to the type of system it is being applied to, is also influenced by the sampling technique we included in the method. According to the original paper, when the lowest energy solution is close to the GS, the conditioning of the EOM matrices deteriorates, worsening the approximation of the energies. Our results indicate that this effect is also influenced by other sources of error included in the method. In our case, the estimation of the EOM matrices through measurement is also affected by the sampling technique.

In our study, we report that simulations with a low number of sampling shots, up to ${10}^4$ shots, were potentially unstable due to the low accuracy of the estimated EOM matrices. As for the energy, we included a line for the chemical precision at $1.6\times {10}^{-3}$ Ha [1] on the plot for comparison purposes. Our results seem to indicate that the point where the line for the chemical precision intersects the estimated energy error does not depend significantly on the number of shots or on the choice of the GS.

These variables may influence the accuracy of the estimated energy, especially in the region below chemical precision. In the low-temperature region, we see some interplay between shot noise and accuracy, but even with limited sampling the estimated energy always falls within chemical precision. At higher values of $\beta$, eigenvalues found using single and double excitation are not enough to accurately estimate the energy, and the approximation falls off.

For this reason, in the high-temperature region, the errors accumulated from the limited sampling and the approximated GS becomes small enough that all curves displayed in Figure 3c,d overlap. While for the energy this may indicate that a higher number of shots may not improve the overall accuracy below a certain temperature, this statement may not hold for all possible quantities of interest, as unlike the energy, in the trace distance differences in accuracy can be seen at higher values of $\beta$ (in the range of values between $\beta =5$ and $\beta =10$).

5. Conclusions

In this work, we investigated the effectiveness of the qEOM algorithm for computing molecular excited states in conjunction with the IC-POVM sampling technique. We applied this method to linear polyenes such as ethylene and butadiene. Measuring through IC-POVMs is particularly convenient because it allows for the estimation of several observables at the same time, using the same dataset of raw measurement data. This procedure is much more efficient than the independent estimation of each observable necessary for building the EOM problem, as their number grows rapidly with the system size. It is also clearly more efficient than the computation of the mean values of these observables through the classical method, which becomes unfeasible for large molecules.

For small molecules like ethylene, we observe robust numerical agreement between the ideal GS and the approximated GS found through the qEOM at any temperature, and for all the number of shots we have considered, which are meaningful for current near-term computers. The estimated average energy of the thermal state of the molecule remains below chemical precision across all temperatures.

For butadiene, the accuracy of our method depends on the number of measurement shots. For a low number of shots (${10}^4$ or lower), we have found numerical instabilities in the solution of the qEOM problem for this molecule. Using ${10}^5$ shots or more, which is doable on near-term quantum computers, the method becomes stable. In this case, we observe a delicate interplay between the number of shots and the accuracy of estimates, particularly at intermediate values of $\beta$. Notably, in this situation, the approximation consistently delivers estimates of energy with precision surpassing chemical precision.

In conclusion, despite the challenges posed by larger molecular systems, our findings underscore the utility of the qEOM algorithm in accurately computing molecular thermal states. As advancements in quantum computing continue, methods that harness the potential of such devices promise to outperform classical computation in a wide class of tasks, including molecular simulations at finite temperatures. In this work we showed that the qEOM method can be run more efficiently through the IC-POVM sampling technique, paving the way for further test and research on the matter. Further work may include using the adaptive IC-POVM scheme as described in [52] and testing the effect of noise on the preparation of the GS. Additionally, different measurement schemes may provide advantages and disadvantages which may be worth investigating. Furthermore, simulations for bigger molecules may be run on quantum hardware.