Generating approximate ground states of strongly correlated quantum many-body systems through quantum imaginary time evolution

We present an overview of the QITE workflow in figure 1. One first considers the formalized problem introduced in section 2.1, a Hamiltonian $\hat{H}$ that is composed of $\tilde{m}$ terms $\hat{h}\left[j\right]$ with $j\in [\tilde{m}]$. The goal of QITE is to imitate the ITE, which requires the preparation of an initial state ∣Ψ_init〉 which is discussed in section 4. The term $\hat{h}\left[j\right]$ can either describe a linear combinations of Pauli strings (see section 2.3.3) or a linear combination of fermionic operators (see section 2.3.4). Depending on whether the underlying Hamiltonian describes a lattice system or a molecular electronic structure problem, we perform an operator expansion on a set D_j of spins, fermionic sites or orbitals. The choice of D_j is crucial for the success of QITE and is guided by the Manhattan distance when we deal with a lattice system (see Section IIC2 a), and by the mutual information between orbitals when we investigate molecular electronic structure systems (see section 2.4), respectively. Finally, section 2.3.2 gives a short introduction to the main idea of the QITE algorithm, with a more detailed discussion on how this is put to practice in spin and fermionic systems in section 2.3.3 and 2.3.4, respectively.

Zoom In Zoom Out Reset image size

Code details We implemented this workflow (see figure 1) using Python. Our calculations rely on state vector simulations of the QITE algorithm on classical hardware. For our initial state preparation, we calculated the Pfaffian of skew-symmetric matrices using the pfapack library[34, 35]. Molecular electronic structure Hamiltonians for different basis sets as well as the studied Fermi-Hubbard system and their corresponding approximate ground states were determined using PySCF 2.2.1. [36–38]. The OpenFermion 1.5.1 [39] library was used for generating the spin and fermionic operators, as well as to translate the obtained covariance matrix of a pure fermionic Gaussian state into a state vector in its Hilbert space representation through the relations we derive in Appendix A.4. The system of linear equations to derive the operator expansion has been solved by a Conjugate Gradient iteration of the SciPy 1.10.1. [40] library. Our calculations of the mutual information were supported by the TeNPy 0.10.0 [41] library, which was used to compute the required reduced density matrices of the ground state more efficiently, by using the Matrix Product State (MPS) representation of the latter for a large bond dimension.

The main idea of the QITE algorithm is to approximate the trotterized imaginary time propagator ${e}^{-{\rm{\Delta }}\tau \hat{h}[f(l)]}$ by a unitary operator ${\hat{U}}_{l}$ at each iteration step, without the need of auxiliary qubits. More specifically, for a fixed discretization step size Δτ and step l ∈ [mn], the QITE algorithm tries to find the unitary ${\hat{U}}_{l}$ which approximates

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{l+1}\rangle \equiv \frac{1}{\sqrt{{c}_{l}}}{e}^{-{\rm{\Delta }}\tau \hat{h}[f(l)]}| {{\rm{\Psi }}}_{l}\rangle \approx {\hat{U}}_{l}| {{\rm{\Psi }}}_{l}\rangle ,\end{eqnarray} \tag{ 8 }$$

by minimizing the distance

$$\begin{eqnarray}&&\mathop{\min }\limits_{{U}_{l}}\unicode{x02016}\frac{1}{\sqrt{{c}_{l}}}{e}^{-{\rm{\Delta }}\tau \hat{h}[f(l)]}| {{\rm{\Psi }}}_{l}\rangle -{\hat{U}}_{l}| {{\rm{\Psi }}}_{l}\rangle \unicode{x02016},\end{eqnarray} \tag{ 9 }$$

where $\parallel | \varphi \rangle \parallel =\sqrt{\langle \varphi | \varphi \rangle }$ is the Euclidean norm of a vector ∣φ〉 and

$$\begin{eqnarray}&&{c}_{l}=\left|\langle {{\rm{\Psi }}}_{l}| {e}^{-2{\rm{\Delta }}\tau \hat{h}[f(l)]}| {{\rm{\Psi }}}_{l}\rangle \right|\end{eqnarray} \tag{ 10 }$$

is a normalization constant. The algorithm is initialized in the state ∣Ψ_l=0〉 = ∣Ψ_init〉 which must have non-zero support on the desired ground state, γ_init > 0. In the following, we will describe how the unitary U_l can be found at each iteration step l, following [6].

Any unitary operator can be written as

$$\begin{eqnarray}&&{\hat{U}}_{l}={e}^{-i{\rm{\Delta }}\tau \hat{A}[l]},\end{eqnarray} \tag{ 11 }$$

where $\hat{A}[l]$ is a Hermitian operator. Therefore, the task of finding ${\hat{U}}_{l}$ is equivalent to that of finding $\hat{A}[l]$. We will assume sufficiently small time steps, such that ${\hat{U}}_{l}={\bf{1}}-i{\rm{\Delta }}\tau \hat{A}[l]+{ \mathcal O }({\rm{\Delta }}{\tau }^{2})$. While the computational cost of finding a general unitary operator acting on L qubits or fermionic modes will scale as ${ \mathcal O }(\exp (L))$, the central finding of [6] is that for lattice systems with a finite correlation length ${ \mathcal C }$, the cost of finding the unitary ${\hat{U}}_{l}$ as well as approximating this unitary in terms of elementary quantum gates scales only as ${ \mathcal O }(\exp ({ \mathcal C }))$. Importantly, for systems with a finite correlation length, while still exponential, this is a finite value which does not increase as one increases the size of the system.

a. Support and domain We define the supportS_j of an operator $\hat{h}[j]$ as the set of indices it acts on non-trivially (i.e. not as the identity), where $j\in [\tilde{m}]$. For spin operators these indices correspond to the spin indices, while for fermionic systems they correspond to fermionic sites or orbitals. The domainD_j similarly defines a set of indices which includes the support S_j, but also includes other indices that are selected based on the specifics of the studied system. In order to describe which indices to include in the domain D_j, one introduces a parameter ν which sets the size of the set (for d-dimensional lattices, D_j ∝ ν^d, while for molecular systems ν will set the number of additional orbitals added to the support). Depending on the studied system, the indices will correspond to spin or fermionic sites, or orbitals.

In figure 2 we give examples of how we choose the support and domain when $\hat{h}[j]$ describes a spin- or fermionic lattice Hamiltonian, or a molecular electronic structure Hamiltonian. Figure 2(a) describes a term $\hat{h}[j]$ which is part of a 18 spins (i.e. 18 qubits) Hamiltonian whose underlying geometry describes a hexagonal spin lattice. For simplicity, we here assume that there are no periodic boundary conditions, contrary to the systems studied later on. We further assume that the interaction strength decreases with increasing distance between two sites, which justifies using the Manhattan distance as a means of picking the domain. In the provided example, the support S_j = {8, 14} (highlighted in yellow) is of size ∣S_j∣ = 2 and the domain is set by the Manhattan distance ν = 1 (ν = 2). In this case, the domain D_j is given by all yellow and purple (and pink) colours, leading to a domain ∣D_j∣ = 7 (∣D_j∣ = 13). Figure 2(b) describes an example where the term $\hat{h}[j]$ is part of a Hamiltonian describing a fermionic lattice of L = 8 sites with short-range interactions, each site containing two possible spin configurations α or β, thus a total of 16 fermionic modes or qubits. Here, the support of $\hat{h}[j]$ contains sites S_j = {3, 4}, which corresponds to 4 fermionic modes or qubits. For the fermionic lattice, ν describes the number of additional sites added, where each site contains two fermionic modes. The domain for ν = 1 contains three sites, but it is not clear whether to pick site 2 or site 5. In this case, we pick one of the two randomly, thus a potential outcome could be D_j = {3, 4, 5}. For ν = 2 (ν = 4), the resulting domain would then include all yellow and purple (and pink) sites.

Zoom In Zoom Out Reset image size

In figure 2(c), $\hat{h}[j]$ is part of a molecular electronic structure Hamiltonian of L = 6 orbitals, i.e. 12 spin-orbitals or qubits. In this example, the support is given by S_j = {1, 5} and the domain is D_j = {1, 3, 5} (D_j = {0, 1, 3, 5}) for ν = 1 (ν = 2). Here, ν describes the number of additional orbitals besides those included in S_j and we choose the additional orbitals based on those which possess the largest mutual information I(p, q) (indicated by the thickness of the connecting lines), where p ∈ S_j and q ∈ [L]\S_j. Note, that the definition and particular form of the mutual information for the studied systems will be detailed in section 2.4 and Appendix G.

