Critical behavior of Ising model by preparing thermal state on quantum computer

TL;DR

This paper demonstrates simulating the critical behavior of the Ising model using a quantum computer to prepare thermal states via variational quantum imaginary time evolution, showing results consistent with exact methods.

Contribution

It introduces a systematic method for designing variational ansatz in QITE and applies it to study critical phenomena in the Ising model on a quantum computer.

Findings

Successful simulation of Ising criticality on small lattices

Quantum results agree with exact diagonalization near critical points

Method demonstrates potential for quantum simulation of phase transitions

Abstract

We simulate the critical behavior of the Ising model utilizing a thermal state prepared using quantum computing techniques. The preparation of the thermal state is based on the variational quantum imaginary time evolution (QITE) algorithm. The initial state of QITE is prepared as a classical product state, and we propose a systematic method to design the variational ansatz for QITE. We calculate the specific heat and susceptibility of the long-range interacting Ising model and observe indications of the Ising criticality on a small lattice size. We find the results derived by the quantum algorithm are well consistent with the ones from exact diagonalization, both in the neighbourhood of the critical temperature and the low-temperature region.

Equations

$$H=-\sum _{i>j\in \Lambda }\frac{J}{r_{ij}^{\alpha }}{Z}_i{Z}_j-h\sum _i{Z}_i,$$

$$r_{ij}=\sum _{d=1}^D\min (|r_d^i-r_d^j|,N_d-|r_d^i-r_d^j|).$$

$$\rho =\frac{1}{Z_{\beta }}{e}^{-\beta H},Z_{\beta }\equiv \text{tr}(e^{-\beta H}).$$

$$⟨O⟩\equiv \text{tr}(\rho O).$$

$$C_v\equiv \frac{1}{|\Lambda |}\frac{\partial ⟨H⟩}{\partial T}=\frac{1}{|\Lambda |T^2}\left[⟨H^2⟩-⟨H{⟩}^2\right],$$

$$⟨M⟩\equiv ⟨Z_{tot}⟩,$$

$$\displaystyle\chi\equiv\frac{1}{|\Lambda|}\frac{\partial\langle M\rangle}{\partial h}\Big{|}_{h=0}=\frac{1}{|\Lambda|T}\left[\langle Z_{tot}^{2}\rangle-\langle Z_{tot}\rangle^{2}\right].$$

$$H^2,H,Z_{tot}^2,Z_{tot}$$

$$|\psi (\tau )⟩=\frac{e^{-\tau H}|\psi (0)⟩}{\sqrt{⟨\psi (0)|e^{-2\tau H}|\psi (0)⟩}},$$

$$|\psi (\tau )⟩⟨\psi (\tau )|=\frac{1}{Z_{2\tau }}{e}^{-2\tau H},Z_{2\tau }\equiv \text{tr}(e^{-2\tau H}).$$

$$\rho (\tau )=\frac{e^{-\tau H}|\psi (0)⟩⟨\psi (0)|e^{-\tau H}}{⟨\psi (0)|e^{-2\tau H}|\psi (0)⟩}.$$

$$\text{tr}(\rho (\tau )O)=⟨\psi (\tau )|O|\psi (\tau )⟩.$$

$$\frac{\mathrm{d}\rho (\tau )}{\mathrm{d}\tau }=\mathcal{L}[\rho (\tau )],$$

$$L^2\equiv {\left|\left|\frac{\mathrm{d}\hat{\rho }(\tau )}{\mathrm{d}\tau }-\mathcal{L}[\hat{\rho }(\tau )]\right|\right|}^2,$$

$$L^2={\left|\left|\sum _{\mu }\frac{\partial \hat{\rho }(\theta )}{\partial {\theta }_{\mu }}{\dot{\theta }}_{\mu }-\mathcal{L}(\hat{\rho })\right|\right|}^2.$$

$$\delta L^2=0\Rightarrow \frac{\partial L^2}{\partial {\dot{\theta }}_{\mu }}=\sum _{\nu }{M}_{\mu \nu }{\dot{\theta }}_{\nu }-V_{\mu }=0,$$

$$M_{\mu \nu }$$

$$V_{\mu }$$

$$\sum _{\nu }{M}_{\mu \nu }{\dot{\theta }}_{\nu }=V_{\mu }.$$

$$\vec{\theta }({\tau }_0+\delta \tau )\simeq \vec{\theta }({\tau }_0)+\dot{\vec{\theta }}\delta \tau ,$$

$$\{|+⟩,|-⟩{\}}^{\otimes 2}=\{|++⟩,|+-⟩,|-+⟩,|--⟩\}.$$

$$X|+⟩=|+⟩,X|-⟩=-|-⟩.$$

$$⟨i(\tau )|O|i(\tau )⟩=\frac{⟨i|e^{-\tau H}O{e}^{-\tau H}|i⟩}{⟨i|e^{-2\tau H}|i⟩}.$$

$$⟨O⟩=\sum _i{p}_i⟨i(\tau )|O|i(\tau )⟩,$$

$$p_i\equiv \frac{⟨i|e^{-2\tau H}|i⟩}{Z_{2\tau }}.$$

$$O=\sum _m{h}_m{\tilde{Z}}_m.$$

$$⟨O⟩=\sum _m{h}_m⟨\widetilde{+}(\tau )|{\tilde{Z}}_m|\widetilde{+}(\tau )⟩,$$

$$|\psi ({\tau }_0+\Delta \tau ))⟩=\frac{e^{-\Delta \tau H}|\psi ({\tau }_0)⟩}{\sqrt{⟨\psi ({\tau }_0)|e^{-2\Delta \tau H}|\psi ({\tau }_0)⟩}}.$$

$$|\psi ({\tau }_0+\Delta \tau )⟩=e^{-i\Delta \tau \hat{A}({\tau }_0)}|\psi ({\tau }_0)⟩,$$

$$\hat{A}({\tau }_0)=\sum _{i_1\dots i_{N_q}}{a}_{i_1\dots i_{N_q}}^{({\tau }_0)}{\sigma }_{i_1}\dots {\sigma }_{i_{N_q}}\equiv \sum _I{a}_I^{({\tau }_0)}{\tilde{\sigma }}_I,$$

$$e^{-i\Delta \tau {\sum }_I{a}_I^{({\tau }_0)}{\tilde{\sigma }}_I}=\prod _I{e}^{-i\Delta \tau a_I^{({\tau }_0)}{\tilde{\sigma }}_I}+\mathcal{O}(\Delta {\tau }^2),$$

Full text

I Introduction

