Ab initio extended Hubbard model of short polyenes for efficient quantum computing

We briefly review the extended Hubbard Hamiltonian of the ab initio downfolding method Aryasetiawan et al. (2004); Nakamura et al. (2021). For explanation, the effective Hamiltonian for molecules is explicitly described, with minor modifications not including the summations of inter-unit cell interactions.

We employ an extended Hubbard Hamiltonian to describe isolated molecules, which is given by,

$\displaystyle\mathcal{H}_{\textrm{ex-Hub}}=\mathcal{H}_{0}+\mathcal{H}_{\rm int}.$

(1)

Here, $\mathcal{H}_{0}$ and $\mathcal{H}_{\rm int}$ consists of one-body and two-body operators, respectively. The one-body term is expressed as:

$\displaystyle\mathcal{H}_{0}=\sum_{i,j,\sigma}(t_{ij}-t_{ij}^{\rm DC})c^{\dagger}_{i\sigma}c_{j\sigma},$

(2)

where $c^{(\dagger)}_{i\sigma}$ denote the annihilation (creation) operators of the $i$-th orbital with spin $\sigma$. Here, $t_{ij}$ is a matrix element of the effective one-electron operator defined by

$\displaystyle t_{ij}=\int\textrm{d}{\bm{r}}\phi_{i}^{*}({\bm{r}})\mathcal{H}_{\textrm{KS}}\phi_{j}({\bm{r}}),$

(3)

where $\phi_{i}({\bm{r}})$ is the $i$-th Wannier orbital. This integral is executed over the crystal volume, and $\mathcal{H}_{\rm KS}$ is the Kohn-Sham (KS) Hamiltonian.

The term $t_{ij}^{\textrm{DC}}$ is introduced to prevent double counting (DC) of two-electron integrals from the one-electron integral $t_{ij}$. It can be defined as

$\displaystyle t_{ij}^{\rm DC}$

$\displaystyle\equiv\begin{cases}\alpha U_{ii}D_{ii}+\sum_{k\neq i}U_{ik}D_{kk}&(i=j)\\ 0&(i\neq j),\end{cases}$

(4)

following the approach in Ref. Kawamura et al. (2017). $\alpha$ is a parameter ranging $0\leq\alpha\leq 1$ for tuning the correction. The DC of the electron-electron interactions arises from the exchange-correlation functional of KS–density functional theory (DFT), and it should be noted that it cannot be eliminated completely. The one-body density matrix $D_{ij}$ can be defined in the KS orbital basis as:

$\displaystyle D_{ij}$

$\displaystyle\equiv\sum_{\sigma}\left\langle c_{i\sigma}^{\dagger}c_{j\sigma}\right\rangle_{\rm KS}.$

(5)

$U_{ij}$ denotes a screened Coulomb integral represented as:

$\displaystyle U_{ij}(\omega)$

$\displaystyle=\int\textrm{d}{\bm{r}}\int\textrm{d}{\bm{r}}^{\prime}\phi_{i}^{*}({\bm{r}})\phi_{i}({\bm{r}})W({\bm{r}},{\bm{r}}^{\prime},\omega)\phi_{j}^{*}({\bm{r}}^{\prime})\phi_{j}({\bm{r}}^{\prime}).$

(6)

In this formulation, the inverse operator $1/|{\bm{r}}-{\bm{r}}^{\prime}|$ of the bare Coulomb interaction is replaced with the frequency-dependent screened Coulomb interaction $W({\bm{r}},{\bm{r}}^{\prime},\omega)$, where $\omega$ is frequency. The $U_{ij}(\omega)$ value can be derived from cRPA Aryasetiawan et al. (2004); Nakamura et al. (2021). For our purposes, we utilize the static limit of the frequency-dependent Coulomb integrals:

$\displaystyle U_{ij}$

$\displaystyle=\lim_{\omega\rightarrow 0}U_{ij}(\omega).$

(7)