In general, the size of the support can vary and is bounded by the locality of the Hamiltonian. For systems with long-range (or also just beyond nearest-neighbor-) interactions, the size of the domain D_j of a given term $\hat{h}[j]$ will be larger than for systems restricted to only nearest-neighbor interactions. As we will see in the following, the computational cost of QITE will scale exponentially in the size of the domain D_j, and is the main limiting factor for the simulations we were able to run (and not the system size L).

b. Exact QITE It was shown in [6], that for d-dimensional lattice spin systems whose Hamiltonian is k-local, the support of the Hermitian operator $\hat{A}[l]$ in equation (11) can be chosen to be the domain D_l of the support S_l of the Hamiltonian term $\hat{h}[l]$ within a Manhattan distance ν for $l\in [\tilde{m}]$. More specifically, it was shown that for spin systems with a finite correlation length ${ \mathcal C }$, in order to approximate the ITE state within precision ε > 0, i.e.

$$\begin{eqnarray}&&\unicode{x02016}\displaystyle \prod _{l=1}^{mn}\frac{1}{\sqrt{{c}_{l}}}{e}^{-{\rm{\Delta }}\tau \hat{h}[f(l)]}| {{\rm{\Psi }}}_{\,\rm{init}\,}\rangle -\displaystyle \prod _{l=1}^{mn}{\hat{U}}_{l}| {{\rm{\Psi }}}_{\,\rm{init}\,}\rangle \unicode{x02016}\,\leqslant \,\varepsilon ,\end{eqnarray} \tag{ 12 }$$

one requires the Manhattan distance to be at least

$$\begin{eqnarray}&&{\nu }_{{ \mathcal C }}=2{ \mathcal C }\mathrm{ln}(2\sqrt{2}nm/\varepsilon ).\end{eqnarray} \tag{ 13 }$$

The proof of the above is based on Uhlmann’s theorem [42] and will not be subject of this work (a detailed proof is given in [6]). For a d-dimensional lattice spin system, the size of the domain of $\hat{h}[f(l)]$ scales as ∣D_f(l)∣ ≤ kν^d for every Hamiltonian term $\hat{h}[l]$. Using equation (13), this means that for a finite correlation length ${ \mathcal C }$, the number of indices in D(l) is independent of the system size L.

Since $\hat{A}[l]$ is unknown, one has to expand it in a basis of operators and determine its components. While in general such a basis would contain a number ${ \mathcal O }(\exp (L))$ of basis operators, a direct consequence of equation (13) is that the basis must only include roughly $\exp (k{\nu }^{d})$ basis operators, which is independent of the system size L. The running time T of the corresponding QITE algorithm for an imaginary time τ = nΔτ is then roughly given by

$$\begin{eqnarray}&&T=mn{e}^{k{\nu }^{d}}.\end{eqnarray} \tag{ 14 }$$

The above results provide the mathematical basis that legitimizes the Ansatz to approximate the non-unitary QITE by a sequence of unitary operators for lattice spin systems with a finite correlation length. It is then conjectured in [6] that a similar result holds for fermionic systems under similar conditions. However, many systems of interests display long-range interactions, diverging correlation lengths, or do not admit to be formulated as a simple lattice theory, such as molecular electronic structure Hamiltonians. Part of the objective of this work is to study how the QITE algorithm used as a heuristic Ansatz described in the following paragraph performs for such systems.

c. Main approximation to speed up ITE and QITE simulation In order to speed up the trotterized ITE, spin- and fermionic-QITE simulations, we are uniting all Hamiltonian terms $\hat{h}[l]$ which possess the same domain D_l for a given ν. However, if there are two domains that are subsets of each other, e.g. D_i = {1, 3} and D_j = {1, 2, 3}, we treat them separately. This reduces the number of Hamiltonian terms $\tilde{m}$ in equation (1) greatly, but also leads to terms $\hat{h}[l]$ which contain a linear combination of different spin and fermionic excitation operators that will be difficult to realize in practice.

d. Inexact QITE It is also possible to use the QITE algorithm as a heuristic algorithm to find an approximate ground state ∣Ψ〉 using a very small domain size, i.e. ν = 1, 2, … , which in many instances will be much smaller than the true correlation length ${ \mathcal C }$ of the studied system. In this case, QITE loses its rigorous performance guarantees from equations (12)–(13). This means that at some point in time, the unitary approximation ${\hat{U}}_{l}$ of the trotterized imaginary time propagator ${e}^{-{\rm{\Delta }}\tau \hat{h}[f(l)]}$ will start to deviate significantly and the resulting state ∣Ψ_l〉 will start to differ largely from its ITE counterpart. Nevertheless, it is still possible, that for small evolution times the evolved state will lead to an improvement over the classical input state ∣Ψ_init〉. We will examine this claim thoroughly in section 5.

If the Hermitian Hamiltonian term $\hat{h}[f(l)]$ describes a spin operator, the Hermitian operator $\hat{A}[l]$ can be determined from an expansion in a spin operator basis, as we explain in the following and refer to as spin-QITE. Let D_l denote the domain of $\hat{h}[l]$ which we assume to be a tensor product of spin operators ${\hat{\sigma }}_{i}^{x}$, ${\hat{\sigma }}_{i}^{y}$, and ${\hat{\sigma }}_{i}^{z}$, and we define ${\hat{\sigma }}_{i}^{0}={\mathbb{1}}$ as the identity on spin index i ∈ [L]. We write ${\{{\hat{\sigma }}^{0},{\hat{\sigma }}^{x},{\hat{\sigma }}^{y},{\hat{\sigma }}^{z}\}}^{\otimes {D}_{l}}$ in order to denote all possible ${4}^{| {D}_{l}| }$ Pauli strings (including the identity) on the indices contained in D_l. A single element of ${\{0,x,y,z\}}^{\otimes {D}_{l}}$ is given by ${\hat{\sigma }}_{{\bf{I}}}={\hat{\sigma }}_{{i}_{1}}^{{\alpha }_{1}}{\hat{\sigma }}_{{i}_{2}}^{{\alpha }_{2}}\cdots {\hat{\sigma }}_{{i}_{| D(l)| }}^{{\alpha }_{| D(l)| }}$, where I = {(i₁, α₁), (i₂, α₂), …, (i_∣D(l)∣, α_∣D(l)∣)}, 1≤i₁ < i₂ < ⋯ < i_∣D(l)∣<L and α_k ∈ {0, x, y, z}. We then define the expansion of an Hermitian operator $\hat{A}[l]$ in the spin basis as

$$\begin{eqnarray}&&\hat{A}[l]=\displaystyle \sum _{{\bf{I}}\in {\{0,x,y,z\}}^{\otimes {D}_{l}}}a{[l]}_{{\bf{I}}}{\hat{\sigma }}_{{\bf{I}}},\end{eqnarray} \tag{ 15 }$$

where $a{[l]}_{{\bf{I}}}\in {\mathbb{R}}$ are real-valued coefficients and ${\hat{\sigma }}_{{\bf{I}}}$ are Hermitian operators. Note, that in principle not all possible Pauli operator strings will contribute to the solution, and their number can be reduced by symmetry considerations and properties of the desired wave function [6]. However, in our spin system simulations we included all ${4}^{| {D}_{l}| }$ operators and leave their potential reduction for future work.

As shown in appendix B, the coefficients a[l]_I (and therefore also $\hat{A}[l]$ and ${\hat{U}}_{l}$), condensed in the vector a[l], can be obtained from a classical solution of the following system of linear equations,

$$\begin{eqnarray}&&{\bf{S}}[l]{\bf{a}}[l]=-{\bf{b}}[l],\end{eqnarray} \tag{ 16 }$$

where we defined the Gram matrix S[l] and vector b[l] with real-valued entries as

$$\begin{eqnarray}&&S{[l]}_{{\bf{I}},{\bf{J}}}=\langle {{\rm{\Psi }}}_{l-1}| \{{\hat{\sigma }}_{{\bf{I}}},{\hat{\sigma }}_{{\bf{J}}}\}| {{\rm{\Psi }}}_{l-1}\rangle ,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&b{[l]}_{{\bf{I}}}=i\langle {{\rm{\Psi }}}_{l-1}| [{\hat{\sigma }}_{{\bf{I}}},\hat{h}[l]]| {{\rm{\Psi }}}_{l-1}\rangle ,\end{eqnarray} \tag{ 18 }$$