With the development of quantum devices and quantum algorithms, it is possible to solve problems on quantum computers that are hard for classical ones. Quantum computers have already been successfully implemented in many fields, including quantum chemistry, condensed matter physics and lattice field theory, see references McArdle_20 ; Klco_2018 ; Klco_2020 ; Ciavarella_2021 ; Funcke_2022 ; Clemente_strategies ; Banuls2020 as some examples. With the growing number of qubits and improved fidelities of quantum devices, more realistic physical models can be tackled, and the potential of quantum computers can be explored. As an example of application, in this article, we prepare the thermal state of the Ising model with a quantum algorithm at various temperatures, including points close to the critical temperature and the low-temperature region. To demonstrate the feasibility of our approach, we compare the quantum simulation results of the chosen physical quantities with the results from classical simulations.

Numerous algorithms have been proposed to enable a quantum computer to prepare a thermal state. These include the quantum thermal dynamic method, where the target system is coupled with a bath at equilibrium Terhal_2000 , variational quantum algorithm based on the thermofield double state Wu_2019 ; 21npj_Sagastizabal , as well as many quantum imaginary time evolution(QITE) algorithms such as the one utilizing Hubbard-Stratonovich transformation Chowdhury_16 , QITE based on variational ansatz (QITE-ansatz) McArdle_19 , QITE based on measurement (QITE-measure) Motta_20 and QITE by performing coordinatewise optimization Benedetti_21 . The scope of our research is to focus on the usage of noisy intermediate-scale quantum (NISQ) devices Preskill_2018 ; Bharti22 . Given the presence of quantum noise, it is necessary to minimize the depth of the quantum circuits. We utilize the QITE-ansatz algorithm to generate thermal states in our research, as it has a relatively shallower circuit depth in comparison to other algorithms mentioned previously. In QITE-ansatz algorithm, the imaginary time evolution is carried out on a prior parameterized quantum circuit, and the parameters are evolved variationally. Thus, the parameterized quantum circuit is usually called variational ansatz. The variational ansatz is designed for ground state preparation in most references utilizing QITE-ansatz, such as McArdle_19 ; Liu2021 ; Zoufal_2021 . Here, for thermal state preparation, we propose to construct a variational ansatz converted from quantum circuits utilized in QITE-measure Motta_20 . The circuit in QITE-measure can also carry out imaginary time evolution, but the circuit depth is quite large. The circuit depth can be much reduced by converting the circuit into a variational ansatz. For example, when simulating the Ising model, the quantum circuits in QITE-measure have $\sim 100$ layers, while the variational ansatz circuits used in this work have less than 10 layers.

In this article, we study the long-range interacting Ising model. Long-range interaction between spins is introduced naturally in trapped-ion spin systems Islam_2013 , and its dynamics can be simulated utilizing quantum simulation algorithms. The long-range interaction also leads to interesting physics such as confinement Liu_2019 and meson scattering vovrosh2022dynamical . Meanwhile, the long-range interaction leads to effective dimensions that impact the system’s critical behavior. Here, we calculate the specific heat of the long-range interacting Ising model near the critical point and in the low-temperature region.

This article is organized as follows. In section II, we introduce the long-range interacting Ising model and the measurement method of relevant physical quantities on a quantum computer. In section III, we discuss the process of thermal state preparation using QITE-ansatz algorithm in detail, especially the method of variational ansatz design. In section IV, we present the numerical results and discuss the observed indications of the criticality. Finally, in section V, we summarize the techniques used in this article and discuss the possible extension for further works.

II Long-range interacting Ising model

We consider the $D=2$ dimensional Ising model on a square lattice $\Lambda$ with long-range interactions. The Hamiltonian reads

[TABLE]

where $Z_{i}$ is the Pauli-$Z$ operator on the $i$th spin. $J$ is the bare coupling strength, and $\alpha$ denotes the range of the interaction. $h$ denotes the strength of the longitudinal external field. The distance $r_{ij}$ is defined by the Manhattan distance under periodic boundary condition(PBC): Assuming the position of spin $i$ on the square lattice is represented by integer vector $\vec{r}^{i}=(r^{i}_{1},\ldots,r^{i}_{D})$ and the volume of the lattice is $|\Lambda|=N_{1}\times\ldots\times N_{D}$, then

[TABLE]

This Hamiltonian is a generalization of the interaction part of the Hamiltonian introduced in reference Liu_2019 . It reduces to the original nearest-neighbor Ising model (NNIM) in the limit $\alpha\rightarrow\infty$. Because in the limit $\alpha\rightarrow\infty$, all the long-range couplings $J/r_{ij}^{\alpha}$ with $r_{ij}>1$ vanish, except for nearest-neighbor ones with $r_{ij}=1$.

The state of the Ising system at a finite temperature is described by the density operator. Its equilibrium state is the Gibbs state of which the density operator reads

[TABLE]

Here $\beta$ is the inverse temperature $\beta\equiv 1/(k_{B}T)$ and we define $K\equiv J\beta$ for later convenience. For an arbitrary observable $O$, its expectation value of the thermal state is given by

[TABLE]

This article targets the case where the expectation values are evaluated for different $K$ and a zero external field $h=0$.

Now we exhibit observables to compute the Ising model’s specific heat and susceptibility. Analyzing these measures allows us to examine the critical behavior of the Ising model. The specific heat is defined by the changing rate of the internal energy in a unit volume when varying the temperature $T$. It can be evaluated by the energy-fluctuation relation:

[TABLE]

where the last expression can be derived by taking the Gibbs state Eq. (3) to evaluate the expectation values.

Similarly, the susceptibility is defined by the changing rate of the magnetization in a unit volume with respect to the external field strength $h$ (evaluated at $h=0$). The total magnetization is given by

[TABLE]

where $Z_{tot}\equiv\sum_{i}Z_{i}$, i.e., the sum of all the spins in the lattice. Then the susceptibility can be evaluated according to the susceptibility-fluctuation relation

[TABLE]

In summary, evaluating the specific heat and susceptibility is equivalent to calculating the expectation values of the corresponding operators. The operators to be measured include

[TABLE]

which can all be reduced to linear combinations of Pauli operators. To evaluate the expectation values of the above operators on quantum computers, we can generate the thermal state utilizing a quantum algorithm and then evaluate the expectation values of the Pauli operators. Notice that for the above operators, the elementary Pauli operators can be written as products of Pauli-$Z$ operators, so they commute and can be measured simultaneously on the quantum computer. Combined with the fact that the Hamiltonian in Eq. (1) consists of only Pauli-$Z$ operators, we can simplify the initial state to be evolved on quantum computer. It enables us to simulate the system on a larger lattice. For general models, such as the Ising model with a transversal field, the simplification does not hold. More details can be found in section (III.1).