For the two-body electron-electron interaction component, $\mathcal{H}_{\rm int}$, we compare the following terms: $\mathcal{H}_{\rm int}^{(1)}$ and $\mathcal{H}_{\rm int}^{(2)}$. These terms are expressed as follows:

	$\displaystyle{\cal H}_{\rm int}^{(1)}$	$\displaystyle=\frac{1}{2}\sum_{i,j}\sum_{\sigma,\rho}\left\{U_{ij}c_{i\sigma}^{\dagger}c_{j\rho}^{\dagger}c_{j\rho}c_{i\sigma}\right.$
		$\displaystyle\quad+\left.J_{ij}\left(c_{i\sigma}^{\dagger}c_{j\rho}^{\dagger}c_{i\rho}c_{j\sigma}+c_{i\sigma}^{\dagger}c_{i\rho}^{\dagger}c_{j\rho}c_{j\sigma}\right)\right\},$		(8)
	$\displaystyle{\cal H}_{\rm int}^{(2)}$	$\displaystyle=\frac{1}{2}\sum_{i,j}\sum_{\sigma,\rho}U_{ij}c_{i\sigma}^{\dagger}c_{j\rho}^{\dagger}c_{j\rho}c_{i\sigma}.$		(9)

Here, $J_{ij}$ is a screened exchange integral represented as:

$\displaystyle J_{ij}(\omega)$

$\displaystyle=\int\textrm{d}{\bm{r}}\int\textrm{d}{\bm{r}}^{\prime}\phi_{i}^{*}({\bm{r}})\phi_{j}({\bm{r}})W({\bm{r}},{\bm{r}}^{\prime},\omega)\phi_{j}^{*}({\bm{r}}^{\prime})\phi_{i}({\bm{r}}^{\prime}),$

(10)

where we also utilize the static limit for exchange integrals and refer to it as $J_{ij}$. The $\mathcal{H}_{\rm int}^{(2)}$ term discards the exchange interactions of the $\mathcal{H}_{\rm int}^{(1)}$ term.

Note that the ab initio downfolding method has scarcely been applied to isolated chemical systems; very recently, Chang et al. reported to apply to vanadocene Chang et al. (2023). Other related studies include the assessment of screened Coulomb interaction based on random phase approximation (RPA) for Nb_xCo ($1\leq x\leq 9$) clusters Peters et al. (2018) and benzene van Loon et al. (2021). Scott et al. recently developed the moment-constrained RPA to evaluate static effective interaction rather than relying on the static limit of cRPA Scott and Booth (2024).

Our approach utilizes the extended Hubbard model for an isolated molecule constructed via the ab initio downfolding method in quantum computation. The conceptual figure comparing our approach with the conventional one is shown in Figure 1.

In this approach, the simplified model Hamiltonian is expected to perform quantum computations with a relatively small and shallow quantum circuit. The two-electron operator of the extended Hubbard model Hamiltonian consists of the second-order tensors ${\bm{U}}$ and ${\bm{J}}$, whose indices belong to a smaller but physically essential space. It becomes a much sparse and compact representation compared to the ab initio Hamiltonian, the fourth-order tensor ${\bm{g}}$, whose indices belong to the entire space of orbitals.

It should be noted that our research direction is closely related to that of dynamical self-energy mapping (DSEM) Dhawan et al. (2021); Daniel et al. (2021). The DSEM procedure parametrizes a sparse Hamiltonian to reproduce the dynamical self-energy of the original molecular Hamiltonian. The sparse Hamiltonian contains at most $\mathcal{O}(N^{2})$ two-body interaction terms, and it makes quantum circuits shallower and increases the feasibility of quantum computation for molecular systems.

A second quantized Hamiltonian can be transformed by fermion-to-qubit mappings, such as Jordan-Wigner mapping, into the linear combination of the Pauli operators as

$\displaystyle\mathcal{H}=\sum_{l=1}^{N_{\rm term}}h_{l}P_{l},$

(11)

where $P_{l}$ is the $l$-th Pauli operator and $h_{l}$ is its coefficient. The summation runs over the number of the Pauli operators, denoted as $N_{\rm term}$.