which correspond to the real- and imaginary parts of the expectation values of $2{\hat{\sigma }}_{{\bf{I}}}{\hat{\sigma }}_{{\bf{J}}}$ and $2{\hat{\sigma }}_{{\bf{I}}}\hat{h}[l]$, respectively. We introduced the notation $[\hat{A},\hat{B}]=\hat{A}\hat{B}-\hat{B}\hat{A}$ for the commutator, and $\{\hat{A},\hat{B}\}=\hat{A}\hat{B}+\hat{B}\hat{A}$ for the anti-commutator of two operators $\hat{A},\hat{B}$. When considering the Gram matrix S[l], we observe that ${\hat{\sigma }}_{{\bf{I}}}{\hat{\sigma }}_{{\bf{J}}}=\alpha ({\bf{I}},{\bf{J}}){\hat{\sigma }}_{{\bf{K}}}$, where α(I, J) ∈ { ± 1, ± i} is a phase factor that depends on the Pauli strings described by the index vectors I, J. Thus, even though the Gram matrix is a ${4}^{| {D}_{l}| }\times {4}^{| {D}_{l}| }$ matrix, it only requires the evaluation of ${4}^{| {D}_{l}| }$ expectation values $\langle {{\rm{\Psi }}}_{l-1}| {\hat{\sigma }}_{{\bf{K}}}| {{\rm{\Psi }}}_{l-1}\rangle$.

At each step of the QITE algorithm, the wave function ∣Ψ_l〉 has to be generated, and the ${4}^{| {D}_{l}| }$ Pauli string operator expectation values in equations (17)–(18) have to be measured on the quantum computer. As previously stated, in practice one still has to decompose each unitary ${\hat{U}}_{l}$ into e.g. one- and two-qubit gate operations, a cost that will also scale as ${ \mathcal O }(\exp | {D}_{l}| )$ [43]. We will ignore this cost (as well as the potential error that comes with its implementation) entirely in this work and will simply assume that U_l can be implemented exactly. We also ignore the error resulting from using a finite number of measurements in order to evaluate equations (17)–(18). This means that the actual performance of the Spin-QITE (as well as the fermionic-QITE formulation introduced in the following section) in an actual quantum computation will most likely perform worse than our algorithmic studies would indicate.

If the Hermitian Hamiltonian term $\hat{h}[f(l)]$ describes a fermionic operator, the Hermitian operator $\hat{A}[l]$ can be determined from an expansion in a fermionic operator basis, as we explain in the following and refer to as fermionic-QITE. We expand the Hermitian operator $\hat{A}[l]$ in terms of products of Majorana operators ${\hat{\gamma }}_{\mu }$ with $\mu \in [2{L}^{{\prime} }]$, where ${L}^{{\prime} }$ denotes the actual number of fermionic modes (for instance, for a Fermi-Hubbard lattice with L sites or a molecular electronic structure Hamiltonian described by a basis of L orbitals, ${L}^{{\prime} }=2L$). Majorana operators are Hermitian operators that fulfill the anticommutator relation $\{{\hat{\gamma }}_{\mu },{\hat{\gamma }}_{\nu }\}=2{\delta }_{\mu \nu }$ and can be expressed in terms of fermionic creation and annihilation operators through ${\hat{\gamma }}_{2j-1}={\hat{c}}_{j}+{\hat{c}}_{j}^{\dagger }$ and ${\hat{\gamma }}_{2j}=-i({\hat{c}}_{j}-{\hat{c}}_{j}^{\dagger })$ for $j\in [{L}^{{\prime} }]$. We further define an ordered set of indices μ = {μ₁, μ₂, …, μ_∣μ∣}, where ${\mu }_{k}\in [2{L}^{{\prime} }]$ and $0\leqslant {\mu }_{1}\lt {\mu }_{2}\lt \cdots \lt {\mu }_{| {\boldsymbol{\mu }}| }\leqslant 2{L}^{{\prime} }$. Since we will only consider systems that conserve parity, we restrict the number of elements in our sets to be even, ∣μ∣ = 2m, where $m\subseteq [{L}^{{\prime} }]$. We introduce the following short-hand notation for an ordered product of Majorana operators

$$\begin{eqnarray}&&{\hat{\gamma }}_{{\boldsymbol{\mu }}}={(-i)}^{\left(\genfrac{}{}{0.0pt}{}{| {\boldsymbol{\mu }}| }{2}\right)}\displaystyle \prod _{\mu \in {\boldsymbol{\mu }}}{\hat{\gamma }}_{\mu },\end{eqnarray} \tag{ 19 }$$

where the prefactor ensures that ${\hat{\gamma }}_{{\boldsymbol{\mu }}}$ is hermitian. Similar to equation (15), for fermionic systems we then choose the following basis expansion,

$$\begin{eqnarray}&&\hat{A}[l]=\displaystyle \sum _{{\boldsymbol{\mu }}}a{[l]}_{{\boldsymbol{\mu }}}{\hat{\gamma }}_{{\boldsymbol{\mu }}}.\end{eqnarray} \tag{ 20 }$$

Here, the sum goes over various sets μ which have yet not been specified (except for them being ordered and containing an even number of elements from $[2{L}^{{\prime} }]$). By choosing $a{[l]}_{{\boldsymbol{\mu }}}\in {\mathbb{R}}$, the operator $\hat{A}[l]$ will be hermitian, and we can use the same linear equation that we derived for the Spin-QITE, equation (B9),

$$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{\mu }}}a{[l]}_{{\boldsymbol{\mu }}}\langle {{\rm{\Psi }}}_{l-1}| \{{\hat{\gamma }}_{{\boldsymbol{\mu }}},{\hat{\gamma }}_{{\boldsymbol{\nu }}}\}| {{\rm{\Psi }}}_{l-1}\rangle =-i\langle {{\rm{\Psi }}}_{l-1}| [{\hat{\gamma }}_{{\boldsymbol{\nu }}},\hat{h}[l]]| {{\rm{\Psi }}}_{l-1}\rangle ,\end{eqnarray} \tag{ 21 }$$

or equivalently

$$\begin{eqnarray}&&{\bf{S}}[l]{\bf{a}}[l]=-{\bf{b}}[l],\end{eqnarray} \tag{ 22 }$$

where the minus sign is convention and we defined the following real-valued entries

$$\begin{eqnarray}&&S{[l]}_{{\boldsymbol{\mu ,\nu }}}=\langle {{\rm{\Psi }}}_{l-1}| \{{\hat{\gamma }}_{{\boldsymbol{\mu }}},{\hat{\gamma }}_{{\boldsymbol{\nu }}}\}| {{\rm{\Psi }}}_{l-1}\rangle ,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&b{[l]}_{{\boldsymbol{\nu }}}=i\langle {{\rm{\Psi }}}_{l-1}| [{\hat{\gamma }}_{{\boldsymbol{\nu }}},\hat{h}[l]]| {{\rm{\Psi }}}_{l-1}\rangle .\end{eqnarray} \tag{ 24 }$$

Since in QITE we expand the operator $\hat{A}[l]$ in the domain of the current Hermitian Hamiltonian term $\hat{h}[l]$, we will define a corresponding set of indices ${D}_{l}^{(\gamma )}$: These are the Majorana indices of the domain D_l. To illustrate, consider the case where the domain contains two indices D_l = {i, j} with $i\lt j\in [{L}^{{\prime} }]$, then ${D}_{l}^{(\gamma )}=\{2i-1,2i,2j-1,2j\}$. Therefore, the sum in equation (20) goes over all ${\boldsymbol{\mu }}\subseteq {D}_{l}^{(\gamma )}$ in the most general case. If one were to include the identity operator 1, one could express every Hermitian operator that acts in the Hilbert-subspace containing D_l modes with equation (20). However, it is often possible to truncate the expansion at lower order polynomials ∣μ∣ < 2∣D_l∣ and use basis expansion operators which respect the symmetry of the underlying system Hamiltonian, similar to the procedure in (unitary) coupled cluster approaches [44–46]. In what follows, we restrict ourselves to expansions of up to order ∣μ∣ = 4, but in principle, higher order expansions are also possible.