III Thermal state preparation with quantum imaginary time evolution

One can use the quantum imaginary time evolution(QITE) algorithm to prepare a thermal state 111We have also tried non-QITE-based variational quantum algorithm proposed in reference Wu_2019 ; 21npj_Sagastizabal . We find its accuracy is not as good as the QITE-ansatz algorithm when the simulation is performed on a noiseless statevector quantum simulator, as demonstrated in previous studies McArdle_19 ; Motta_20 . This section provides an explanation of the QITE-ansatz algorithm. QITE-ansatz algorithm is designed to evolve an $N_{q}$-qubit quantum state $\ket{\psi(0)}$ to

[TABLE]

where $\tau$ is a real number denoting imaginary time. The denominator is a normalization factor to guarantee the evolution’s unitarity. Assuming we have the quantum circuit to carry out the unitary evolution, then by choosing the initial state to be the maximally mixed state (defined as the density operator) $\ket{\psi(0)}\bra{\psi(0)}=\mathbf{I}/\mathbf{d}$ 222Here we abuse bra-ket notation for the mixed state. It will be explained in detail in the subsection (III.1) ($\mathbf{I}$ is the identity operator of the $\mathbf{d}\equiv 2^{N_{q}}$ dimensional Hilbert space), one finds the final state is the thermal state with inverse temperature $\beta=2\tau$

[TABLE]

The QITE-ansatz algorithm was proposed in references McArdle_19 ; Yuan_2019 . This technique is originally used to project out the ground state of the Hamiltonian according to Eq. (9). It has been successfully implemented in the field of quantum chemistry, quantum field theory and machine learning, see e.g. McArdle_20 ; Liu2021 ; Zoufal_2021 .

Following Yuan_2019 , we first review the QITE-ansatz algorithm within the density operator formalism. The density operator of Eq. (9) reads

[TABLE]

The mathematical description of a quantum state with the density operator is equivalent to that with the pure state. In particular, the expectation values of any observable $O$ coincide

[TABLE]

The imaginary time evolution of the density operator follows the von-Neumann equation Yuan_2019

[TABLE]

where $\mathcal{L}$ is the Liouville operator defined by $\mathcal{L}(\rho)=-\{H,\rho\}+2\tr(\rho H)\rho$ with anti-commutator $\{H,\rho\}=H\rho+\rho H$. As the Hilbert space of the whole $N_{q}$ qubits is hard to be explored by a quantum circuit, we utilize a density operator $\hat{\rho}(\tau)=\ket{\phi(\tau)}\bra{\phi(\tau)}$ to approximate the target density $\rho(\tau)$. The approximation $\hat{\rho}(\tau)$ satisfies the following requirements: (1) It has the same initial state $\hat{\rho}(0)=\rho(0)=\ket{\psi(0)}\bra{\psi(0)}$. (2) The evolution of $\hat{\rho}(\tau)$ approximately satisfies the von-Neumann equation $\mathrm{d}\hat{\rho}(\tau)/\mathrm{d}\tau-\mathcal{L}[\hat{\rho}(\tau)]=0$.

The approximation $\hat{\rho}(\tau)$ is generated with a variational ansatz $\ket{\phi(\vec{\theta}(\tau))}=U(\vec{\theta}(\tau))\ket{\psi(0)}$, where $\vec{\theta}$ is a real variational parameter vector with $N$ components. $U(\vec{\theta})=U_{N}(\theta_{N})\ldots U_{1}(\theta_{1})$ is a series of parameterised unitary quantum gates. According to the first requirement mentioned above, $U(\vec{\theta}(0))$ should be the identity operator $\mathbf{I}$. With the variational ansatz, the evolution of the quantum state is converted to the evolution of the variational parameters $\vec{\theta}$. However, as the variational ansatz cannot explore the whole Hilbert space, $\ket{\phi(\vec{\theta}(\tau))}$ can not fulfill the von-Neumann equation exactly. Instead, we demand that the von-Neumann equation is fulfilled sufficiently well according to the second requirement. The violation of the von-Neumann equation is measured by the McLachlan distance $L^{2}$, which is defined by

[TABLE]

where $||A||^{2}=\tr(A^{\dagger}A)$ represents Frobenius norm. According to the differential chain rule, we have

[TABLE]

So that the McLachlan distance is a quadratic function of the time derivatives of the variational parameters $\dot{\theta}_{\mu}\equiv\partial\theta_{\mu}/\partial\tau$. $L^{2}$ can be minimized with the variational principle, which leads to

[TABLE]

where

[TABLE]

Here $M$ is a $N\times N$ matrix while $V$ is a $N$ dimensional vector. Following McArdle_19 ; Funcke2021dimensional , one can construct some specific quantum circuits to measure $M$ and $V$, which cost $\mathcal{O}(N^{2})$ quantum device calls and one additional ancilla qubit.

After deriving $M$ and $V$, we can construct the following linear equations

[TABLE]

Then one can solve for the time derivative of the variational parameters $\dot{\theta}_{\nu}|_{\tau=\tau_{0}}$ at a given imaginary time $\tau_{0}$, utilizing methods such as pseudo-inverse McArdle_19 . The variational parameters at the next time slice $\tau_{0}+\delta\tau$ are given according to the Euler method

[TABLE]

where $\dot{\theta}_{\nu}=\sum_{\mu}M^{-1}_{\nu\mu}V_{\mu}$.

The computational complexity of the QITE-ansatz grows polynomially with the number of variational parameters $N$. In each time slice, the time complexity of solving linear equations grows polynomially with $N$, while the matrix $M$ and vector $V$ can also be evaluated using quantum computers within polynomial time. Thus as long as $N$ grows polynomially with the system size $N_{q}$, the time complexity of the QITE-ansatz grows polynomially with $N_{q}$ and can be extended to large-scale quantum systems. The following subsections will introduce how to prepare the maximally mixed state and choose an appropriate variational ansatz.

III.1 Initial state preparation

Here we introduce how to prepare the initial state as the maximally mixed state $\mathbf{I}/\mathbf{d}$. Quantum circuits are suitable for generating pure states. We need some strategies to generate mixed states utilizing pure states. As discussed in White2009 , there are two strategies: ancilla pair state (APS) and classical product state (CPS). Both strategies can be used to prepare maximally mixed state $\mathbf{I}/\mathbf{d}$. However, preparing $\mathbf{I}/\mathbf{d}$ with APS doubles the number of qubits to $2N_{q}$ Zoufal_2021 . It also introduces some complexities in variational ansatz design to evolve the pair state.