$L^{1}$-norm of the coefficient vector of the fermion-to-qubit mapped Hamiltonian $\mathcal{H}$ is defined as the sum of the absolute values of the coefficients

$\displaystyle\lambda=\sum_{l=1}^{N_{\rm term}}|h_{l}|.$

(12)

Given the crucial role that the parameter $\lambda$ of a given Hamiltonian plays in determining the scalability of various quantum algorithms Lee et al. (2021); Koridon et al. (2021), evaluating the value of $\lambda$ is essential for demonstrating the feasibility of quantum computing. For example, the scaling of the phase estimation algorithms based on the qubitization method using the single factorization and the tensor hypercontraction of the Coulomb operator is $\tilde{\mathcal{O}}(N^{3/2}\lambda/\varepsilon)$ Berry et al. (2019) and $\tilde{\mathcal{O}}(N\lambda/\varepsilon)$ Lee et al. (2021), respectively. Here, $N$ is the number of spin orbitals, and $\varepsilon$ is the target precision. The qDRIFT algorithm Campbell (2019), a randomized compiler for Hamiltonian simulation, has the $\mathcal{O}(\lambda^{2}/\varepsilon^{2})$ scaling.

Measurement is one of the most troublesome processes on variational quantum eigensolver (VQE) Peruzzo et al. (2014); Tilly et al. (2022), and the number of measurements of VQE also depends on the $L^{1}$-norm $\lambda$ of the Hamiltonian Wecker et al. (2015); Rubin et al. (2018). The total number of measurements $M$ is given by the sum of each number of measurements $M_{l}$ for each Pauli operator $P_{l}$ of the Hamiltonian as

$\displaystyle M=\sum_{l=1}^{N_{\rm term}}M_{l}.$

(13)

The optimal $M_{l}$ is proportional to $|h_{l}|$, which is derived through the method of Lagrange multipliers Rubin et al. (2018). The optimal number of measurements is represented as

$\displaystyle M=\frac{1}{\varepsilon^{2}}\left(\sum_{l=1}^{N_{\rm term}}|h_{l}|\sigma_{l}\right)^{2}\leq\frac{\lambda^{2}}{\varepsilon^{2}},$

(14)

where $\sigma_{l}$ is the intrinsic standard deviation of the Pauli operator $P_{l}$. Here, $\varepsilon$ is the sampling error for the expectation value of the Hamiltonian. The right inequality of Eq. (14) is derived from the condition of the intrinsic variance $\sigma_{l}^{2}\leq 1$. Hence, the total number of measurements of VQE is bounded using the $L^{1}$-norm $\lambda$.

The Wannier functions $\{\phi_{i}\}$ of our models are shown in Figure 2.

The modeled active space is represented as ($m$e, $n$o), where $(m,n)$ is the number of electrons and orbitals. The shape of these Wannier functions is reasonable because they have the $\pi$-orbital character, which is a key to understanding the lower electron excitations of the polyenes. In the (4e, 4o) and (6e, 6o) cases of hexatriene, the shape of the Wannier functions is reasonable but different. The difference is that the former uses two fewer $\pi$-type KS orbitals for localization.

In particular, the Wannier functions of ethylene (2e, 2o), butadiene (4e, 4o), and hexatriene (6e, 6o) have shapes reminiscent of a $p_{z}$-orbital. This characteristic suggests that the model Hamiltonian with such a basis closely resembles the second-quantized Pariser-Parr-Pople (PPP) Hamiltonian Schulten and Karplus (1972); Tavan and Schulten (1986), which is also the extended Hubbard-type (analogous to the model with $\mathcal{H}_{\rm int}=\mathcal{H}_{\rm int}^{(2)}$) and composed of the set of the orthonormal $p_{z}$-orbital basis.