Number-conserving Ansatz Since for many fermionic systems, the particle number operator $\hat{N}={\sum }_{p=1}^{L}\left({\hat{c}}_{p,\alpha }^{\dagger }{\hat{c}}_{p,\alpha }+{\hat{c}}_{p,\beta }^{\dagger }{\hat{c}}_{p,\beta }\right)$ and other operators, such as the spin projection operator ${\hat{S}}_{z}=\frac{1}{2}{\sum }_{p=1}^{L}\left({\hat{c}}_{p,\alpha }^{\dagger }{\hat{c}}_{p,\alpha }-{\hat{c}}_{p,\beta }^{\dagger }{\hat{c}}_{p,\beta }\right)$ commute with the system Hamiltonian, one can reduce the number of terms in equation (20) by restricting the operator basis ${\hat{\gamma }}_{{\boldsymbol{\mu }}}$ to operators that respect those symmetries. We define the orbital excitation operator

$$\begin{eqnarray}&&{\hat{E}}_{pq}={\hat{c}}_{p,\alpha }^{\dagger }{\hat{c}}_{q,\alpha }+{\hat{c}}_{p,\beta }^{\dagger }{\hat{c}}_{q,\beta }.\end{eqnarray} \tag{ 25 }$$

This allows us to define the single orbital excitation operator

$$\begin{eqnarray}&&i({\hat{E}}_{pq}-{\hat{E}}_{pq}^{\dagger })\end{eqnarray} \tag{ 26 }$$

and double excitation operator

$$\begin{eqnarray}&&i({\hat{E}}_{pq}{\hat{E}}_{rs}-{\hat{E}}_{rs}^{\dagger }{\hat{E}}_{pq}^{\dagger }),\end{eqnarray} \tag{ 27 }$$

where we assume p < r and q < s. Equations (26)–(27) will replace the more general Ansatz γ_μ of equation (19) in our fermionic-QITE algorithm for both, fermionic lattice and molecular electronic structure Hamiltonians.

In principle, one would have to include not just the single and double excitation operators in equations (26)–(27) in γ_μ, but also up to n_el-th order excitation operators, where n_el denotes the number of particles, in order to be able to exactly reproduce the ITE behavior with QITE. Thus the single- and double restriction introduces an additional error, which is independent of the Trotterization error.

Note, that while fermionic-QITE can be formulated in the fermionic picture and then transformed into Pauli operators by means of e.g. the Jordan-Wigner transformation, one could directly implement it on a fermionic quantum processing unit [47], provided sufficient coherent control and the means of decomposing the unitaries ${\hat{U}}_{l}$ into the basic underlying fermionic operations is possible at acceptable cost.

We summarize the main steps of the trotterized ITE (steps 1+4a), exact ITE (steps 1+4b), spin-QITE (steps 1+2+3+5a) and fermionic-QITE (steps 1+2+3+5b), as well as the main sources of approximation errors (highlighted by red boxes) in figure 3.

Zoom In Zoom Out Reset image size

The Heisenberg model describes an effective model for Mott insulators, and follows as a special case of the Fermi-Hubbard Hamiltonian—described below in section 3.1.3—in the large U / t limit. The Heisenberg Hamiltonian is given by [62]

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{i\lt j}{J}_{ij}\left({\hat{\sigma }}_{i}^{x}{\hat{\sigma }}_{j}^{x}+{\hat{\sigma }}_{i}^{y}{\hat{\sigma }}_{j}^{y}+{\hat{\sigma }}_{i}^{z}{\hat{\sigma }}_{j}^{z}\right)+B\displaystyle \sum _{i}{\hat{\sigma }}_{i}^{z}.\end{eqnarray} \tag{ 29 }$$

Here, J_ij > 0 describes antiferromagnetic (AFM) coupling and B the strength of the external magnetic field. We implicitly assume an enumeration i = 1, 2, …, L of the underlying lattice. Equation (29) is not restricted to only describing one-dimensional, but can also describe higher dimensional lattice systems. For NN (NNN) interactions, J_ij ≠ 0 only when the indices i, j are next to each other in the lattice (two nearest-neighbor hops apart), which is indicated by 〈ij〉 (〈〈ij〉〉). Then, J_ij reduces to J_〈ij〉 = J. The Hamiltonian which describes a Heisenberg model restricted to J₁ = J_〈ij〉 and J₂ = J_{〈〈ij〉〉} interactions will also be studied, and is known as the J₁J₂-model [60]. We also consider long-range interactions, which have been realized in experimental setups [63] and are described by

$$\begin{eqnarray}&&{J}_{ij}=\frac{1}{| i-j{| }^{\alpha }},\end{eqnarray} \tag{ 30 }$$

where α is a free parameter which allows one to tune the system from a uniform (α = 0) to a NN (α → ∞) interaction. The range 0 < α < 1 in one dimensional systems is considered as the strong-long-range regime. Here, as well as in the TFIM below, we assume that ∣i − j∣ denotes the Manhattan distance of the sites i and j in the lattice (thus the minimal amount of nearest-neighbor hops in the lattice connecting two sites rather than the actual physical distance).

Chromium Triiodide While ferromagnetism in two-dimensional van der Waals (vdW) materials has been known for decades [64, 65], reports of isolated, atomically thin magnetic sheets through successful exfoliation from the bulk have only recently emerged [66], such as Cr₂Ge₂Te₆ [66], FePS₃ [67], VSe₂ [68], and MnSe [69]. One of the earliest and most studied 2D magnets is CrI₃ [70], which exhibits a honeycomb lattice where the magnetic Cr³⁺ ions within each layer are interconnected by bridging iodine ligands [71, 72]. The implication of integrating such magnets in vertical vdW heterostructures are tremendous, allowing controlled tuning of magnetic properties through tailored interlayer interactions or external fields to explore different magnetic order or to investigate new physical phenomena. A successful deployment in novel devices could potentially unlock their applications in spintronics where the electron spin is exploited as a further degree of freedom in addition to their charge, for example to store or transfer data, or for computation. To this end, spintronic devices have been explored extensively, including, e.g., spin valves and spin filters [73, 74].

From a theoretical perspective, while the Mermin-Wagner theorem excludes order in two dimensions at finite temperatures, the introduction of any anisotropy can lead to a gapping of the low-energy modes and allows the formation of long-range ordering. In practice, the source of such anisotropies can stem, e.g., from spin-orbit coupling or lattice distortions, leading to persistent magnetism with finite critical temperatures [68]. By mapping the relevant intersite interactions onto effective lattice models, real 2D magnetic vdW materials thus serve as an ideal platform to theoretically investigate emerging physical behaviors with their respective spin Hamiltonians. Depending on the spin degree of freedom, the isotropic Heisenberg model (in the absence of magnetic anisotropy), the XY-model (easy-plane anisotropy) or the Ising model (easy-axis anisotropy) are appropriate to study the underlying physics.

In addition to empirically fitting the effective model parameters to experimental measurements, electronic structure methods have been employed to extract them from first principles calculation, in particular by computing the exchange coupling J_ij that describe the interactions between the spins on site i and j of the lattice (see 3.1.2). Depending on the specific method employed and whether bulk or monolayer CrI₃ is considered, the exact values of J_ij may somewhat differ [59, 75–77]. In Solovyev et al [77], the intralayer parameters are denoted with ${J}_{1}^{\,\rm{Solovyev}\,}$, ${J}_{3}^{\,\rm{Solovyev}\,}$, and ${J}_{6}^{\,\rm{Solovyev}\,}$, whereas the interlayer interactions are denoted with ${J}_{2}^{\,\rm{Solovyev}\,}$, ${J}_{4}^{\,\rm{Solovyev}\,}$, and ${J}_{5}^{\,\rm{Solovyev}\,}$. In this work, we approximate a single layer CrI₃ with a two-dimensional hexagonal lattice by neglecting all interlayer exchange interactions, and by only taking into account the NN and NNN intralayer parameters ${J}_{1}^{\,\rm{Solovyev}\,}$ and ${J}_{3}^{\,\rm{Solovyev}\,}$, respectively, which in our notation are renamed to ${J}_{1}={J}_{\langle ij\rangle }={J}_{1}^{\,\rm{Solovyev}\,}$ and ${J}_{2}={J}_{\langle \langle ij\rangle \rangle }={J}_{3}^{\,\rm{Solovyev}\,}$ at the beginning of this section.