Instead, we can prepare the maximally mixed state via CPS, which reduces the required qubits to $N_{q}$. The maximally mixed state $\mathbf{I}/\mathbf{d}$ describes that the probabilities of sampling every basis vector from a given orthogonal basis are the same, where each basis vector is a pure state. As the maximally mixed state is unitarily invariant $U(\mathbf{I}/\mathbf{d})U^{-1}=\mathbf{I}/\mathbf{d}$, the orthogonal basis can be chosen arbitrarily. To generate the thermal state, it is recommended in White2009 to use a basis formed by classical product states, such as $\{\ket{+},\ket{-}\}^{\otimes N_{q}}$, where $\{\cdot\}^{\otimes N_{q}}$ represents a set generated by the $N_{q}$ times tensor product of each element in $\{\cdot\}$. For example,

[TABLE]

Here $\ket{+}$,$\ket{-}$ represent the eigenvectors of the Pauli-$X$ operator

[TABLE]

If we use the classical product state as the initial state, the thermal expectation value $\langle O\rangle$ can not be measured straightforwardly due to the normalization factor in Eq. (9). Assume that we take the orthogonal basis as $\{\ket{i}\}$. Evolving all basis vectors $\ket{i}$ for imaginary time $\tau$, one gets the expectation values of an observable $O$, which read

[TABLE]

Usually, the denominators would be different for different basis vectors $\ket{i}$. To derive the thermal expectation value $\langle O\rangle$ in Eq. (4), we should multiply the above expectation values with coefficients $\{p_{i}\}$

[TABLE]

where $p_{i}$ is defined by

[TABLE]

Here $\{p_{i}\}$ can be treated as a probability distribution, as they are all positive and satisfy the normalization condition $\sum_{i}p_{i}=1$. To evaluate the thermal expectation value of the operator $O$, as mentioned in Motta_20 , we do not need to calculate all the $\{p_{i}\}$(which would be impossible to calculate, as the number of $p_{i}$ grows exponentially with the number of qubits). With the minimally entangled typical thermal state(METTS) algorithm proposed by Stoudenmire and White Stoudenmire_2010 , one can sample $\{\ket{i}\}$ according to the distribution $\{p_{i}\}$. The thermal expectation value $\langle O\rangle$ is the average of the expectation of $O$ with the time-evolved sampled vectors. In conclusion, though imaginary time evolution with CPS as initial states requires the number of qubits equal to the system size, one has to evolve different initial states $\ket{i}$ to acquire statistics. On the other hand, imaginary time evolution with APS as an initial state doubles the number of qubits while evolving only one initial state.

However, the situation gets simplified when we consider the classical Ising model and the observables in Eq. (8), which consist of Pauli-$Z$ operators. The observables can be generally expressed as

[TABLE]

Here $\tilde{Z}_{m}$ represents the tensor product of $Z$ operators at some sites and identity operators at the others, such as $\tilde{Z}_{m}=Z_{N_{q}-1}\ldots I_{1}Z_{0}$. In Appendix A, we prove that the thermal expectation value of $O$ can be calculated according to

[TABLE]

where $\ket{\tilde{+}(\tau)}$ is imaginary time evolved state according to Eq. (9). The state is initialized as $\ket{\tilde{+}(0)}=\ket{\tilde{+}}$, where $\ket{\tilde{+}}\equiv\ket{+}^{\otimes N_{q}}$ is the $N_{q}$-fold tensor product of $\ket{+}$ in Eq. (21). Thus for the Ising model, we only need to calculate the imaginary time evolution with the initial state $\ket{\tilde{+}}$.

In this work, we use $\ket{\tilde{+}}$ as the initial state to present our results. For general models, such as the Ising model with a transversal field, the above simplification does not hold. We need to sample the classical product states using the METTS algorithm or utilize the ancilla pair state.

III.2 Variational ansatz design

Choosing a proper variational ansatz is a cornerstone for the success of the QITE-ansatz algorithm Bharti22 . In most literature on QITE-ansatz, the variational ansatz is designed to prepare the ground state of a Hamiltonian, and it is suitable to evolve some specific initial states, such as the unitary coupled cluster ansatz evolving Hartree-Fock states McArdle_20 . Focusing on thermal state preparation and the initial state introduced in the previous section, we propose to construct a variational ansatz converted from quantum circuits utilized in the QITE-measure algorithm proposed by Motta et. al. Motta_20 .

We briefly introduce how to construct the quantum circuits used in the QITE-measure algorithm. The goal of QITE-measure is also evolving an initial state $\ket{\psi(0)}$ according to Eq. (9). Consider evolving the state $\ket{\psi(\tau_{0})}$ for a small time slice $\Delta\tau$

[TABLE]

As this transformation is unitary, we can always find a Hermitian operator $\hat{A}(\tau_{0})$ such that

[TABLE]

and $\hat{A}(\tau_{0})$ can be expanded in a complete Pauli basis

[TABLE]

where the expansion coefficients $a_{i_{1}\ldots i_{N_{q}}}^{(\tau_{0})}$ are real due to the Hermicity of $\hat{A}(\tau_{0})$, and $\sigma_{i_{j}}=I,X,Y,Z$ corresponding to $i_{j}=0,1,2,3$ is the single-qubit Pauli operator on the site $j$, and we call the tensor product of the single-qubit Pauli operator, $\tilde{\sigma}_{I}$ as Pauli string. For this reason, the single-qubit Pauli operator is sometimes called Pauli letter Oliver_22 . For each imaginary time $\tau_{0}$, one can calculate all the expansion coefficients $a_{I}^{(\tau_{0})}$ by evaluating the expectation values of some observables with respect to the quantum state $\ket{\psi(\tau_{0})}$. The observables are the composition of Pauli strings and the Hamiltonian (See more details in Motta_20 ). Notice that the transformation in Eq. (28) can be approximated by

[TABLE]

where the product consists of several Pauli exponentials which have the form $e^{-i\theta\tilde{\sigma}_{I}}$, and the Pauli exponential can be realized with quantum gates in a standard way Nielsen2000 . Thus, the whole quantum circuit used in the QITE-measure can be constructed using several Pauli exponentials for each time slice. In the last time slice, the circuit depth is proportional to the final imaginary time $\tau$.

Notice that if a system has $N_{q}$ qubits, the total number of Pauli strings on these qubits is $4^{N_{q}}$. Thus the number of Pauli exponentials required for evolving each time slice seems exponential as a function of system size according to Eq. (29). However, the situation gets simplified when the Hamiltonian $H$ consists of some local interaction terms

[TABLE]