Note that the PPP model is a semi-empirical molecular orbital method for $\pi$-conjugated systems Pariser and Parr (1953a, b); Pople (1953). The PPP-multireference double excitation configuration interaction (PPP-MRD-CI) calculation, a configuration interaction (CI) calculation based on PPP, had achieved qualitative success for polyenes of relatively long conjugation length, even though the molecular integrals were given empirically. The PPP-MRD-CI calculations provided a qualitatively reasonable explanation for spectroscopic findings for the energetic order of the optically allowed single excitation and the forbidden double excitation Hudson and Kohler (1972, 1984), although the evaluated state order of butadiene is now known to be reversed Watson and Chan (2012).

The major difference between the second-quantized PPP and our extended Hubbard models is how the molecular integrals are determined. Whereas the former is given empirically, the latter is derived by the ab initio downfolding approach. The integral parameters of our models are summarized in Table 1.

These $U_{ij}$ and $J_{ij}$ parameters include the electron correlation effects from outside of the active space. In the next section, we examine the incorporation of the electron correlation effects by increasing the number of bands and extrapolating the continuum limit.

We investigate the excitation energies of the models varying the number of bands $N_{\rm band}$. The values of excitation energy $\Delta E$ are extrapolated to confirm the reliability of our models via estimating the continuum limit. For extrapolation, the following equation is used:

$\displaystyle f(N_{\rm band})=\Delta E_{\infty}+b\times{\rm exp}(-N_{\rm band}/c),$

(15)

where $\Delta E_{\infty}$, $b$, and $c$ are real fitting parameters. $\Delta E_{\infty}$ corresponds to the extrapolated value of the excitation energy. Fitting is performed based on the eigenvalues of the model Hamiltonian with $\mathcal{H}_{\rm int}=\mathcal{H}_{\rm int}^{(1)}$ obtained by H$\Phi$ Kawamura et al. (2017); Ido et al. (2023). In this section, we focus on the Hamiltonian employing $\mathcal{H}_{\rm int}^{(1)}$ rather than both $\mathcal{H}_{\rm int}^{(1)}$ and $\mathcal{H}_{\rm int}^{(2)}$. As shown in Table 1, the magnitude of $J_{ij}$ is smaller than $U_{ij}$, and it is considered that the Hamiltonian with $\mathcal{H}_{\rm int}^{(2)}$ exhibits a similar trend to that with $\mathcal{H}_{\rm int}^{(1)}$.

For further confirmation of the convergence behavior, we also vary the cutoff parameter set by comparing the $(a,E_{\textrm{cut}}^{\psi},E_{\textrm{cut}}^{\varepsilon})=(17.0,30.0,3.0)$ condition with the $(13.0,50.0,5.0)$ condition, where the units of the lattice constant and cutoff parameters are Å and Ry, respectively. Due to our computational resources, it is difficult to increase the value of cutoff parameters while maintaining the lattice constant $a$ at 17.0 Å. However, it is important to note that this lattice constant is sufficiently large to treat the molecule as isolated, ensuring the relevance and accuracy of our results within this computational framework.

Figure 3 shows the $N_{\rm band}$ dependency of our systems and the result of extrapolation. The state characterization is conducted in the way shown in Appendix A.

As a result, the excitation energy data varying on $N_{\rm band}$ are successfully fitted using Eq. (15), as shown in Figure 3. Comparing the two $(a,E_{\textrm{cut}}^{\psi},E_{\textrm{cut}}^{\varepsilon})$ conditions, the shapes, the intersection point, and the convergence behavior to the large $N_{\textrm{band}}$ limit of the curves fitted to the $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ states are similar. It suggests that our excitation energy results do not vary significantly with the setting of these hyperparameters.