The current state-of-the-art classical method for simulating spin systems is based on Matrix Product States (MPS) paired with the Density Matrix Renormalization Group (DMRG) method [78]. The amount of entanglement an MPS can describe is bounded by the logarithm of its bond dimension χ, and therefore it is only efficient in one-dimensional systems, since their ground states follow an area law which can be efficiently represented by an MPS in one dimension [42, 79]. In higher dimensions, the area law conjecture states that constant-gapped and locally interacting lattice quantum systems also satisfy the area law, thus leading to the requirement that the bond dimension of the MPS scales exponentially with system size⁷ [2, 80]. While MPS are one-dimensional representatives of the family of tensor networks, their higher-dimensional extensions, such as projected entangled pair states [81, 82] or states generated by entanglement renormalization [83], can describe two-dimensional systems that satisfy an area law. Unfortunately, unlike for MPS in one dimension, computing expectation values with Projected Entangled Pair States (PEPS) in two dimensions is in general inefficient and approximate methods for the tensor contractions have to be employed [84, 85]. Therefore, spin systems in two or higher dimensions are natural candidates for the application of ground state preparation quantum algorithms such as QITE.

We consider the Transverse-Field Ising Model (TFIM), which is a paradigmatic model for the appearance of quantum phase transition at a critical point at zero temperature for short- and long-range interactions [86–88] and has recently been studied at scale with a neutral atom quantum simulator [89]. A similar model on a triangular lattice has been studied using QITE in [90]. Its Hamiltonian is given by

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{i\lt j}{J}_{ij}{\hat{\sigma }}_{i}^{x}{\hat{\sigma }}_{j}^{x}+B\displaystyle \sum _{i}{\hat{\sigma }}_{i}^{z},\end{eqnarray} \tag{ 31 }$$

where B sets the strength of the transverse magnetic field. The TFIM is an extension of the classical Ising model (that does not include a transversal field) and, unlike its classical version, displays a quantum phase transition at criticality even in one dimension. For weak long-range interactions in equation (30) with α ≥ 1, the one-dimensional TFIM can be described well by a classical GHF method, while in the long-range regime α < 1 the GHF solutions start deviating from DMRG results as α gets closer to the uniform interaction limit [91]. Thus, the TFIM is a model where we can apply QITE to a problem where we can control the quality of the initial state ∣Ψ_init〉 simply by tuning the parameter α.

The Hubbard model was introduced to describe the Mott transition from a conductor to an insulator and is one of the simplest many-particle theories that cannot be reduced to a single particle theory without resorting to some sort of approximation. We consider a minimal one-dimensional Hubbard model with L sites for a band structure in the atomic limit at half filling, n_el = L, whose Hamiltonian is given by [62]

$$\begin{eqnarray}&&\hat{H}=-t\displaystyle \sum _{\langle p,q\rangle }{\hat{E}}_{pq}+\frac{U}{2}\displaystyle \sum _{p=1}^{L}\left({\hat{E}}_{pp}^{2}-{\hat{E}}_{pp}\right)\end{eqnarray} \tag{ 32 }$$

with periodic boundary conditions, where t describes the hopping and U the onsite interaction parameter, and we introduced the fermionic orbital excitation operators as defined in equation (25). Here, the indices p, q describe the sites of the lattice containing L fermionic sites and 〈p, q〉 describes NN sites in the lattice. As we can see from the definition of the fermionic excitation operators ${\hat{E}}_{pq}$ in equation (25), the spin is implicitly contained in the operators. While the Hamiltonian in equation (32) is exactly solvable in one dimension [92], it serves as first test for fermionic-QITE applied to a local fermionic lattice Hamiltonian, before moving on to the more complicated molecular electronic structure Hamiltonians in section 3.2.

Section 4.1 describes a special cases of Generalized Hartree–Fock theory (GHF) [126], which is the subject of this section. GHF describes a procedure for finding the Fermionic Gaussian State (FGS) that possesses the lowest energy expectation value of a given Hamiltonian. A pure FGS of L fermionic modes can be written as

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle ={\hat{U}}_{{\bf{Q}}}| \,\rm{vac}\,\rangle ,\end{eqnarray} \tag{ 37 }$$

where Q ∈ O(2L) is a matrix that fully characterizes the FGS, and ${\hat{U}}_{{\bf{Q}}}$ is the generator of the FGS which can be written as the exponential of a quadratic polynomial of Majorana operators whose coefficients depend on Q, as detailed in Appendix A.3. An FGS is fully described by its respective (2L × 2L) covariance matrix Γ which is related to equation (37) through

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mu \nu }=\frac{i}{2}\langle {{\rm{\Psi }}}_{\,\rm{FGS}\,}| [{\hat{A}}_{\mu },{\hat{A}}_{\nu }]| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle ,\end{eqnarray} \tag{ 38 }$$

where μ, ν ∈ [2L] and $\hat{{\bf{A}}}={({\hat{A}}_{1},\ldots ,\,{\hat{A}}_{2L})}^{T}={({\hat{a}}_{1},{\hat{a}}_{3},\ldots ,\,{\hat{a}}_{2L-1},{\hat{a}}_{2},{\hat{a}}_{4},\ldots ,\,{\hat{a}}_{2L})}^{T}$ is a vector of Majorana operators that are defined as ${\hat{a}}_{2p-1}={\hat{c}}_{p}^{\dagger }+{\hat{c}}_{p}$ and ${\hat{a}}_{2p}=i({\hat{c}}_{p}^{\dagger }-{\hat{c}}_{p})$ for p ∈ [L]. As shown in [18] and detailed in Appendix A.1-A.2, one can find the covariance matrix Γ which describes the state in equation (37) with minimal energy w.r.t. a given fermionic Hamiltonian at a classical computational cost of ${ \mathcal O }({L}^{3})$ (attributed to the evaluation of the Pfaffian of the associated covariance matrix). The iterative method for finding the minimum is based on ITE of the covariance matrix and requires the evaluation of the energy expectation value of the fermionic Hamiltonian w.r.t. a given FGS, as well as its gradient, which both can be computed efficiently on a classical computer through its covariance matrix using Wick’s theorem [127]. Given an optimized covariance matrix Γ describing a pure FGS, we derive the explicit form of ${\hat{U}}_{{\bf{Q}}}$ in equation (37) for both possible cases, Γ describing an odd (even) parity FGS, i.e. Pf(Γ) = + 1 (−1). This then enables one to construct the initial FGS ∣Ψ_init〉 = ∣Ψ_FGS〉 on a quantum computer using ${ \mathcal O }({L}^{2})$ gates and a circuit depth that scales as ${ \mathcal O }(L)$ [16].

In order to apply the GHF to a general spin Hamiltonian, one has to first map the spin Hamiltonian to a fermionic Hamiltonian through e.g. the Jordan-Wigner transformation [128]. This can lead to fermionic Hamiltonians that are polynomials of Majorana operators of degree 2L. In addition, the resulting fermionic Hamiltonian does necessarily have to conserve the fermionic particle number, which is the case for the TFIM studied in this work. In case the resulting fermionic Hamiltonian does in fact conserve the particle number, one can also restrict equation (37) to covariance matrices that conserve the fermionic particle number. These number-conserving FGS are just the Slater determinants we introduced in equation (34) and GHF reduces to HF theory. We note, that GHF can be applied to lattice systems in arbitrary dimension.

a. HM I The first system we consider is identical to a system studied in [6] and is used in order to see if the QITE implementation reproduces similar results. As shown in table 2, HM I describes a one-dimensional spin chain with periodic boundary conditions, see figure 5(a), whose Hamiltonian is given by equation (29) with NN interactions J_〈ij〉 = 1 and no magnetic field, B = 0. Figure 6 shows the energies and fidelities of the spin-QITE and trotterized ITE. At the smallest Manhattan distance ν = 1 (i.e. ${\rm{\max }}(| {D}_{l}| )=4$), spin-QITE starts deviating from trotterized ITE after only a handful of steps, displaying an oscillatory behavior in both energy and fidelity far from the exact ground state energy. This behavior seems to be independent of the step size, as can be seen by the similar evolution for Δτ = 0.01 and Δτ = 0.001. We find that already at ν = 2 (${\rm{\max }}(| {D}_{l}| )=6$), spin-QITE is able to reproduce the trotterized ITE evolution and the energy saturates at a value E_QITE > E_ITE which is above the converged trotterized energy value E_ITE for step size Δτ = 0.1. The discrepancy in energy between QITE and trotterized ITE is due to finite value ν, and they should coincide at ν* = 4, where ${\rm{\max }}(| {D}_{l}| )=L$. One can also see that the trotterized ITE lies slightly above the exact ground state energy, which is due to the trotterization error discussed in equation (5), and can be removed by either using a smaller step size for the evolution, or using a higher order Trotter formula. The results agree well with the results reported in [6].