where each $H_{m}$ acts on a local set of qubits, and the number of $H_{m}$ is polynomial as a function of system size. For example, $H_{m}\propto Z_{i}Z_{j}$ and the number of $H_{m}$ is $\mathcal{O}(N_{q}^{2})$ in case of long-range interacting Ising model. Though the local terms $H_{m}$ may not commute, the imaginary time evolution $e^{-\Delta\tau H}$ can be decomposed by

[TABLE]

Then the previous steps in QITE-measure can be implemented for each $e^{-\Delta\tau H_{m}}$. As shown in Motta_20 , when the Hamiltonian consists of local terms and the correlation length of the system is finite, the expansion in Eq. (29) for each $H_{m}$ can be implemented with Pauli strings on a support constantly larger than the support of $H_{m}$ (Support of a Pauli string is defined by the set of qubits on which the Pauli letters are not identity). The correlation length of a system is finite when its Hamiltonian is outside the critical region. Thus the support of the Pauli strings has no dependence on the system size, and the total number of Pauli exponentials $e^{-i\theta\tilde{\sigma}_{I}}$ is a polynomial function of the system size at least when the Hamiltonian is sufficiently far away from the critical point.

Compared with the QITE-ansatz, the precision of the QITE-measure is not limited by the variational ansatz. However, the circuit depth grows linearly with the evolution time $\tau$. Thus this algorithm would be very sensitive to coherent or incoherent noise in real quantum devices and can only be applied to small spin systems Aydeniz_2020 .

Quantum circuits constructed in QITE-measure can be naturally converted into a variational ansatz with the following steps: (1) using all the necessary Pauli exponentials at one time slice as one layer of the variational ansatz; (2) sequentially repeating the layer several times in the quantum circuit; (3) converting all the expansion coefficients $a_{I}^{(\tau_{0})}$ into undetermined parameters, which are initially zero and to be evolved according to the QITE-ansatz algorithm. Times of repetition for one layer is called the depth of the variational ansatz, also called the number of layers.

The behavior of this variational ansatz can be analyzed with the help of QITE-measure. Assuming we have the same quantum circuit layers for the variational ansatz in QITE-ansatz and the quantum circuits in QITE-measure. Because the states prepared in QITE-measure can all be explored by the variational ansatz, one can expect QITE-ansatz using this circuit to behave at least better than QITE-measure. The systematic error of the QITE-measure circuit is of the first-order Trotter type, i.e., $\mathrm{error}\sim\mathcal{O}(\Delta\tau)$ Motta_20 . By equalizing the longest circuit depth used in QITE-measure and the depth in variational ansatz, it can be deduced that in the worst case, the variational ansatz leads to an error of $\mathcal{O}(1/\mathrm{L})$, where $L$ is the number of layers.

In the numerical simulations, we find that the circuit depth required in QITE-ansatz is much smaller than that required in the QITE-measure. For example, in our numerical simulation of the Ising model, if the imaginary time of the final state is $\tau=0.5$, with step size $\Delta\tau=0.002$, QITE-measure requires the number of layers $\tau/\Delta\tau=250$. In contrast, to reach a sufficiently good precision using the variational ansatz, we find the number of layers required is at most $L=N_{d}$ for the 2-D nearest neighbor Ising model where $N_{d}$ is the side length of the Ising lattice system. More details on the number of circuit layers required are shown in Appendix B.

The variational ansatz can be simplified due to some special structures of the Hamiltonian and the initial state. In the numerical simulations, we notice that some of the variational parameters are always zero during the whole evolution, which corresponds to the same set of Pauli strings over all the layers. We call the Pauli string in this set irrelevant, and the other Pauli strings corresponding to non-zero variational parameters are relevant. As the irrelevant Pauli exponentials are identity, they can be removed a priori when constructing the variational ansatz. These irrelevant Pauli strings can be identified according to the symmetry and some special structures of the Hamiltonian and the initial state. For example, if all the entries in the Hamiltonian and the initial state are real, then the corresponding unitary operator $e^{-i\Delta\tau\hat{A}}$ should also be real. Thus all Pauli strings with an even number of Pauli-$Y$ letters are irrelevant.

We demonstrate the above construction of variational ansatz using an example of a two-qubit ($N_{q}=2$) Ising system. There are $4^{2}=16$ Pauli strings on the two-qubit system. Assume we have an initial state $\ket{++}$ and the system Hamiltonian $H=-Z_{1}Z_{0}$. Because all the entries in the Hamiltonian and the initial state are real, eliminating Pauli strings with an even number of Pauli-$Y$ letters leaves 6 Pauli strings: $I_{1}Y_{0},X_{1}Y_{0},Y_{1}I_{0},Y_{1}X_{0},Z_{1}Y_{0},Y_{1}Z_{0}$. Evolving one layer with these 6 Pauli strings using QITE-ansatz, we further find 4 Pauli strings are irrelevant. It leaves only two relevant Pauli strings for the imaginary time evolution

[TABLE]

One can verify that

[TABLE]

with expansion coefficients

[TABLE]

In the QITE-measure algorithm, to evolve the initial state to an arbitrary time $\tau$, the quantum circuit is shown in figure 2a. It has $\tau/\Delta\tau$ layers. The variational ansatz with $L$ layers for the two-qubit Ising system is constructed as shown in figure 2b. In this circuit, $\{\theta_{1},\theta_{1}^{\prime}\ldots\theta_{L},\theta_{L}^{\prime}\}$ are all variational parameters, taking zero as initial values, and to be evolved according to the QITE-ansatz algorithm.

IV Numerical results

In this section, we apply the previous variational ansatz design procedure to the long-range interacting Ising model, where we prepare CPS $\ket{\tilde{+}}$ as the initial state. Equipped with the thermal state, we can calculate the specific heat $C_{v}$ and susceptibility $\chi$ as a function of $K\equiv J\beta$. The measurement process for $C_{v}$ and $\chi$ is shown schematically in Figure 1. The QITE-ansatz algorithm for Ising model is compactly shown in Algorithm 1. Our numerical simulations are carried out on the Qiskit noiseless statevector quantum simulator Qiskit .

The initial state and variational ansatz are chosen as described in section III. To calculate the thermal expectation values of the Ising model, we only need to calculate the imaginary time evolution of the product state $\ket{\tilde{+}}$. With the initial state, and for every local interaction term in Ising model $Z_{i}Z_{j}(\forall i,j\in\Lambda$), we have the corresponding relevant Pauli strings

[TABLE]