It is found that there is a qualitative difference in the $N_{\rm band}$ dependence among the molecules. In the ethylene (2e, 2o) case, as shown in Figure 3 (a), the excitation energies to the $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ states decrease in a similar way as $N_{\rm band}$ increases. Conversely, in the other cases, the excitation energy to the $1{}^{1}{\rm B_{u}}$ state decreases, but that to the $2{}^{1}{\rm A_{g}}$ state increases as $N_{\rm band}$ increases, as shown in Figure 3 (b)–(d). In the small $N_{\rm band}$ region, lacking the dynamical correlation, the $2{}^{1}{\rm A_{g}}$ state is lower than the $1{}^{1}{\rm B_{u}}$ state. As $N_{\rm band}$ increases, the order is reversed, and the excitation energies become saturated eventually. It is crucial to take a sufficiently large $N_{\rm band}$ to reasonably discuss the energetic order of the $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ states and its gap width within the framework of cRPA.

In Table 2, we compare the extrapolated $\Delta E_{\infty}$ values with the $\Delta E$ values of the models, which employ $\mathcal{H}_{\rm int}=\mathcal{H}_{\rm int}^{(1)}$ and $N_{\textrm{band}}=N_{\textrm{band}}^{\textrm{max}}$.

These values are similar in each molecule, and it is considered that the models at the $N_{\rm band}^{\rm max}$ points are constructed by incorporating the electron correlation effects almost the same as the effects in the continuum limit. The values of $U_{ij}$ and $J_{ij}$ at the $N_{\rm band}^{\rm max}$ points are shown in Table 1.

We compare the excitation energies of our models with the results of other quantum chemical calculations and experiments shown in Table 3.

Hereafter, the models that employ $\mathcal{H}_{\rm int}^{(1)}$ and $\mathcal{H}_{\rm int}^{(2)}$ together with the parameters in Table 1 are referred to as Model 1 and Model 2, respectively.

Overall, Model 2 practically yields a better result than Model 1. Model 1 tends to underestimate the excitation energies compared to Model 2, which is considered to be due to the exchange interactions in Model 1.

The excitation energies of Models 1 and 2 become smaller as the conjugation length becomes longer. This tendency is similar to the results of quantum chemical calculations and experiments. In particular, the excitation energies of Model 2 are in good agreement with the experimental values in either molecule. We employ the (4e, 4o) and (6e, 6o) active spaces for hexatriene, and the former model has the excitation energies closer to the experimental values. This suggests that a larger active space does not necessarily lead to a better result in the construction of the ab initio downfolding model.

The magnitude of the excitation energies to the lower excited states of polyenes has been an important research topic in traditional quantum chemistry Schulten and Karplus (1972); Tavan and Schulten (1986); Nakayama et al. (1998); Kurashige et al. (2004). In butadiene and hexatriene, it is known that the singly excited $1{}^{1}{\rm B_{u}}$ and doubly excited $2{}^{1}{\rm A_{g}}$ states are energetically close to each other. It can be said that the gap of our models is small — the excitation energy gaps of Model 2 of butadiene and hexatriene are 0.16 eV and 0.31 eV, respectively. These narrow gaps are attributed to the inclusion of dynamical electron correlation from outside of the active space. Comparison between CASCI and CASPT2 results also indicates that the dynamical electron correlation reduces the energy gap between $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$.

We compare the results of our models with those of CC3 and CASPT2 calculations. In our CC3 calculations, the excitation energies get close to the models’ and experimental results in each molecule. The states of CC3 are characterized using the oscillator strengths and the molecular orbital shapes related to the major single and double amplitudes. The excitation energy gaps of butadiene and hexatriene are 0.43 eV and 0.38 eV, comparable to those in Model 2. In our CASPT2 calculations, the excitation energies to the $1{}^{1}{\rm B_{u}}$ state capture the trend of the results of Model 2, experiment, and CC3. The excitation energy gaps of butadiene and hexatriene are 0.31 eV and $-$0.35 eV in CASPT2(4e, 4o) — The excitation energy to the $2{}^{1}{\rm A_{g}}$ state is lower than that of the $1{}^{1}{\rm B_{u}}$ state in hexatriene.

Highly precise quantum chemical calculations have provided insight into the energetic order of the $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ states. It has been known that the $1{}^{1}{\rm B_{u}}$ state is lower in butadiene Watson and Chan (2012); Chien et al. (2018). Our models of butadiene reproduce the order of the $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ states with a reasonably small energy gap.