Zoom In Zoom Out Reset image size

b. HM II In figure 7 we show the results obtained for the one-dimensional long-range Hamiltonian that describes system HM II, where the exponent in equation (30) that sets the strength of the interaction is α = 1. Note, that since the Hamiltonian is long-range, the resulting domain sizes are larger than for NN Hamiltonians for the same value ν. Our numerical simulations were limited to ν = 1 not just for this model, but also for all lattice models that either included beyond-NN interactions or a lattice dimension larger than one, due to the large computational demands that the scaling ${4}^{| {D}_{l}| }$ leads to. We probe the QITE evolution for discretization step sizes Δτ = 0.1 and Δτ = 0.01, respectively. We find that for both realizations, QITE almost immediately deviates strongly from the trotterized ITE, indicating that, just as with the HM I model, ν = 1 is not sufficient. A QITE simulation with ν = 2 would here lead to ${\rm{\max }}(| {D}_{l}| )=10$, thus a 256-fold increase in computational cost (both on the classical side, and regarding the number of operators that need to be measured by the quantum computer). The long-term behavior shows an oscillatory behavior with a fidelity that quickly drops off to zero, see figure 19.

Zoom In Zoom Out Reset image size

c. HM III In figure 8 we show the results obtained for the Heisenberg model on a two-dimensional triangular lattice with periodic boundary conditions in one direction and two rows of six spins each, leading to a total of L = 12, as shown in figure 5(b). Triangular lattices are of particular interest when studying antiferromagnetic systems, since they can lead to frustration [89]. Like HM I, the system HM III only considers antiferromagnetic NN interactions. Unlike HM I, QITE is already stable for ν = 1 which corresponds to a domain size of ${\rm{\max }}(| {D}_{l}| )=7$ ⁸ . The discrepancy between the converged QITE energy E_QITE and the converged trotterized ITE energy E_QITE can be explained by the small domain size.

Zoom In Zoom Out Reset image size

a. TFIM I In figure 9 we present results of the simulation of the TFIM I system, which describes the transverse-field Ising model with long-range interactions in the strong-long-range regime α = 0.3 [87]. For the one-dimensional chain of a system displaying long-range interactions, ${\rm{\max }}(| {D}_{l}| )=6$ for ν = 1. Similar to the HM II model which also describes long-range interactions, one can observe that for ν = 1 the QITE quickly deviates from the monotonic behavior of trotterized ITE, independent of the time step size being Δτ = 0.1 or Δτ = 0.01. Thus, as observed in the long-range HM II, a larger ν is required in order to be able to improve the initial state overlap. The long-time behavior up to τ = 10 is presented for completeness in figure 20.

Zoom In Zoom Out Reset image size

b. TFIM II System TFIM II using the same lattice geometry as TFIM I considers the even stronger long-range regime α = 0.1, reflected in the lower initial state overlap. The step sizes we consider here are Δτ = 0.01 for the trotterized ITE and Δτ = ∈ {0.01, 0.001} for QITE. The smaller step sizes are chosen in order to be able to resolve if an improvement over the initial state overlap can be achieved in the first part of the evolution. The results display in figure 10 show that QITE does indeed improve until τ ≈ 0.1, before quickly dropping off to zero similar to figure 20.

Zoom In Zoom Out Reset image size

The last spin system we consider is the two-dimensional J₁J₂-model on a hexagonal lattice as depicted in figure 5(c) with NN and NNN interactions. Similar to the previous NN models shown in figures 6 and 8, the results for this system in figure 11 display an initially stable behavior of QITE regarding the energy, until τ ≈ 0.5, even for the smallest non-trivial Manhattan distance ν = 1, which corresponds to ${\rm{\max }}(| {D}_{l}| )=7$. Interestingly, the fidelity continues to improve well beyond the point where the energy starts growing again, and peaks at around τ ≈ 2.5, before worsening again. We observe that the fidelity can be increased from around 40% to 90% even though we are only able to simulate the smallest Manhattan distance ν = 1. We also note that this presents heuristic evidence, that the behavior of the energy might not be a reliable tool for determining the performance of QITE, whose main purpose is to amplify the overlap γ, and not necessarily the improvement of the energy.

Zoom In Zoom Out Reset image size

Note, that this is only considering a simplified model of the CrI₃ description introduced in section 3.1.1 with periodic boundary in one direction. In order to be able to simulate the larger lattice 5(d) of CrI₃ with periodic boundary conditions in both directions and include a study of Manhattan distance of ν = 2, the code base would have to be improved significantly since this would lead to domain sizes of ${\rm{\max }}(| {D}_{l}| )=14$.

As a proof of principle, we study the one-dimensional FHM system described in table 2, whose Hamiltonian is given by equation (32) and allows for hopping between NN sites. The ITE and QITE results are displayed in figure 12. Here, ν represents the number of additional sites to be included in the domain of a Hamiltonian term (in addition to the Hamiltonian term’s support). Each site contains two fermionic modes, so the system we consider with L = 10 corresponds to 20 fermionic modes. For the FHM, ν describes the number of additional sites added to the domain, which are picked according to their Manhattan distance from the support sites, as explained in section 2.3.2.

Zoom In Zoom Out Reset image size

For ν = 1, we observe an early deviation of QITE from ITE similar to the HM I system in figure 6. This can be an effect that is due to the fact that ν = 1 considers only one of the two neighboring sites (chosen randomly) of the Hamiltonian term support in figure 2(b). For ν = 2, both the left- and right-hand side sites are included, which corresponds to a Manhattan distance radius of size one. The behavior of ν = 2 already leads to an almost exact replication of the trotterized ITE. Interestingly, as shown in figure 21, we find that the long-time behavior is less stable for larger domain sizes ν = 4 than for ν = 2 and ν = 3.

The first electronic structure system we study is the Neon atom described in section 3.2. The AS is characterized by (n_el, n_orb) = (8, 8), and the mutual information of its ground state is shown in figure 30. Figure 13 shows the results of fermionic-QITE for a range of domain sizes which are characterized by the number of additional orbitals ν included to the orbitals of the support. One can see that at step size Δτ = 0.01 all ν initially lead to an improvement of the fidelity and energy, but only the exact fermionic-QITE simulation ν^* = 7, where the domain size is always identical to the number of AS orbitals ∣D_l∣ = L, is stable. One also observes that the discretization step size Δτ = 0.01 leads to an almost exact replication of the trotterized ITE behavior for ν^*, compared with a more coarse grained discretization step size Δτ = 0.1, where the energy is decreasing at a lower rate. Since at ν^* there is no trotterization error anymore, the only difference between QITE and ITE is the error due to the truncated operator expansion and the error from the taylor expansion of the imaginary time propagator. The long-time behavior of the evolution is shown in figure 22. The low diagnostic value Z_s(1) = 0.0441 in table 3 indicates that a single Slater determinant could provide an initial state with a large fidelity, which is confirmed by the large fidelity F_ROHF ≈ 97.9%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now consider the fermionic-QITE results for the (5, 5) AS of the Fe(III)-NTA (0) system described in table 1 and section 3.2. The mutual information used for picking the orbitals of the domain are given in figure 31. The initial state has a fidelity of F_ROHF = 89% and the fidelity only gradually increases. Unlike the results for Ne (see section 5.2.1), all approximate QITE simulations with ν < ν^* are at the beginning following the trotterized ITE evolution and start to deviate from it around τ = 0.5. Interestingly, a finer step size Δτ = 0.01 leads to an earlier deviation of the approximate QITE from trotterized ITE in comparison to Δτ = 0.1. The only stable behavior can be observed for exact QITE at ν^* = 4.