Then we can construct the variational ansatz for the target Ising Hamiltonian, i.e., for finite $\alpha$, since the lattice sites are all-to-all coupled with the $Z_{i}Z_{j}$ term, the Pauli exponentials gates $e^{-i\theta Z_{i}Y_{j}}$ and $e^{-i\theta Y_{i}Z_{j}}$ are added with all-to-all connection. While for infinite $\alpha$, only nearest-neighbor couplings are considered in the Hamiltonian. So the gates are added only between pairs of nearest neighbor sites. An example of a variational ansatz for nearest-neighbor Ising chain under periodic boundary conditions is shown in figure 2c. Each layer of the variational ansatz consists of one layer of ZY-Pauli exponentials and one layer of YZ-Pauli exponentials, as shown in the dashed box. The Pauli exponentials are ladder-arranged, to increase the correlation between sites that can be generated by the variational ansatz. In the figure, we show the case of layers $L=1$. In the following numerical simulations, we use $L=2$ if not specified otherwise.

We assume the imaginary time evolution of each local interaction term $e^{\tau Z_{i}Z_{j}}$ can be realized with the Pauli exponentials $e^{-i\theta Z_{i}Y_{j}}e^{-i\theta^{\prime}Y_{i}Z_{j}}$, which have the same support of $Z_{i}Z_{j}$. These two Pauli exponentials are enough in the 2-qubit case as indicated by Eq. (34), but are not when the system size is large and when the system approaches the critical point, as explained in the previous section. It means that the expressivity of this variational ansatz is not sufficiently good to carry out the whole imaginary time evolution $e^{-\tau H}$. Limited expressivity leads to systematic errors, which will affect the numerical results.

First, we present the numerical results of the nearest-neighbor Ising model(NNIM), i.e., taking the limit $\alpha\rightarrow\infty$ in Eq. (1). With the nearest-neighbor interaction, there are $N=2D|\Lambda|L$ parameters in the variational ansatz. In two and three-dimensional NNIMs, there is a second-order phase transition in the infinite volume limit, where the critical points are $K_{c}=\ln(1+\sqrt{2})/2\approx 0.441$ Schultz_64 and $0.222$ Sonsin_2015 for dimension $D=2,3$, respectively. The specific heat and susceptibility would hence diverge near the critical point in the infinite volume limit. Figure 3 shows the specific heat and susceptibility for various $K$ values obtained via QITE-ansatz. The lattice size is $2\times 2,3\times 3,4\times 4$ for the 2-D system, marked by triangular-down, circle and triangular-up, respectively, and $2\times 2\times 2,3\times 3\times 2$ for the 3-D system, with results marked by triangular-down and circle respectively. In the evolution of the variational parameters, we use the Euler method with step length $\delta\tau=0.002$ as in Eq. (19), which is chosen such that further shrinking the step length has no impact on the numerical results (We will take this step length also throughout the following simulations.). We see that the QITE results converge well with the results from exact diagonalization(ED) when the system size is small for both 2-D and 3-D systems. For $4\times 4$ and $3\times 3\times 2$ lattices, the specific heat curves deviate from the ED curves near the critical point, which result from the limitation of the variational ansatz expressivity. The expressivity can be improved by increasing the number of ansatz layers and using longer Pauli strings for each local interaction term beyond $Z_{i}Y_{j},Y_{i}Z_{j}$. More detailed error analyses are shown in Appendix B.

Indications of the Ising criticality can be observed in figure 3. The critical temperatures of 2-D and 3-D systems in the infinite volume limit are denoted by the black dashed line. Near the critical points, values of the specific heat and susceptibility increase, and there are peaks in the specific heat as a function of $K$. For 2-D NNIM with volume $N_{d}\times N_{d}$, we denote the position of the peak as $K_{c}(N_{d})$. For larger system sizes, $K_{c}(N_{d})$ moves slowly towards the infinite volume critical point $K_{c}$. To guide the eye of this movement, we draw the grey solid line in the upper left panel of figure 3. The analytic expression of the grey solid line is inspired by the finite size scaling Pascal_lecture.

Figure 4 presents the behavior of the specific heat for the 2-D long-range interacting Ising model with finite $\alpha$. Compared with the nearest neighbor interaction, the long-range model introduces more $Z_{i}Z_{j}$ interactions and requires more variational parameters. There are $N=|\Lambda|(|\Lambda|-1)L$ parameters in the variational ansatz. The system size in the figure is $|\Lambda|=3\times 3$, with $\alpha=1,2,3$ and the nearest neighbor case $\alpha=\infty$, marked by the triangular-up, cross, triangular-down and circle, respectively. We see that for various $\alpha$ and $K$, the QITE-ansatz results and the ED results are consistent. Moreover, the peak of the specific heat shifts to the direction of high temperature(smaller $K$) for a larger interaction range(smaller $\alpha$). This behavior is reasonable since the long-range interaction effectively raises the system’s dimension, and a higher system dimension leads to a higher critical temperature, e.g., 3-D NNIM critical temperature is higher than that of 2-D NNIM.

V Discussion

This work discuss the possibility of using the imaginary time evolution algorithm to prepare the thermal state of the Ising model on NISQ devices. We numerically calculate the specific heat and susceptibility of the long-range interacting Ising model with the prepared thermal state. We find that the results using the quantum algorithm are consistent with the ones from exact diagonalization for various temperatures, including the critical and low-temperature regions.

We present a systematic procedure to design a variational ansatz for the thermal state preparation. This ansatz is inherited from the quantum circuits used in QITE-measure algorithm. We show that it out-performs the original circuit designed using QITE-measure. This variational ansatz can be further simplified according to the symmetries of the Hamiltonian and the initial state.

In our numerical simulation results, indication of critical behavior can be observed in the calculation of heat capacity and susceptibility of 2-D and 3-D Ising model. The universality properties of Ising model including the critical exponents can be extracted from these quantities in the thermal dynamic limit, where larger Ising system should be simulated to approach the limit. Larger Ising system simulations resort to more advanced quantum devices with more qubits and less error, which are hopefully to be experimentally realized in the near future.

The ideas proposed in this work can be applied to study the critical behavior of other classical models, such as the $Q$-state Potts model, which would be difficult to simulate using the Monte-Carlo algorithm when $Q$ is very large. Additionally, according to the correspondence of the $D$ dimensional quantum model to the $D+1$ dimensional classical model Hsieh2012FromDQ , the algorithm can also be used to study quantum phase transition.

VI Acknowledgements