In hexatriene, high-level quantum chemical calculations show that the $2{}^{1}{\rm A_{g}}$ state is much closer to $1{}^{1}{\rm B_{u}}$ in energy and slightly lower Nakayama et al. (1998); Kurashige et al. (2004). For instance, the result of the semi-stochastic heat-bath CI (SHCI) reported by Chien et al. Chien et al. (2018) is shown in Table 3, and the $2{}^{1}{\rm A_{g}}$ state is slightly lower than and accidentally degenerates to the $1{}^{1}{\rm B_{u}}$ state in hexatriene. The state order of our models is reversed to the SHCI result, yet the energy gap remains narrow. The order of the CC3 hexatriene result is also reversed in the order of these states. Note that the state order from experiments is similar to the results of our models and CC3. There seems to be room for a more detailed discussion about the energy gap in hexatriene. For example, the geometry optimization calculation of the ground state at a much higher level might refine our discussion.

We also compare the result of our models with that of the CI calculation using the PPP model, termed PPP-CI Tavan and Schulten (1986). The single excitation energy of PPP-CI for butadiene closely matches the experimental value, whereas the double excitation energy is smaller than the single excitation energy. On the other hand, our models reproduce the relationship between $1{}^{1}{\rm B_{u}}$ and $2{}^{1}{\rm A_{g}}$ in butadiene well. For hexatriene, PPP-CI qualitatively reproduced the order that the $2{}^{1}{\rm A_{g}}$ state is lower than the $1{}^{1}{\rm B_{u}}$, but the gap seems large ($-$0.69 eV). In our models, the $1{}^{1}{\rm B_{u}}$ state is lower than the $2{}^{1}{\rm A_{g}}$ state, albeit with a small gap.

Finally, we discuss the TD-DFT result. TD-DFT reproduces the excitation energy to the $1{}^{1}{\rm B_{u}}$ state well, but in the butadiene and hexatriene cases, the excitation energy to the $2{}^{1}{\rm A_{g}}$ state is about 1 eV larger than that to the $1{}^{1}{\rm B_{u}}$ state. This feature differs from our models and other established quantum chemical calculations. Conventional TD-DFT is known to encounter challenges in describing double electron excitations. While it can reasonably predict the single excitation to the $1{}^{1}{\rm B_{u}}$ state, TD-DFT struggles with representing the double excitation to the $2{}^{1}{\rm A_{g}}$ state, as highlighted in Ref. Cave et al. (2004).

Here, we analyze the extent to which the extended Hubbard model Hamiltonian can contribute to efficient quantum computation. The $L^{1}$-norm $\lambda$ and the number of terms $N_{\rm term}$ of the fermion-to-qubit mapped Hamiltonians are shown in Table 4.

Our model Hamiltonians have the $L^{1}$-norm value comparable to the active space Hamiltonian in each molecule. Since Model 2 does not include the exchange interaction terms, the $\lambda$ values of Model 2 are reduced compared to Model 1, yet their differences are negligible. On the contrary, the first-principles Hamiltonians of the full space have a huge $L^{1}$-norm value. This contrast clearly shows that the quantum computing cost of our extended Hubbard Hamiltonian is similar to that of the active space Hamiltonian and significantly cheaper than that of the full space first-principles Hamiltonian. Recalling that our models include dynamical electron correlations from outside of the active space, as discussed in Section IV.3, they have a significant advantage in the efficiency of quantum computation compared to the active space Hamiltonians, which do not include such correlations.

Since the required active space is very small for these targeted molecules, we could not show numerically that the extended Hubbard models have the advantage regarding the $L^{1}$-norm value against the corresponding active space Hamiltonian. The $L^{1}$-norm values are, however, expected to decrease in much larger active space cases because the extended Hubbard and active space Hamiltonians have the square and quartic number of the electron-electron interaction terms, respectively. This effect may be seen in the hexatriene (6e, 6o) case and is expected to be more remarkable in a chemical system requiring a larger active space, e.g., larger $\pi$-conjugated systems or transition metal complexes.