In figure 15 we show the results for Fe(III)-NTA with three unpaired electrons. The mutual information is depicted in figure 32 and shows significantly less overall correlation than in the case of only one unpaired electron of figure 31. One can see that approximate QITE of system Fe(III)-NTA (3) is much more stable than Fe(III)-NTA (1), and only starts to deviate at a time τ ≈ 5 (compared to τ ≈ 0.5 for Fe(III)-NTA (1)), see figure 24. In addition, Z_s(1) = 0.049 and F_ROHF ≈ 97.6% for three unpaired electrons, while for one unpaired electron Z_s(1) = 0.1839 and F_ROHF ≈ 89%. Therefore, the high spin state has significantly less static correlation than the low spin state.

Zoom In Zoom Out Reset image size

We now consider the three molecular forms of oxygen introduced in section 3.2. Figure 16 shows the short-term behavior of singlet oxygen O₂. Similar to the strongly correlated Fe(III)-NTA complex with one unpaired electron, trotterized ITE produces a state whose fidelity increases only at a very slow rate. The only stable QITE realization is exact QITE which includes up to ν^* = 4 orbitals to the support. In order to check if the ITE actually reaches the ground state, we simulate the long-time behavior of QITE and trotterized ITE in figure 25 until an evolution time τ = 10 and perform a simulation of the exact (un-trotterized) ITE up to τ = 100 in figure 26 to confirm that the ground state of O₂(0) is reached in the long-time limit of the ITE. The fidelity of the ROHF solution w.r.t. the ground state F_ROHF ≈ 46% is by far the lowest across all studied molecular instances, see table 3, and the multi-reference diagnostic also gives the highest value Z_s(1) = 0.2607 across all chemical systems, indicating the strong multi-reference character of singlet oxygen. It should be noted, that for this system, we also find that the first excited state of the AS Hamiltonian possesses a fidelity of approximately F = 45% w.r.t. the ROHF solution. Figure 17 shows the results for triplet oxygen, O₂ (2). Even the domain size at ν = 1 is able to improve the fidelity of the evolved state from initially F_ROHF ≈ 93% to F ≈ 98%. The overall mutual information of the ground states of the singlet oxygen (see figure 33) also shows that it is significantly larger than that of triplet oxygen (see figure 34), similar to the behavior of the Fe(III)-NTA systems. We want to stress that we are not implying that the two systems, Fe(III)-NTA and oxygen are in any way similar, only that the qualitative behavior of their (Q)ITE is similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The last system we consider is singlet molecular ozone in a (12, 9) AS in figure 18. Its respective long time behavior and mutual information are shown in figures 28 and 4, respectively. The multi-reference diagnostic and ROHF fidelity are Z_s(1) = 0.186 and F_ROHF  ≈ 81%, indicative of a strong multi-reference character of the ground state. This is similar to what has been observed for the smaller AS of the Fe(III)-NTA (1) system in section 5.2.2. In the O₃ system, the smaller step size Δτ = 0.01 simulations of QITE were systematically better performing than the more coarse grained size Δτ = 0.1 simulations, and even the smallest domains at ν = 1 initially improves the fidelity until τ ≈ 0.5. The fact that exact QITE, here at ν^* = 8, still displays a gap to the exact ITE is due to the first-order approximation in terms of the step size Δτ when setting up the linear system of equations in equation (22). This can be seen more clearly in figure 29, where we zoom in on the long-time behavior of ITE and exact QITE.

Zoom In Zoom Out Reset image size

The even case has been described e.g. in [31], and is included for completeness. For even parity sectors, we have

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle ={\hat{U}}_{{\bf{R}}}| \,\rm{vac}\,\rangle ={\hat{U}}_{{\bf{R}}}| {0}_{1},\ldots ,{0}_{L}\rangle \end{eqnarray} \tag{ A14 }$$

where Q = R ∈ SO(2L), i.e. $\det ({\bf{R}})=1$ and

$$\begin{eqnarray}&&{\hat{U}}_{{\bf{R}}}={e}^{\frac{1}{4}{\hat{{\bf{A}}}}^{T}\mathrm{log}({\bf{R}})\hat{{\bf{A}}}}.\end{eqnarray} \tag{ A15 }$$

Here, $\mathrm{log}({\bf{R}})$ denotes the matrix logarithm of R and we introduced the p-q ordering of the covariance matrix (meaning it is ordered as $\hat{{\bf{A}}}={({\hat{a}}_{1},{\hat{a}}_{3},\ldots ,\,{\hat{a}}_{2L-1},{\hat{a}}_{2},{\hat{a}}_{4},\ldots ,\,{\hat{a}}_{2L})}^{T}$, which we will indicate by a tilde over the covariance matrix, $\tilde{{\rm{\Gamma }}}$). Using the relation

$$\begin{eqnarray}&&{\hat{U}}_{{\bf{R}}}^{\dagger }{\hat{A}}_{\mu }{\hat{U}}_{{\bf{R}}}=\displaystyle \sum _{\nu =1}^{2L}{R}_{\mu \nu }{\hat{A}}_{\nu },\end{eqnarray} \tag{ A16 }$$

we can compute the following relation

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\rm{\Gamma }}}}_{\mu \nu } & = & \frac{i}{2}\langle {{\rm{\Psi }}}_{\,\rm{FGS}\,}| [{\hat{A}}_{\mu },{\hat{A}}_{\nu }]| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle \\ & = & \frac{i}{2}\langle \,\rm{vac}\,| {\hat{U}}_{{\bf{Q}}}^{\dagger }[{\hat{A}}_{\mu },{\hat{A}}_{\nu }]{\hat{U}}_{{\bf{Q}}}| \,\rm{vac}\,\rangle \\ & = & i\displaystyle \sum _{\lambda ,\sigma =1}^{2L}{Q}_{\mu \sigma }{Q}_{\nu \lambda }\langle \,\rm{vac}\,| {\hat{A}}_{\sigma }{\hat{A}}_{\lambda }| \,\rm{vac}\,\rangle .\end{array}\end{eqnarray} \tag{ A17 }$$

We will now divide the possible values for σ, λ ∈ [2L] into four sectors and analyze them separately. Note, that ${\tilde{{\rm{\Gamma }}}}_{\sigma \sigma }=0$ and we therefore assume σ ≠ λ in the following.

• 1.  
  Sector σ, λ ∈ [1, …, L]: Here, due to the p-q-ordering, we have both ${\hat{A}}_{\sigma }$ and ${\hat{A}}_{\lambda }$ of the form $({\hat{c}}_{p}^{\dagger }+{\hat{c}}_{p})$, where p ∈ [L]. We therefore have to evaluate terms such as

  $$\begin{eqnarray}\begin{array}{rcl}\langle \,\rm{vac}\,| ({\hat{c}}_{p}^{\dagger }+{\hat{c}}_{p})({\hat{c}}_{q}^{\dagger }+{\hat{c}}_{q})| \,\rm{vac}\,\rangle & = & \langle \,\rm{vac}\,| {\hat{c}}_{p}^{\dagger }{\hat{c}}_{q}| \,\rm{vac}\,\rangle \\ & & +\langle \,\rm{vac}\,| {\hat{c}}_{p}{\hat{c}}_{q}^{\dagger }| \,\rm{vac}\,\rangle \\ & = & 0.\end{array}\end{eqnarray} \tag{ A18 }$$

  since we assumed that σ ≠ λ which here translated to p ≠ q. Thus, $\langle \,\rm{vac}\,| {\hat{A}}_{\sigma }{\hat{A}}_{\lambda }| \,\rm{vac}\,\rangle =0$ for σ, λ ∈ [L].
• 2.  
  Sector σ, λ ∈ [L + 1, …, 2L]:Similar arguments apply and we see that $\langle \,\rm{vac}\,| {\hat{A}}_{\sigma }{\hat{A}}_{\lambda }| \,\rm{vac}\,\rangle =0$ for σ, λ ∈ [L + 1, …, 2L].
• 3.  
  Sector σ ∈ [1, …, L] and λ ∈ [L + 1, …, 2L]:Here we get terms of the form

  $$\begin{eqnarray}&&\langle \,\rm{vac}\,| ({\hat{c}}_{p}^{\dagger }+{\hat{c}}_{p})i({\hat{c}}_{q}^{\dagger }-{\hat{c}}_{q})| \,\rm{vac}\,\rangle =i{\delta }_{pq}.\end{eqnarray} \tag{ A19 }$$