We thank Xiao Yuan, Jinzhao Sun, Lena Funcke, Stefan Kühn and Yahui Chai for helpful discussions. X.W. and X.F. were supported in part by NSFC of China under Grants No. 12125501, No. 12070131001, and No. 12141501, and National Key Research and Development Program of China under No. 2020YFA0406400. PS acknowledges support from: ERC AdG NOQIA; Ministerio de Ciencia y Innovation Agencia Estatal de Investigaciones (PGC2018-097027-B-I00/10.13039/501100011033, CEX2019-000910-S/10.13039/501100011033, Plan National FIDEUA PID2019-106901GB-I00, FPI, QUANTERA MAQS PCI2019-111828-2, QUANTERA DYNAMITE PCI2022-132919, Proyectos de I+D+I “Retos Colaboración” QUSPIN RTC2019-007196-7); MICIIN with funding from European Union NextGenerationEU(PRTR-C17.I1) and by Generalitat de Catalunya; Fundació Cellex; Fundació Mir-Puig; Generalitat de Catalunya (European Social Fund FEDER and CERCA program, AGAUR Grant No. 2021 SGR 01452, QuantumCAT  U16-011424, co-funded by ERDF Operational Program of Catalonia 2014-2020); Barcelona Supercomputing Center MareNostrum (FI-2022-1-0042); EU Horizon 2020 FET-OPEN OPTOlogic (Grant No 899794); EU Horizon Europe Program (Grant Agreement 101080086 — NeQST), National Science Centre, Poland (Symfonia Grant No. 2016/20/W/ST4/00314); ICFO Internal “QuantumGaudi” project; European Union’s Horizon 2020 research and innovation program under the Marie-Skłodowska-Curie grant agreement No 101029393 (STREDCH) and No 847648 (“La Caixa” Junior Leaders fellowships ID100010434: LCF/BQ/PI19/11690013, LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012, LCF/BQ/PR21/11840013). Views and opinions expressed in this work are, however, those of the author(s) only and do not necessarily reflect those of the European Union, European Climate, Infrastructure and Environment Executive Agency (CINEA), nor any other granting authority. Neither the European Union nor any granting authority can be held responsible for them.

Appendix A Simplification of thermal state preparation in classical field theory

The Hamiltonian of a classical field theory is naturally diagonalized and can be written as a linear combination of Pauli-$Z$ operators, such as the Ising model considered in the main text and $Q$-state Potts model. Such Hamiltonian has energy eigenstates that can be encoded on the computational basis of qubits, and all the Pauli-$Z$ operators commute with each other. To compute the expectation values of such Hamiltonian’s thermal state, we only need imaginary time evolution on an initial state $\ket{\tilde{+}}\equiv\ket{+}^{\otimes N_{q}}$ where $N_{q}$ is the number of system’s qubits and $\ket{+}=(\ket{0}+\ket{1})/\sqrt{2}$. A similar idea is also proposed in the tensor network algorithm targeting on classical Ising model Verstraete_06 . The above statement is proved as follows.

The thermal expectation values $\langle O\rangle$ as defined in Eq. (4) can be expanded with an arbitrary orthogonal basis $\{\ket{i}\}$

[TABLE]

where

[TABLE]

We choose the orthogonal basis of Pauli-$X$ operators $\{\ket{i}\}=\{\ket{+},\ket{-}\}^{\otimes N_{q}}$. Notice that all vectors in the set can be generated by applying Pauli-$Z$ operators on one basis vector $\ket{\tilde{+}}$. For example

[TABLE]

The Hamiltonian consists of Pauli-$Z$ operators, so it commutes with all the Pauli-$Z$ operators. Thus, all terms in the partition function are equal

[TABLE]

for all $\ket{i}\in\{\ket{+},\ket{-}\}^{\otimes N_{q}}$, and we have $Z_{2\tau}=2^{N_{q}}\bra{\tilde{+}}e^{-2\tau H}\ket{\tilde{+}}$. Further, notice that all the observables concerning specific heat and susceptibility in Eq. (8) consist of Pauli-$Z$ operators, which can be formally written as

[TABLE]

where $\tilde{Z}_{m}$ denotes the tensor product of $Z$ operators at some sites and identity operators at others. Similar to Eq. (40), all terms in the numerator of Eq. (37) are equal

[TABLE]

for all $\ket{i}\in\{\ket{+},\ket{-}\}^{\otimes N_{q}}$. Thus we have

[TABLE]

In conclusion, the thermal expectation value of an observable $O=\sum_{m}h_{m}\tilde{Z}_{m}$ with a thermal state of a classical Hamiltonian can be derived according to imaginary time evolution on initial state $\ket{\tilde{+}}$,

[TABLE]

where $\ket{\tilde{+}(\tau)}$ is imaginary time evolved state according to Eq. (9). The state is initialized as $\ket{\tilde{+}(0)}=\ket{\tilde{+}}$. Thus we prove the statement in Eq. (26).

Appendix B Error analysis and circuit layers estimation

There are four main sources of errors when implementing the QITE-ansatz algorithm on real quantum devices LiYing2017

• •

The variational ansatz has limited expressivity. The imaginary time evolution proceeds on the manifold expanded by the variational ansatz. Thus the evolved wave function deviates from the true wave function in Eq. (9), and leads to the systematic error of the expectation values of the observables.

• •

Errors arise from the numerical integration using the Euler method as in Eq. (19).

• •

Noisy quantum gates and readout processes in quantum devices result in systematic errors when evaluating expectation values and estimating $M$ and $V$ (See Eq. (17)).

• •

Finite number of shots results in statistical errors in evaluating expectation values, $M$ and $V$.

Errors from the first and the second items are specific to the QITE-ansatz algorithm. The third and forth errors exist in general for any quantum algorithms. In the following contents, we will discuss these errors in detail.

The errors from the limited variational ansatz expressivity have been shortly discussed in the main text. There are two ways to improve expressivity. The first is by increasing the number of ansatz layers, and the second is by considering longer Pauli strings expansion in Eq. (29) for each local interaction term in the Hamiltonian. It is not hard to see that by extending number of layers to infinity and taking the expansion on the whole system, the variational ansatz can carry out the evolution $e^{-\tau H}$ exactly. In the following text, we numerically investigate how these two aspects affect the performance in calculating the specific heat of 2-D NNIM.

The limitation of finite ansatz layers can be observed by tuning the number of layers $L$. In figure 5, we compute the average absolute error of 2-D NNIM specific heat as a function of $L$, in case of lattice volumes $|\Lambda|=3\times 3,4\times 4$. The average absolute error is defined by

[TABLE]