Table 4 also shows the results of ER localization of active space Hamiltonians in the butadiene and hexatriene cases. The maximally localized Wannier function is localized orbital rather than canonical molecular orbital (CMO). Compared to the canonical orbital cases, the $L^{1}$-norm values in the ER-localized basis have slightly increased. Koridon et al. have previously reported that localized orbitals can decrease the $L^{1}$-norm values Koridon et al. (2021), but our numerical investigation does not clearly show this trend. It can be just because the necessary active space for our system is very small.

Regarding $N_{\textrm{term}}$, Model 2 is the smallest within the Hamiltonians of the same active space in either molecule. The effect of neglecting the exchange terms appears in comparing Model 1 and Model 2, even though our systems are small. The maximally localized Wannier functions do not reflect the high symmetry of the molecules. However, the extended Hubbard model Hamiltonians have a significant advantage in $N_{\textrm{term}}$, whose second-quantized representation does not have the electron-electron interaction terms of the three and four kinds of distinct indices. It is expected to facilitate the execution of quantum computations.

On the contrary, in the cases of active space Hamiltonians in an ER-localized basis, the number of terms $N_{\rm term}$ has almost doubled compared to that in a CMO basis. Since the all-trans polyenes are highly symmetrized, the fermion-to-qubit mapped active space Hamiltonian in a CMO basis is sparse, and the number of terms is reduced. Orbital localization leads to the lower symmetry of the orbitals, and it can be considered that the number of the terms increases.

The small number of terms in our models is considered to have advantages in quantum computations. For example, it can be robust against errors from the noises of the noisy intermediate-scale quantum devices Preskill (2018) because the number of measurements to evaluate the expectation values of Pauli operators is expected to become smaller. For quantum algorithms such as quantum phase estimation, the quantum circuits become smaller and shallower.

We further apply our extended Hubbard models to the classical simulation of quantum computations. We estimate the expectation values of the models and the ER-localized active space Hamiltonians by sampling. The quantum circuits were prepared for the ground and excited states using VQD calculations.

The VQD algorithm is an extension of VQE for calculating excited states, proposed by Higgott et al. Higgott et al. (2019). VQD’s cost function $F({\bm{\theta}_{k}})$ is defined as

$\displaystyle F({\bm{\theta}_{k}}):=\braket{\Psi({\bm{\theta}_{k}})}{\mathcal{H}}{\Psi({\bm{\theta}_{k}})}+\sum_{i=0}^{k-1}\beta_{i}|\braket{\Psi({\bm{\theta}_{k}})}{\Psi({\bm{\theta}_{i}})}|^{2},$

(16)

where $|\Psi({\bm{\theta}_{k}})\rangle$ is the $k$-th target excited state that can be obtained by optimizing the quantum circuit parameters ${\bm{\theta}_{k}}$. The first term of the cost function corresponds to the expectation value of the energy of the $k$-th excited state. The second term is the penalty term that is introduced as the sum of the overlap between the $k$-th state and $i$-th ansatz states ($i=0,\cdots,k-1$). The parameter $\beta_{i}$ is the weight factor of each overlap term.

Table LABEL:tab:smpl shows the estimated energy differences and the errors of the (4e, 4o) models and ER-localized active space Hamiltonians of butadiene and hexatriene.

The energy difference is between the energy estimated by sampling and the exact ground-state energy. The energy is estimated by repeating the $10^{4}$ shots energy calculation $10^{4}$ times.

The estimated energy differences using our models are in good agreement with those of exact diagonalization. Those using the ER-localized active space Hamiltonians are similarly estimated well but overestimate the excitation energies due to the absence of electron correlation related to the outside of the active space. The errors are similar for both our models and ER-localized active space Hamiltonians. Hence, we validate the applicability of the ab initio downfolded model to compute excited states using quantum devices through numerical calculations using our models.