• 4.  
  Sector σ ∈ [L + 1, …, 2L] and λ ∈ [1, …, L]:Here we get terms of the form

  $$\begin{eqnarray}&&\langle \,\rm{vac}\,| i({\hat{c}}_{p}^{\dagger }-{\hat{c}}_{p})({\hat{c}}_{q}^{\dagger }+{\hat{c}}_{q})| \,\rm{vac}\,\rangle =-i{\delta }_{pq}.\end{eqnarray} \tag{ A20 }$$

By defining the even-parity vacuum covariance matrix in the p-q-ordering,

$$\begin{eqnarray}{\tilde{{\boldsymbol{\Upsilon }}}}^{(+1)}=\left(\begin{array}{cc}{{\bf{0}}}_{L} & {{\bf{1}}}_{L}\\ -{{\bf{1}}}_{L} & {{\bf{0}}}_{L}\end{array}\right),\end{eqnarray} \tag{ A21 }$$

we can write equation (A17) as

$$\begin{eqnarray}&&\tilde{{\boldsymbol{\Gamma }}}=-{\bf{Q}}{\tilde{{\boldsymbol{\Upsilon }}}}^{(+1)}{{\bf{Q}}}^{T}\end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray}&&=-{e}^{i{\boldsymbol{\xi }}}{\tilde{{\boldsymbol{\Upsilon }}}}^{(+1)}{e}^{-i{\boldsymbol{\xi }}},\end{eqnarray} \tag{ A23 }$$

where ${\boldsymbol{\xi }}=-i\mathrm{log}({\bf{R}})$ is a purely imaginary and skew-symmetric matrix.

Let us now consider the case of odd parity. Here, Q ∈ O(2L) with $\det ({\bf{Q}})=-1$. In this case we have

$$\begin{eqnarray}\begin{array}{rcl}| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle & = & {\hat{U}}_{{\bf{Q}}}| \,\rm{vac}\,\rangle \\ & = & {\hat{U}}_{\bar{{\bf{R}}}}{X}_{L}| \,\rm{vac}\,\rangle \\ & = & {\hat{U}}_{\bar{{\bf{R}}}}| {1}_{L}\rangle ,\end{array}\end{eqnarray} \tag{ A24 }$$

where we use the simplified notation ∣1_L〉 = ∣0₁, 0₂, …, 0_L−1, 1_L〉. Here, $\bar{{\bf{R}}}\in \,\rm{SO}\,(2L)$, where

$$\begin{eqnarray}&&{\bar{R}}_{\mu \nu }=\left\{\begin{array}{l}{Q}_{\mu \nu },\,\,\rm{if}\,\,\nu \ne 2L\\ -{Q}_{\mu \nu },\,\,\rm{if}\,\,\nu =2L.\end{array}\right.\end{eqnarray} \tag{ A25 }$$

We know from equation (A16) that for special orthogonal matrices we have

$$\begin{eqnarray}&&{\hat{U}}_{\bar{{\bf{R}}}}^{\dagger }{\hat{A}}_{\mu }{\hat{U}}_{\bar{{\bf{R}}}}=\displaystyle \sum _{\nu =1}^{2L}{\bar{R}}_{\mu \nu }{\hat{A}}_{\nu }.\end{eqnarray} \tag{ A26 }$$

Similar to the even parity sector calculation, we now compute

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\rm{\Gamma }}}}_{\mu \nu } & = & \frac{i}{2}\langle {{\rm{\Psi }}}_{\,\rm{FGS}\,}| [{\hat{A}}_{\mu },{\hat{A}}_{\nu }]| {{\rm{\Psi }}}_{\,\rm{FGS}\,}\rangle \\ & = & \frac{i}{2}\langle {1}_{L}| {\hat{U}}_{\bar{{\bf{R}}}}^{\dagger }[{\hat{A}}_{\mu },{\hat{A}}_{\nu }]{\hat{U}}_{\bar{{\bf{R}}}}| {1}_{L}\rangle \\ & = & i\displaystyle \sum _{\lambda \ne \sigma }^{2L}{\bar{R}}_{\mu \sigma }{\bar{R}}_{\nu \lambda }\langle {1}_{L}| {\hat{A}}_{\sigma }{\hat{A}}_{\lambda }| {1}_{L}\rangle .\end{array}\end{eqnarray} \tag{ A27 }$$

We again consider four different sectors for the values of σ, λ ∈ [2L]. Following the same approach as before, we find that by defining the (L × L)-diagonal matrix

$$\begin{eqnarray}{\bar{{\bf{1}}}}_{L}=\left(\begin{array}{ccccc}1 & 0 & \cdots \, & 0 & 0\\ 0 & 1 & \cdots \, & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots \, & 1 & 0\\ 0 & 0 & \cdots \, & 0 & -1\end{array}\right)\end{eqnarray} \tag{ A28 }$$

which is identical to the identity matrix except for the flipped sign in the last row, and defining the odd-parity vacuum covariance matrix in the p-q-ordering

$$\begin{eqnarray}{\tilde{{\boldsymbol{\Upsilon }}}}^{(-1)}=\left(\begin{array}{cc}{{\bf{0}}}_{L} & {\bar{{\bf{1}}}}_{L}\\ -{\bar{{\bf{1}}}}_{L} & {{\bf{0}}}_{L}\end{array}\right)\end{eqnarray} \tag{ A29 }$$

we can write

$$\begin{eqnarray}\begin{array}{rcl}\tilde{{\boldsymbol{\Gamma }}} & = & -\bar{{\boldsymbol{R}}}{\tilde{{\boldsymbol{\Upsilon }}}}^{(-1)}{\bar{{\boldsymbol{R}}}}^{T}\\ & = & -{e}^{i\bar{{\boldsymbol{\xi }}}}{\tilde{{\boldsymbol{\Upsilon }}}}^{(-1)}{e}^{-i\bar{{\boldsymbol{\xi }}}},\end{array}\end{eqnarray} \tag{ A30 }$$

where $\bar{{\boldsymbol{\xi }}}=-i\mathrm{log}(\bar{{\bf{R}}})$ is a purely imaginary and skew-symmetric matrix.

By similar considerations as in [91], the energy expectation value of the Heisenberg Hamiltonian in equation (29) of a fermionic Gaussian state described by the covariance matrix Γ is given by

$$\begin{eqnarray}\begin{array}{rcl}E({\boldsymbol{\Gamma }}) & = & \displaystyle \sum _{i\lt j}{J}_{ij}\left[{(-1)}^{j-i}\left\{\,\rm{Pf}\,\left({\left.{\boldsymbol{\Gamma }}\right|}_{2i,2i+1,\ldots ,\,2j-2,2j-1}\right)\right.\right.\\ & & \left.-\,\rm{Pf}\,\left({\left.{\boldsymbol{\Gamma }}\right|}_{2i-1,2i+1,2i+2,\ldots ,2j-3,2j-2,2j}\right)\right\}\\ & & \left.+\,\rm{Pf}\,\left({\left.{\boldsymbol{\Gamma }}\right|}_{2i-1,2i,2j-1,2j}\right)\right]\\ & & -\frac{B}{2}\displaystyle \sum _{i}\left({{\rm{\Gamma }}}_{2i-1,2i}-{{\rm{\Gamma }}}_{2i,2i-1}\right).\end{array}\end{eqnarray} \tag{ A40 }$$

We can compute the mean-field matrix ${\bf{F}}=4\frac{\partial E({\boldsymbol{\Gamma }})}{\partial {\boldsymbol{\Gamma }}}$ of equation (A40) using a Pfaffian gradient identity from equation (A11).

The energy expectation value of the transverse field Ising model Hamiltonian in equation (31) of a fermionic Gaussian state is given by [91]

$$\begin{eqnarray}\begin{array}{rcl}E({\boldsymbol{\Gamma }}) & = & \displaystyle \sum _{i\lt j}{J}_{ij}{(-1)}^{j-i}\,\rm{Pf}\,\left({\left.{\boldsymbol{\Gamma }}\right|}_{2i,2i+1,\ldots ,\,2j-2,2j-1}\right)\\ & & -\frac{B}{2}\displaystyle \sum _{i}\left({{\rm{\Gamma }}}_{2i-1,2i}-{{\rm{\Gamma }}}_{2i,2i-1}\right).\end{array}\end{eqnarray} \tag{ A41 }$$

The mean-field matrix F can be derived from equation (A11).