where $C_{v}$ is specific heat from the quantum simulator, and $C_{v}^{ED}$ is from exact diagonalization. Here we take the integration range $[K_{\min},K_{\max}]=[0,1]$. The errors of specific heat decrease rapidly as $L$ increase and saturate to a platform after a certain layer $L^{*}$. We will analyze this transition layers $L^{*}$ after a while. When $L>L^{*}$, the remaining average absolute error of specific heat is mainly from the finite length of Pauli strings expansion.

Here provide an empirical explanation of the transition layer $L^{*}$, as observed in figure 5. It also helps to estimate how many layers we need when constructing variational ansatz for simulating NNIM. As shown in figure 6, variational ansatz generates correlation in the spin system. In the best case, the correlation between two neighboring spins is generated by one unitary transformation such as $e^{-i\theta ZY}$ in the Ising case; In the worst case, we need a whole layer of the variational ansatz such as $e^{-i\theta ZY}e^{-i\theta^{\prime}YZ}$ to generate such correlation. The transition layer $L^{*}$ indicates the lowest number of circuit layers to generate correlation between the two most distant spins in the D-dimensional nearest neighbor lattice system. Thus for D-dimensional NNIM with volume $N_{d}^{D}$ and PBC, as the Manhatten distance of the most remote two spins is $DN_{d}/2$ (Equal to the number of the yellow arrows in figure 6, where $D=2,N_{d}=3,4$ respectively.), the transition layer would be in the range

[TABLE]

which corresponds to the best case and worst case mentioned above. Here $G$ is the number of Pauli exponentials in one layer, i.e., the number of relevant Pauli operators for some local interaction terms. The transition layers in figure 5 are in accord with this range, i.e., $N_{d}/2\leq L^{*}\leq N_{d}$, and we see larger number of layers have almost no improvement to the average absolute error of the specific heat. Thus we say $L^{*}$ layers are enough for variational ansatz to simulate NNIM. This estimation on the number of ansatz layers can be generalized to more complicated short-range interacting models.

Comparing the required number of layers of the variational ansatz provided by Eq. (46) and the layers of quantum circuits used in QITE-measure, one finds the former is much less than the latter. It can be partially explained using the example of the two-qubit Ising system shown in the main text. For the QITE-measure circuit figure 2a, due to the commutability of relevant Pauli operators $[ZY,YZ]=0$, it is equivalent to the circuit shown in figure 7, which consists of only two Pauli exponentials where the rotation angles are the summation of all the coefficients of the corresponding Pauli exponentials in figure 2a. Therefore, if we use one layer of the circuit in figure 2b, and $\theta_{1}=a_{1}^{(0)}+\ldots+a_{1}^{(\tau-\Delta\tau)},\theta_{1}^{\prime}=a_{2}^{(0)}+\ldots+a_{2}^{(\tau-\Delta\tau)}$, the QITE-measure circuit could be rephrased without loss of the precision. Thus compared with the QITE-measure circuit, the number of variational ansatz layers used in our simulation can be significantly reduced.

Numerical integration errors can be controlled via a more elaborate numerical integration algorithm. In the main text, we use the Euler method that accumulates a global error of $\mathcal{O}(\delta\tau)$ at the final step. One could use a more elaborate numerical algorithm such as the 4th-order Runge-Kutta method to control the systematic error, which accumulates a global error of $\mathcal{O}(\delta\tau^{4})$ at the final step. In our simulations, as the numerical integration error is not the dominate systematic error, the Euler method is sufficiently good.

Errors from noisy quantum gates and readout processes lead to systematic deviations of the measurement results to the noiseless ones. For NISQ devices, there are many error mitigation techniques to reduce these deviations. For example, errors from two-qubit gates can be mitigated by zero-noise extrapolation Temme_2017 ; LiYing2017 and quasi-probability decomposition LiYing2017 ; Suguru2018 . The readout error can be mitigated by classical bit-flip correction Lena2022 ; Berg2022 . The error mitigation techniques reduce the systematic deviations of the noisy results to the noiseless ones and shed light on the real applications of NISQ devices Kim2023 ; Cao2023 .

Finite number of shots error is a statistical error, which can be suppressed by increasing the number of shots. In the measurement procedure, the observables are split into a weighted sum of Pauli operators and each can be measured separately, at the cost of many shots. The number of shots can be reduced by collecting mutually commuting Pauli operators together before measuring all operators within a collection simultaneously Kandala2017 ; Crawford2021efficientquantum . In the measurement process of the Ising model in Figure 1, the weighted Pauli operators are mutually commuting Pauli-$Z$ strings. They can be measured simultaneously in computational basis.

Appendix C Execution time of the nearest-neighbor Ising model

In this appendix, we estimate the execution time to study the critical behavior of the nearest-neighbor Ising model (NNIM) using the QITE-ansatz algorithm.

The execution time of the QITE-ansatz algorithm can be estimated by the number of steps of the imaginary time evolution, times the number of expectation values evaluated per step, times the number of two-qubit quantum gates (Assume no parallelization of the quantum gates, and the two-qubit gates dominate the execution time), i.e.,

[TABLE]

Among the three factors, if the critical temperature of the system is of $\mathcal{O}(1)$, then the total evolution should have the same order, and the number of evolution steps is a constant, which is the case of simulating NNIM.

For each step of the time evolution, one needs to estimate the expectation values of entries of the $M$ matrix and $V$ vector. Thus

[TABLE]

where $N$ is the dimension of the $M$ matrix, equals the number of variational parameters. For the variational ansatz proposed in this work, the number of variational parameters is $N=2D|\Lambda|L$ for NNIM, where $D$ is the spatial dimension of the NNIM lattice $\Lambda$, $|\Lambda|$ is the lattice size, $L$ is the number of ansatz layers. In Appendix B of the manuscript, we estimate the number of layers required for NNIM, where $L\sim\mathcal{O}(DN_{d})$, $N_{d}$ is the side length of $\Lambda$. The lattice size and the side length are related by $|\Lambda|=N_{d}^{D}$. Thus, the number of expectation evaluations reads

[TABLE]

To count the number of two-qubit quantum gates, we assume that the CNOT gate is the basic two-qubit gate that can be realized on quantum devices, and the CNOT gate can be applied to every two-qubit pair. Then, the number of CNOT gates in the proposed ansatz is proportional to its number of variational parameters $N$, so that

[TABLE]

In summary, the execution time of the QITE-ansatz algorithm for D-dimensional NNIM reads

[TABLE]

which is a polynomial function of the system size.

There are some improvement methods of the QITE-ansatz algorithm to reduce the execution time, such as the DualQITE algorithm proposed in reference gacon2023variational . Using this algorithm, one can reduce the number of expectation evaluations in Eq. (48) from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N)$, and the total execution time is correspondingly reduced to

[TABLE]