# Orbital transformations to reduce the 1-norm of the electronic structure Hamiltonian for quantum computing applications

Emiel Koridon, Saad Yalouz, Bruno Senjean, Francesco Buda, Thomas O'Brien, Lucas Visscher

<figure><img src="image_1.png" /></figure>

## To cite this version:

Emiel Koridon, Saad Yalouz, Bruno Senjean, Francesco Buda, Thomas O'Brien, et al., Orbital transformations to reduce the 1-norm of the electronic structure Hamiltonian for quantum computing applications. Physical Review Research, 2021, 3 (3), 10.1103/PhysRevResearch.3.033127. [https://doi.org/10.1103/PhysRevResearch.3.033127](https://doi.org/10.1103/PhysRevResearch.3.033127)

HAL Id: hal-03365850

https://hal.science/hal-03365850

Submitted on 1 Jun 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Orbital transformations to reduce the 1-norm of the electronic structure
Hamiltonian for quantum computing applications

Emiel Koridon●,¹,² Saad Yalouz,²,¹,\* Bruno Senjean,³,¹,† Francesco Buda,⁴
Thomas E. O’Brien,⁵,¹ and Lucas Visscher²

¹Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, Netherlands

²Theoretical Chemistry, Vrije Universiteit, De Boelelaan 1083, NL-1081 HV, Amsterdam, Netherlands

³ICGM, Univ Montpellier, CNRS, ENSCM, 34090 Montpellier Cedex 5 - France

⁴Leiden Institute of Chemistry, Leiden University, Einsteinweg 55, P.O. Box 9502, 2300 RA Leiden, Netherlands

⁵Google Research, 80636 Munich, Germany

<figure><img src="image_3.png" /></figure>

(Received 18 April 2021; accepted 7 July 2021; published 6 August 2021)

Reducing the complexity of quantum algorithms to treat quantum chemistry problems is essential to demonstrate an eventual quantum advantage of noisy-intermediate scale quantum devices over their classical counterpart. Significant improvements have been made recently to simulate the time-evolution operator $U(t) = e^{i\hat{H}t}$, where $\hat{H}$ is the electronic structure Hamiltonian, or to simulate $\hat{H}$ directly (when written as a linear combination of unitaries) by using block encoding or qubitization techniques. A fundamental measure quantifying the practical implementation complexity of these quantum algorithms is the so-called 1-norm of the qubit representation of the Hamiltonian, which can be reduced by writing the Hamiltonian in factorized or tensor-hypercontracted forms, for instance. In this paper, we investigate the effect of classical preoptimization of the electronic structure Hamiltonian representation, via single-particle basis transformation, on the 1-norm. Specifically, we employ several localization schemes and benchmark the 1-norm of several systems of different sizes (number of atoms and active space sizes). We also derive a formula for the 1-norm as a function of the electronic integrals and use this quantity as a cost function for an orbital-optimization scheme that improves over localization schemes. This paper gives more insights about the importance of the 1-norm in quantum computing for quantum chemistry and provides simple ways of decreasing its value to reduce the complexity of quantum algorithms.

DOI: 10.1103/PhysRevResearch.3.033127

I. INTRODUCTION

Quantum chemistry has been identified as the killer ap-
plication of quantum computers, which promises to solve
problems of high industrial impact that are not tractable for
their classical counterparts [1–4]. Unfortunately, quantum de-
vices are noisy and only shallow circuits can be implemented
with a relatively accurate fidelity. Hence, even though we
know how to solve the electronic structure problem with,
for instance, the quantum phase estimation [5–9] algorithm,
this remains out of reach in practice and one has to come
up with new algorithms that can actually be run within this
noisy intermediate-scale quantum (NISQ) era [10]. Reducing
circuit depth has been achieved by interfacing quantum and
classical devices, thus leading to hybrid quantum-classical
algorithms such as the variational quantum eigensolver (VQE)

[11]. Since then, variational algorithms have been successfully
applied to the estimation of ground-state (see Refs. [12–15]
and references therein) and excited-state energies [16–26], as
well as molecular properties [27–30].

However, it remains unclear if these algorithms can provide
a clear quantum advantage in the long run, especially due
to the difficulty in circuit optimization [31] and to the large
overhead in the number of measurements required to achieve
sufficient accuracy [32,33], though significant progress has
been made recently [33–45]. This led to the development of
strategies for reducing the resources required to implement
quantum algorithms (such as gate complexity) for both the
fault-tolerant and NISQ era. For simulating the time-evolution
operator $U(t) = e^{i\hat{H}t}$, where $\hat{H}$ is the electronic Hamilto-
nian, Campbell has shown that using a randomized compiler
is better suited than the Trotter-Suzuki decomposition, the
so-called quantum stochastic drift protocol (qDRIFT) [46].
One can also directly simulate the Hamiltonian by using a
linear combination of unitaries (LCU), generalized by the
block encoding or qubitization formalism [47–50]. Most of
these algorithms have a gate complexity that scales with re-
spect to the parameter $\lambda_Q = \sum_i |h_i|$, where $\hat{H} = \sum_i h_i \hat{P}_i$ is
the qubit Hamiltonian and $\hat{P}_i$ are Pauli strings [51]. Low-
ering the value of $\lambda_Q$—usually referred to as the 1-norm
of the Hamiltonian—has a significant impact on quantum

Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.

<footnote>*yalouzsaad@gmail.com</footnote>

<footnote>†bruno.senjean@umontpellier.fr</footnote>
algorithms, and techniques such as double-factorization [4], tensor hypercontraction [52], and n-representability constraints [33] have proven successful in this respect.

In this paper, we investigate how the value of the 1-norm and its scaling with respect to the number of orbitals can be reduced by single-particle basis rotations. As a starting point, we consider the use of localized orbitals that are commonly available in most quantum chemistry codes. Generating such orbitals presents a simple and effortless way of reducing gate complexity by classical preoptimization of the Hamiltonian representation.

We show that the use of localized orbitals has a significant impact on the 1-norm of all systems studied, from simple one-dimensional hydrogen chains (also studied in Refs. [33,52]) to much more complex organic and inorganic molecules. Additionally, we connect the expressions of the 1-norm before and after the fermion-to-qubit mapping. This expression is used to define a cost function for single-particle basis rotations, and is shown to reduce the value of $\lambda_Q$ even more than standard localization schemes.

The paper is organized as follows. After introducing the electronic structure problem in Sec. II A, the impact of the 1-norm on quantum algorithms is detailed in Sec. II B. Then, we review the localization schemes applied in this paper in Sec. II C followed by our so-called 1-norm orbital-optimization (OO) method in Sec. II D. Finally, preceded by the computational details in Sec. III, the scaling of the 1-norm when increasing the number of atoms is investigated in Sec. IV A for hydrogen and alkane chains, followed by a benchmark on several systems in Sec. IV B and a study of the effect of increasing the active space size in Sec. IV C. Conclusions and perspectives are given in Sec. V.

## II. THEORY

### A. Electronic structure Hamiltonian in second quantization

In quantum chemistry, the second quantization formalism is usually employed to describe the electronic properties of molecules. Under the Born-Oppenheimer approximation, $N_e$ electrons can rearrange around $N_a$ nuclei by occupying a restricted set of $N$ spatial molecular orbitals (MOs). The orthonormal MO basis $\{\phi_p(\mathbf{r})\}$ is built from linear combination of atomic orbitals (AOs) and is usually obtained from an inexpensive mean-field calculation, e.g., with the Hartree-Fock (HF) method [53]. In this basis, the molecular Hamiltonian can be expressed in a spin-free form that reads (in atomic units)

$$ \hat{\mathcal{H}} = \sum_{pq}^{N} h_{pq} \hat{E}_{pq} + \frac{1}{2} \sum_{pqrs}^{N} g_{pqrs} \hat{e}_{pqrs}, \quad (1) $$

where $\hat{E}_{pq} = \sum_{\sigma} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma}$ and $\hat{e}_{pqrs} = \sum_{\sigma, \tau} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{\tau\sigma}^{\dagger} \hat{a}_{s\tau} \hat{a}_{q\sigma}$ are the one- and two-body spin-free operators. In Eq. (1), the coefficients

$$ h_{pq} = \int \phi_p^*(\mathbf{r}_1) \left( -\frac{1}{2} \nabla_{\mathbf{r}_1}^2 - \sum_{A=1}^{N_a} \frac{Z_A}{\mathbf{r}_{1A}} \right) \phi_q(\mathbf{r}_1) d\mathbf{r}_1 \quad (2) $$

are the so-called one-electron integrals which encode, for each individual electron, the associated kinetic energy and

Coulombic interaction with the $N_a$ nuclei of the molecule (with $\mathbf{r}_{iA} = |\mathbf{r}_i - \mathbf{R}_A|$ and $Z_A$ the distance between electron i and nuclei A, and the atomic number of atom A). The coefficients

$$ g_{pqrs} = (pq|rs) = \iint \phi_p^*(\mathbf{r}_1) \phi_r^*(\mathbf{r}_2) \frac{1}{\mathbf{r}_{12}} \phi_q(\mathbf{r}_1) \phi_s(\mathbf{r}_2) d\mathbf{r}_1 d\mathbf{r}_2 \quad (3) $$

are the so-called two-electron integrals encoding the Coulombic repulsion between each pair of electrons (with $\mathbf{r}_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ the distance between electrons i and j). Solving the electronic structure problem consists of solving the time-independent Schrödinger equation

$$ \hat{\mathcal{H}} |\Psi_{\ell}\rangle = E_{\ell} |\Psi_{\ell}\rangle, \quad (4) $$

with $|\Psi_{\ell}\rangle$ an electronic eigenstate with corresponding energies $E_{\ell}$. This is only possible in the case of small molecules and small basis sets due to the exponential scaling of the computational cost with respect to the size of the system. A commonly employed approach is the so-called active space approximation which divides the MO space into three parts: the core (frozen occupied orbitals), the active, and the virtual spaces (deleted unoccupied orbitals), such that only the electrons inside the active space are treated explicitly. This approximation will be used throughout this paper and is equivalent to finding the eigenstates of an effective Hamiltonian $\hat{\mathcal{H}}^{\text{FC}}$ also called the frozen-core Hamiltonian (see Appendix A).

The electronic Hamiltonian can be mapped onto an appropriate representation for quantum computers by doing a fermion-to-qubit transformation (such as Jordan-Wigner [51]), resulting in a linear combination of Pauli strings $\hat{P}_j \in \{I, X, Y, Z\}^{\otimes N}$,

$$ \hat{\mathcal{H}}_Q = \sum_j^S h_j \hat{P}_j. \quad (5) $$

Here, $S$ denotes the sparsity of the Hamiltonian and generally scales as $O(N^4)$, but sometimes $O(N^2)$ for a sufficiently large system and a localized basis [54].

### B. The 1-norm in quantum computing

The term 1-norm in quantum computing refers to a norm induced on a Hilbert-space operator by its decomposition as a sum of simpler terms. To calculate a 1-norm, we write an operator $\hat{\mathcal{H}}$ (e.g., a Hamiltonian) in the form

$$ \hat{\mathcal{H}} = \sum_j b_j \hat{B}_j, \quad (6) $$

where the $\hat{B}_j$ are operators on $\mathbb{C}^{2N}$ and $b_j$ are complex numbers. Typically, the operators $\hat{B}_j$ are chosen to be either unitary ($\hat{B}_j^\dagger \hat{B}_j = 1$) or unital ($\hat{B}_j^\dagger \hat{B}_j < 1$). We define the 1-norm of $\hat{\mathcal{H}}$ (induced by this decomposition) to be the 1-norm on the vector $\vec{b}$ (with components $b_j$):

$$ \lambda_B = \lambda_B(\hat{\mathcal{H}}) = \sum_j |b_j|. \quad (7) $$
(It is common to suppress the dependence on the Hamiltonian $\hat{\mathcal{H}}$ when $\lambda_B$ is unambiguous.) For example, under the above decomposition into Pauli operators [Eq. (5)], we have

$$ \lambda_Q \equiv ||\vec{h}||_1 = \sum_j |h_j|. \qquad (8) $$

In general, $\lambda_B$ is a norm on the space spanned by the set $\{\hat{B}_j\}$ as long as this set is linearly independent, as the mapping $\hat{\mathcal{H}} \leftrightarrow \vec{h}$ is then bijective and linear. The Pauli operators $\hat{P}_j$ used above make a natural choice for $\hat{B}_j$ as they span the set of $2^{2N} \times 2^{2N}$ matrices, and so induce a norm for any operators on the Hilbert space. Pauli operators are also unitary, and may be easily implemented in a quantum circuit [57], making them a common choice for Hamiltonian decompositions. However, unlike other operator size measures, 1-norms may be heavily optimized by a change of operator basis (whereas e.g., the 2-norm or infinity-norm are both invariant under unitary rotation). As operator 1-norms play a role in determining the cost of many quantum computing algorithms, this makes them a key target of study in quantum algorithm research.

A versatile method for implementing an arbitrary (non-Hermitian) operator on a quantum device is the LCU method [47], which involves a decomposition over a set $\{\hat{B}_j\}$ of strictly unitary operators. As $\hat{\mathcal{H}}$ is typically not unitary, the original LCU technique requires postselection, but this forms the basis for unitary methods such as qubitization [49]. These, in turn, underpin many proposals for quantum computing on a quantum computer [4,52,56], and the dependence on $\lambda_B$ is pulled directly through. We give a brief description here of the key identity in LCU methods to show where the dependence on $\lambda_B$ appears. LCU methods all require a control register with states $|j\rangle$ that encode the index to the operator $\hat{B}_j$. (One need only encode those operators $\hat{B}_j$ that are relevant to the target $\hat{\mathcal{H}}$.) Then, these methods encode the coefficients of $\hat{\mathcal{H}}$ on the register by preparing the state $\frac{1}{\sqrt{\lambda_B}} \sum_j \sqrt{\hat{b}_j|j\rangle}$. Then, consider implementing the joint unitary $\sum_j |j\rangle\langle j| \otimes \hat{B}_j$ (that is, applying the unitary operator $\hat{B}_j$ to a system conditional on the control register being in the state $j$). Postselecting the control register returning to the initial state performs the LCU operation to the register:

$$ \frac{1}{\sqrt{\lambda_B}} \sum_j \sqrt{\hat{b}_j \langle j|} \left[ \sum_j |j\rangle \langle j| \hat{B}_j \right] \frac{1}{\sqrt{\lambda_B}} \sum_j \sqrt{\hat{b}_j |j\rangle} \\ = \frac{1}{\lambda_B} \sum_j b_j \hat{B}_j = \frac{1}{\lambda_B} \hat{\mathcal{H}}. \qquad (9) $$

We see here that $\lambda_B$ emerges naturally as the normalization constant of the LCU unitary. This dependence is passed on to any LCU-based method, suggesting that minimizing $\lambda_B$ is key to optimizing any such techniques.

The 1-norm can also play a role in methods where Hamiltonian simulation is implemented stochastically. For instance, in the qDRIFT [46] method, Hermitian terms $\hat{B}_j$ in the Hamiltonian are chosen at random, weighted by $|b_j|/\lambda_B$, and exponentiated on a system as $e^{i\tau\hat{B}_j}$ ($\lambda_B$ appears here immediately as the normalization of the probability distribution). To

first order, this implements the channel on a state $\hat{\rho}$ [46],

$$ \mathcal{E}(\hat{\rho}) = \hat{\rho} + i \sum_j \frac{b_j \tau}{\lambda_B} (\hat{B}_j \hat{\rho} - \hat{\rho} \hat{B}_j) + O(\tau^2), \qquad (10) $$

which can be observed to be the unitary evolution of $\hat{\rho}$ under $e^{i\hat{\mathcal{H}}dt}$ for $dt = \tau/\lambda_B$. To approximate time evolution for time $t$ to error $\epsilon$ then requires repeating this process $O(2\lambda_B^2 t^2/\epsilon)$ times. This has been extended upon recently [58], yielding a more complicated $\lambda_B$ dependence.

A final situation where the 1-norm of an operator plays a role is in tomography. One may, in principle, measure the expectation value of an arbitrary Hermitian operator $\hat{\mathcal{H}}$ on a state $\rho$ by repeatedly preparing $\hat{\rho}$, rotating into the eigenbasis of $\hat{\mathcal{H}}$, and reading out the device. The variance in averaging $M$ shots of preparation and measurement is directly given by the variance of $\hat{\mathcal{H}}$ on $\hat{\rho}$:

$$ \mathrm{Var}[\langle\hat{\mathcal{H}}\rangle] = \frac{1}{M} [\langle\hat{\mathcal{H}}^2\rangle - \langle\hat{\mathcal{H}}\rangle^2]. \qquad (11) $$

However, in practice (especially in the NISQ era), arbitrary rotations into an unknown basis may be costly. As expectation values are linear, given a decomposition of $\hat{\mathcal{H}}$ over $\{\hat{B}_j\}$, one may instead perform $M_j$ different preparations of $\hat{\rho}$ and measurements in the basis of each $\hat{B}_j$ and sum the result,

$$ \langle\hat{\mathcal{H}}\rangle = \sum_j b_j \langle\hat{B}_j\rangle, \qquad (12) $$

with variance

$$ \mathrm{Var}[\langle\hat{\mathcal{H}}\rangle] = \sum_j \frac{|b_j|^2}{M_j} [\langle\hat{B}_j^2\rangle - \langle\hat{B}_j\rangle^2] =: \sum_j \frac{|b_j|^2}{M_j} V_j. \qquad (13) $$

Typically, the $V_j$ are not known in advance, but if $\hat{B}_j$ are chosen to be unital (by a suitable rescaling of $b_j$), then one can bound $V_j < 1$. One can confirm (e.g., via the use of Lagrange multipliers [33]) that the least number of measurements required to bound $\mathrm{Var}[\langle\hat{\mathcal{H}}\rangle]$ below some error $\epsilon^2$ can be found by choosing $M_j = M|b_j|/\lambda_B$, and that this yields a total number of measurements

$$ M = \sum_j M_j = \frac{\lambda_B^2}{\epsilon^2}, \qquad (14) $$

first mentioned by [32]. State tomography is essential in hybrid quantum-classical algorithms such as a VQE. In the VQE, the Hamiltonian is usually expressed as a linear combination of Pauli operators, meaning $\lambda_B = \lambda_Q$. As the number of measurements is one of the biggest overheads in the practical implementation of VQE, lowering the 1-norm would be highly beneficial. One way to do so is to add constraints to the Hamiltonian based on the fermionic n-representability conditions, which has led to a reduction of up to one order of magnitude for hydrogen chains and the H$_4$ ring molecule [33]. Another hybrid algorithm that would benefit for a decrease in the 1-norm is the constrained-VQE algorithm, as this would decrease the higher-bound used in the penalty term, thus improving convergence [59].

As the above considerations imply, many of the most competitive algorithms have a strong dependence on 1-norms for various choices of the basis decomposition $\{\hat{B}_j\}$, as shown in
<table><caption>TABLE I. Lowest asymptotic scaling of different quantum algorithms involving the $\lambda$ parameters. $\lambda$ denotes the 1-norm of the qubit Hamiltonian in Eq. (15), $\lambda_{SF}$ (denoted $\lambda_W$ in Ref. [50]) the one of the singly factorized fermionic Hamiltonian, $\lambda_{DF}$ the one of the doubly factorized Hamiltonian and $\lambda_\zeta$ the one of the nonorthogonal tensor hypercontracted Hamiltonian representation. S is the sparsity of the electronic Hamiltonian and $\Xi$ the average rank of the second tensor factorization [4,56]. See text for further details.</caption><thead><tr><th>Method</th><th>Toffoli/T complexity</th></tr></thead><tbody><tr><td>Database method [55]</td><td>O(N<sup>4</sup>λ/ε)</td></tr><tr><td>qDRIFT [46]</td><td>O(λ<sup>2</sup>/ε<sup>2</sup>)</td></tr><tr><td>Qubitization (sparse method) [50]</td><td>O((N + √S)λ/ε)</td></tr><tr><td>Qubitization (single factorization) [50]</td><td>O(N<sup>3/2</sup>λ<sub>SF</sub>/ε)</td></tr><tr><td>Qubitization (double factorization) [4]</td><td>O(Nλ<sub>DF</sub>√Ξ/ε)</td></tr><tr><td>Tensor hypercontraction [52]</td><td>Ō(Nλ<sub>ζ</sub>/ε)</td></tr><tr><td>VQE</td><td>No. of measurements</td></tr><tr><td>State tomography [11,12,32]</td><td>M ≈ λ<sup>2</sup>/ε<sup>2</sup></td></tr><tr><td>Basis rotation grouping [41]</td><td>M ≈ λ'<sub>SF</sub>/ε<sup>2</sup></td></tr></tbody></table>

Table I. The 1-norm induced by a Pauli decomposition, $\lambda_Q$ [Eq. (8)], scales with the number of orbitals somewhere in between $O(N)$ and $O(N^3)$. This depends on the system, and whether $N$ increases because the number of atoms increases for a fixed number of basis functions per atom or because the number of basis functions per atom increases while fixing the number of atoms. As far as we are aware, this scaling has only been numerically benchmarked for the $H_4$ ring Hamiltonian [33,52] and one-dimensional hydrogen chains [52], which remain quite far from realistic chemistry problems encountered in enzymatic [1] and catalytic reactions [4]. In this paper, we perform a similar benchmark on a larger set of organic and inorganic molecules, from hydrogen and carbon chains to the FeMoco and ruthenium complexes.

The 1-norm used in many of the aforementioned algorithms in Table I takes the following form [50,52]:

$$ \lambda = \lambda_T + \lambda_V, \quad (15) $$

where

$$ \lambda_T = \sum_{pq}^{N} \left| h_{pq} + \frac{1}{2} \sum_{r}^{N} (2g_{pqrr} - g_{prrq}) \right| \quad (16) $$

and

$$ \lambda_V = \frac{1}{2} \sum_{pqrs}^{N} |g_{pqrs}|, \quad (17) $$

with $\lambda_V \gg \lambda_T$. Note that the difference between Eqs. (16) and (17) and Eqs. (A10) of Ref. [52] comes from a different convention when writing the electronic Hamiltonian in Eq. (1) and the fact that $N$ refers to spatial orbitals in this work instead of spin-orbitals. This norm of the Hamiltonian expressed as a LCU is used in the database method of Babbush et al. [55], the qubitized simulation by the sparse method of Berry et al. [50] (further improved by Lee et al. [52]) and the qDRIFT protocol of Campbell [46].

Further attempts have been made to reduce the number of terms in the Hamiltonian. For instance, one may perform

a low rank decomposition of the Coulomb operator $\hat{V}$ such that the 1-norm of the—now singly factorized—Hamiltonian reads [50,52]

$$ \lambda_{SF} = \frac{1}{4} \sum_{\ell=1}^{L} \left( \sum_{p,q}^{N} |W_{pq}^{(\ell)}| \right)^2, \quad (18) $$

where $W_{pq}^{(\ell)}$ are obtained from a Cholesky decomposition of the $g_{pqrs}$ tensor. Note that $\lambda_{SF}$ is a higher-bound of $\lambda_V$, and so this decomposition tends to increase the 1-norm. To improve over this low-rank factorization, one can write the Hamiltonian in a doubly-factorized form [4] by rotating the single-particle basis to the eigenbasis of the Cholesky vectors, such that the corresponding 1-norm $\lambda_{DF}$ is much lower than after a single factorization. One can also use the tensor hyper-contraction representation of the Hamiltonian [52]. However, applying the qubitization algorithm on this Hamiltonian directly is not efficient, as its associated 1-norm $\lambda_{THC}$ is even larger than $\lambda_V$. To bypass this issue, Lee et al. provided a diagonal form of the Coulomb operator by projection into an expanded nonorthogonal single-particle basis, thus leading to a nonorthogonal THC representation of the Hamiltonian with an associated 1-norm $\lambda_\zeta$ that scales better than any prior algorithm [52].

As readily seen in this section, the $\lambda$ norm depends on how you represent your Hamiltonian. However, there is a unique way to write the Hamiltonian as a sum of unique Pauli strings, or equivalently as a sum of unique products of Majorana operators (see Appendix C). Doing so allows us to express the qubit 1-norm $\lambda_Q$ in Eq. (8) as a function of the electronic integrals,

$$ \lambda_Q = \lambda_C + \lambda_T + \lambda'_V, \quad (19) $$

where

$$ \lambda_C = \left| \sum_{p}^{N} h_{pp} + \frac{1}{2} \sum_{pr}^{N} g_{purr} - \frac{1}{4} \sum_{pr}^{N} g_{prrp} \right|, \quad (20) $$

$$ \lambda_T = \sum_{pq}^{N} \left| h_{pq} + \sum_{r}^{N} g_{pqrr} - \frac{1}{2} \sum_{r}^{N} g_{prrq} \right|, \quad (21) $$

$$ \lambda'_{V} = \frac{1}{2} \sum_{p>r,s<q}^{N} |g_{pqrs} - g_{psrq}| + \frac{1}{4} \sum_{pqrs}^{N} |g_{pqrs}|. \quad (22) $$

The first term $\lambda_C$ corresponds to the absolute value of the coefficient of the identity term that emerges when rearranging the Hamiltonian in terms of unique Majorana operators, see Eq. (C8). It is invariant under orbital rotations and can be added to the energy as a classical constant together with the nuclear repulsion energy and, if one employs the frozen core approximation (see Appendix A), the mean-field energy of the frozen core. Thus, this term (apart from being small compared to $\lambda_T$ and $\lambda'_V$) will not be important for quantum algorithms and we will leave it out of our results, redefining the 1-norm as $\lambda_Q = \lambda_T + \lambda'_V$. $\lambda_T$ represents the absolute values of the coefficients of the quadratic term in Majorana operators and $\lambda'_V$ of the quartic term. Notice that $\lambda'_V$ is slightly different from $\lambda_V$ of Ref. [52] also given in Eq. (17) (which added a slight correction to the one of Ref. [50]). In this paper, we add another correction that comes from the swapping of two
Majorana operators. We confirm that this is the correct form by directly comparing $\lambda_Q$ with the norm of the qubit Hamiltonian after doing a qubit-to-fermion transformation in Eq. (8). Note that $\lambda'_V \le \lambda_V$ (see Appendix C for more information). By expressing the 1-norm of the qubit Hamiltonian in Eq. (5) with respect to the electronic integrals, we can compute its value before doing any fermion-to-qubit transformation, which can be costly for large systems.

### C. Localized orbitals

As discussed previously, different approaches have already been considered to minimize the 1-norm of a given Hamiltonian $\hat{H}$. In this paper, we choose to tackle the problem from a chemistry point of view by focusing on the use of orbital transformations (i.e., single-particle states rotations). More precisely, we investigate the use of orbital localization techniques as a classical preoptimization method to express the electronic structure Hamiltonian $\hat{H}$ in a basis (see Appendix B for details about electronic integral transformations) that presents natural advantages for quantum computing. Note that exploiting spatial locality to reduce the 1-norm has already been mentioned in Ref. [52], or to reduce the number of significant integrals in Ref. [54].

In computational chemistry, localization schemes represent state-of-the-art orbital-rotation techniques employed in various situations. For example, localized orbitals (LOs) are regularly used to alleviate the computational cost of numerical simulations in post-HF methods such as second-order Møller Plesset [60–63], coupled cluster [64–67], and multireference methods [68–71]. They can also be used to partition a system in spatially localized subsystems that are treated at different levels of theory [72–84]. In the context of quantum algorithms for the NISQ era, one may, for instance, consider performing a calculation with a classical mean-field method to produce localized orbitals and using the orbitals localized in the spatial region of interest as a basis for a calculation with a quantum algorithm. In the current paper, we demonstrate that LOs can also be of significant help beyond isolating the spatial region that is of chemical interest by reducing the qubit 1-norm $\lambda_Q$ of the electronic Hamiltonian $\hat{H}$ after a fermion-to-qubit transformation.

In the following, we introduce the orbital localization schemes considered in this paper where the notations $\tilde{\psi}_p(\mathbf{r})$ and $\chi_\mu(\mathbf{r})$ are used to denote orthogonal LOs and nonorthogonal AOs, respectively.

#### 1. Lowdin orthogonal atomic orbital

The first approach we investigate to generate LOs is the orthogonalization of atomic orbitals method (OAO). In practice, several techniques exist to produce orthogonal AOs [85,86]. Here, we focus on Löwdin’s method [87,88], known to generate orthogonal LOs with a shape that is the closest to the original AOs (in the least square sense). In practice, orthogonal Löwdin orbitals $\tilde{\psi}_p(\mathbf{r})$ are built via a linear combination of the N original AOs,

$$ \tilde{\psi}_p(\mathbf{r}) = \sum_{\mu=1}^{N} \chi_{\mu}(\mathbf{r}) \tilde{C}_{\mu p}, \quad (23) $$

where the orbital coefficient matrix $\tilde{C}$ takes a very simple form like

$$ \tilde{C} = \mathbf{S}^{-1/2}, \quad (24) $$

and $S_{\mu\nu} = \int d\mathbf{r} \chi_\mu^*(\mathbf{r}) \chi_\nu(\mathbf{r})$ is the overlap matrix encoding the overlap between different nonorthogonal AOs.

From a practical point of view, this method represents one of the simplest and numerically cheapest localization methods, where the computational cost is dominated by the exponentiation of the overlap matrix and typically scales as $O(N^3)$. However, this approach is defined with respect to the full AO space and cannot be straightforwardly applied to the practical case of active space calculations where only a restricted set of MOs—formed by linear combination of AOs—is used.

#### 2. Molecular orbital localization schemes

Fortunately, other localization schemes can be used in any circumstances (i.e., with or without active space approximation). Among the possible approaches, we focus on three of them: the Pipek-Mezey [89] (PM), Foster-Boys [90] (FB), and Edmiston-Ruedenberg [91] (ER) methods. From a practical point of view, these methods fundamentally differ from OAO as they generate LOs out from MOs and not AOs. In practice, all these methods start from a set of orthogonal canonical MOs (CMOs) $\{\phi_s\}_{s=1}^N$ obtained with an initial mean-field calculation (e.g., HF or density functional theory). The CMOs form linear combinations of AOs as

$$ \phi_s(\mathbf{r}) = \sum_{\mu}^{N} \chi_{\mu}(\mathbf{r}) C_{\mu s}, \quad (25) $$

with $\tilde{C}$ the CMO-coefficient matrix. A set of LOs $\{\tilde{\psi}_p\}_{p=1}^N$ is then generated by applying a unitary transformation matrix $\tilde{U}$ (with $\tilde{U}^\dagger\tilde{U} = \tilde{U}\tilde{U}^\dagger = \mathbf{1}$) to express each LO as a linear combination of the original CMOs,

$$ \tilde{\psi}_p(\mathbf{r}) = \sum_{s}^{N} \phi_s(\mathbf{r}) U_{sp}, \quad (26) $$

or, in a more compact matrix form,

$$ \tilde{C} = \mathbf{CU}, \quad (27) $$

where $\tilde{C}$ is the LO-coefficient matrix. In practice, the shape of the LOs is numerically determined by modifying the unitary $\tilde{U}$ to optimize a cost function $\mathcal{L}$ based on a relevant localization criterion (depending on the scheme considered). This localization can be prone to many local minima and the OO process is usually realized with different numerical techniques (e.g., Jacobi rotations [92,93], gradient descent [94], Newton, and trust-region methods [95–97], etc.). In chemistry, the question of how to realize an efficient OO process in practice still constitutes an active topic of research which goes way beyond the scope of our paper (we refer the interested reader to Ref. [98] and references within). Let us now focus on the different criteria used in the FB, PM, and ER methods together with their respective numerical costs.

In the FB scheme [90], the localization criterion is the square of the distance separating two electrons $\mathbf{r}_{12}^2 = |\mathbf{r}_2 - \mathbf{r}_{12}|$.
$\mathbf{r}_1|^2$ and the set of LOs is obtained by minimizing the following cost function

$$ \mathcal{L}_{\text{FB}} = \sum_{p}^{N} \iint |\tilde{\psi}_{p}(\mathbf{r}_{1})|^2 \mathbf{r}_{12}^2 |\tilde{\psi}_{p}(\mathbf{r}_{2})|^2 d\mathbf{r}_{1} d\mathbf{r}_{2}, \quad (28) $$

which represents the average value of the distance criterion over the set of orbitals to be optimized. In practice, a single estimation of the cost function $\mathcal{L}_{\text{FB}}$ requires a decomposition in the AO basis that generates five nested sums, thus leading to a scaling of $O(N^5)$. However, specific manipulations (see Refs. [90] and [98]) can reduce the total computational cost of the FB method to $O(N^3)$ multiplied by the number of times $\mathcal{L}_{\text{FB}}$ is called (which intrinsically depends on the optimization algorithm considered in practice). Note that extensions of the FB scheme exist and are based on the orbital variance [99] or the fourth-moment method [100], resulting in more localized orbitals, especially for basis sets augmented by diffuse functions.

In the case of the ER method, the inverse distance $1/\mathbf{r}_{12}$ [proportional to the two-body electronic repulsion operator, see Eq. (3)] is used as a criterion, and the LOs are obtained by maximizing

$$ \mathcal{L}_{\text{ER}} = \sum_{p}^{N} \iint |\tilde{\psi}_{p}(\mathbf{r}_{1})|^2 \frac{1}{\mathbf{r}_{12}} |\tilde{\psi}_{p}(\mathbf{r}_{2})|^2 d\mathbf{r}_{1} d\mathbf{r}_{2}. \quad (29) $$

The computational cost of a single estimation of $\mathcal{L}_{\text{ER}}$ scales as $O(N^5)$ (multiplied by the number of times $\mathcal{L}_{\text{ER}}$ is called), though it can be reduced up to $O(N^3)$ by the use of density fitting or the Cholesky decomposition.

Finally, the PM scheme adopts a very different point of view. PM orbitals are generated by maximizing the Mulliken charge of each orbital [89],

$$ \mathcal{L}_{\text{PM}} = \sum_{A=1}^{N_A} \sum_{p}^{N} (Q_A^p)^2, \quad (30) $$

where

$$ Q_A^p = \sum_{\mu \in A} \sum_{v}^{N} \tilde{C}_{\mu p} S_{\mu v} \tilde{C}_{vp} \quad (31) $$

is the contribution of orbital $p$ to the Mulliken charge of atom A. In this case, the numerical cost of one call of $\mathcal{L}_{\text{PM}}$ essentially depends on the triple sums over the MOs and AOIs in Eqs. (30) and (31). As a result, the computational cost of the global PM method scales as $O(N^3)$ (multiplied by the number of times $\mathcal{L}_{\text{PM}}$ is called). Note that the traditional PM method is ill-defined since Mulliken charges are basis set sensitive and do not have a basis set limit. Various partial charge estimates can be used instead of Mulliken charges [101], such as intrinsic orbitals that are basis set insensitive and lead to a cheaper and better behaved localization procedure than PM localization [102,103].

In summary, the construction of OAOs is the cheapest approach but can't be used when considering active space. Somewhat more expensive are PM and FB, which share an equivalent scaling, and, finally, ER, which is the most computationally demanding method.

### D. 1-norm orbital-optimization

One of the main contributions of this paper is the use of an OO process specifically dedicated to the minimization of the qubit 1-norm $\lambda_Q$. To proceed, we introduce a unitary operator

$$ \mathbf{U}^{\text{OO}} = e^{-\mathbf{K}} \quad (32) $$

with an anti-Hermitian generator

$$ \mathbf{K}^T = -\mathbf{K} = \begin{pmatrix} 0 & K_{12} & K_{13} & \dots & K_{1N} \\ -K_{12} & 0 & K_{13} & \dots & K_{2N} \\ -K_{13} & -K_{13} & 0 & \dots & K_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -K_{1N} & -K_{2N} & -K_{3N} & \dots & 0 \end{pmatrix}. \quad (33) $$

This operator is then used to transform a reference MO basis into a new basis denoted by $\{\phi'_q\}$ with $\mathbf{C}' = \mathbf{CU}^{\text{OO}}$, for which the optimal orbitals are obtained by minimizing the following cost function:

$$ \mathcal{L}_{\text{OO}} = \lambda_Q(\{\phi'_q\}) \quad (34) $$

by varying the $N(N-1)/2$ off-diagonal independent elements of the matrix $\mathbf{K}$. This optimization can be carried out for a reference MO basis that consists of canonical HF orbitals, but could equally well be carried out on a subspace of the full MO-basis, e.g., only the set of fractionally occupied orbitals in an active space-type calculation, or only a subset of localized orbitals resulting from a preliminary localization using one of the methods discussed above.

In practice, to realize the OO of the 1-norm, we choose to use brute-force optimization algorithms which autonomously estimate the local derivatives of the cost function (i.e., gradients and Hessians). This choice is motivated by the expression of $\lambda_Q$ in Eq. (19) that contains many absolute values and for which analytical derivatives are clearly nontrivial to estimate. Numerically, the main bottleneck of the OO method is linked to the repetitive transformation of the two-electron integrals realized at each step of the process. As a result, the core algorithm scales as $O(N^5)$ (similar to the ER localization scheme). However, with the computational effort deployed to estimate the numerical gradient and Hessian of the complicated cost function [given in Eq. (34)], this scaling should increase more especially when treating large systems.

### III. COMPUTATIONAL DETAILS

The geometries of the small systems considered in this paper were optimized using the ADF program of the Amsterdam Modeling Suite [104] with a quick universal force-field optimization, sufficient for our purposes. These geometries are provided in the Supplemental Material [105]. The first geometry in Ref. [1] with a charge of +3 is used for FeMoco and the geometries in Ref. [4] (all with a charge of +1) were used for the Ruthenium metal complexes. The electronic integrals of the Hamiltonian were computed using the restricted HF method from the PySCF package [106]. For large molecules, the frozen-core approximation was invoked according to Appendix A. All the localization schemes used to transform the Hamiltonian were already implemented in the PySCF package. For our 1-norm OO scheme, the SLSQP optimizer was used from the SciPy package [107] and the Python version
of the L-BFGS-B optimizer of the OPTIMPARALLEL package [108]. The choice of these optimizers has been motivated by their capacity to automatically approximate gradients (and Hessians) for minimization. This represents an interesting tool when no evident analytical gradient/Hessian can be determined as for the 1-norm. A fast algorithm to estimate this $\lambda_Q$ [Eq. (19)] was implemented in the OPENFERMION package [109], allowing users to calculate 1-norms of large molecular Hamiltonians without employing expensive fermion-to-qubit mappings.

## IV. RESULTS

In this section, we study the scaling of $\lambda_Q$ with respect to the number of orbitals. First, we fix the basis set and increase the number of atoms, in the spirit of Refs. [33,52]. Concerning the orbital localization schemes, we allow rotations between the active occupied and virtual spaces as this can lead to a better localized orbital basis.

Then we benchmark the value of $\lambda_Q$ for several different chemical systems and active space sizes, using the different localization schemes and our 1-norm orbital-optimized scheme. When considering an active space, we always choose our active space based on the CMOs by considering the $N_{\text{act}}$ orbitals around the Fermi level (i.e., around the highest occupied and lowest unoccupied MOs). We localize only inside the active space (i.e., the localizing unitary in Eq. (26) only has indices corresponding to active orbitals) such that the subspace spanned by the active space remains invariant under the unitary rotations. Although this means that the expectation values of observables remain the same inside the active space when an exact solver is used, one can converge to different expectation values when approximate solvers are considered, such as the truncated unitary coupled cluster ansatz (that becomes exact if not truncated [110]). However, rotating from the CMO to the LMO basis will not necessarily deteriorate the results of the (approximate) simulation, as LMOs have significant importance in local correlation treatments in post-HF methods like second-order Møller Plesset [60–63], coupled cluster [64–67], embedding approaches [72–84], and multireference methods [68–71] on classical computers. Hence, using localized orbitals could lead to similar advantages in quantum computing simulations, and could also inspire new ansatz proposals based on embedding strategies. This is left for future work.

Finally, to study the scaling of $\lambda_Q$ for large molecules, we increase the size of the active space $N_{\text{act}}$ while fixing the basis set and the number of atoms.

### A. Hydrogen and alkane chains: Scaling of the 1-norm by increasing the number of atoms

In this section, we investigate the performance of the localization schemes together with our brute-force OO method, studying the scaling of $\lambda_Q$ with respect to the number of orbitals by increasing the number of atoms in small systems. Inspired by Refs. [33,54], we consider linear chains of hydrogens with a spacing of $r = 1.4$ Å and linear alkane chains (for which the geometry has been optimized). Results are shown in Fig. 1 and Table II.

<figure><img src="image_4.png" /><figcaption>FIG. 1. The 1-norm $\lambda_Q$ (in Hartree) with respect to (a) the number of atoms in a linear hydrogen chain (same as the number of orbitals) and (b) the number of orbitals in a linear alkane chains (ethane, propane, butane, etc. up to decane), for different orbitals in the STO-3G basis. The constant term in the Hamiltonian is ignored. Assuming a polynomial scaling of $\lambda_Q = O(N^\alpha)$, we fit $\log_{10} \lambda_Q = \alpha \log_{10} N + \beta$ and show the plot with log-log axes. The fitting and regression coefficients are given in Table II.</figcaption></figure>

The orbital-optimizer method is not included for the linear hydrogen chains, because the gradient-estimating optimizer has trouble finding the approximate gradient in the first step of the optimization. The reason for this is that the localization schemes already come so close to the optimal solution that the optimizer estimates the gradient to be zero or positive in each direction. The STO-3G basis results in the localized orbitals resembling just 1s AOs on every nucleus which are very localized and therefore have a minimal 1-norm. As linear hydrogen chains are very artificial systems and give limited information on the scaling of $\lambda_Q$ for actual interesting chemical systems, we thought the results of the conventional localization schemes to be sufficient here.

As readily seen, all localization schemes perform comparably well in reducing the 1-norm compared to canonical MOs. Indeed, the scaling of the 1-norm decreases significantly from $O(N^{2.31})$ to $O(N^{1.34})$ for the hydrogen chains, and from $O(N^{2.21})$ to $O(N^{1.38})$ for the alkane chains. In terms of concrete values, for the largest $N$ studied in Fig. 1, we can see a reduction of $\lambda_Q$ by more than a factor of 13 for the hydrogen
<table><caption>TABLE II. Fitting and regression coefficients from Fig. 1.</caption><thead><tr><th>Hydrogen chain</th><th>α</th><th>β</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>CMO</td><td>2.31</td><td>-1.54</td><td>0.9990</td></tr><tr><td>PM</td><td>1.47</td><td>-0.70</td><td>0.8832</td></tr><tr><td>FB</td><td>1.34</td><td>-0.61</td><td>0.9967</td></tr><tr><td>ER</td><td>1.34</td><td>-0.54</td><td>0.9643</td></tr><tr><td>OAO</td><td>1.46</td><td>-0.99</td><td>0.9846</td></tr><tr><td>Alkane chain</td><td>α</td><td>β</td><td>R<sup>2</sup></td></tr><tr><td>CMO</td><td>2.21</td><td>-0.82</td><td>0.9978</td></tr><tr><td>PM</td><td>1.41</td><td>0.93</td><td>0.9965</td></tr><tr><td>FB</td><td>1.38</td><td>1.08</td><td>0.9976</td></tr><tr><td>ER</td><td>1.39</td><td>1.13</td><td>0.9926</td></tr><tr><td>OAO</td><td>1.43</td><td>0.92</td><td>0.9954</td></tr><tr><td>OO</td><td>1.44</td><td>0.76</td><td>0.9964</td></tr></tbody></table>

and 5 for the alkane chains, when passing from the CMO basis to the LO ones. In the alkane chains, we do not see any significant distinction between the different schemes. While α for the orbital-optimization scheme is slightly higher than the best localization scheme, it is not so informative here as the λ_Q associated to the OO has the lowest absolute value for any orbital space. In the hydrogen chains, it seems that the PM scheme is slightly less efficient. Surprisingly, the AO-based Löwdin orthogonalization method does not provide any further improvement over other localized orbitals, which might be due to the simplicity of the systems for which localization schemes can generate extremely localized orbitals.

## B. Benchmarking $\lambda_Q$ for a variety of molecules and active spaces

To benchmark the impact of the localization schemes on the 1-norm, we calculated $\lambda_Q$ for a variety of molecules and active spaces. Whenever we consider the full space of orbitals, we include the value of $\lambda_Q$ obtained from the Hamiltonian expressed in the OAO basis. When an active space is considered, we use the frozen-core approximation and we localize

the orbitals inside the active space. The resulting $\lambda_Q$ values are tabulated in Table III, where the lowest ones obtained from standard localization schemes for each system are in bold. Apart from the different localization schemes, results of our brute force optimizer are also shown. Note that we start the optimization with already localized orbitals, such that it always gives a lower 1-norm than the best localization scheme. On most molecules, we had the best experience with the L-BFGS-B optimizer. Unfortunately, this optimizer was not able to estimate the gradient of $\lambda_Q$ with respect to $U^{OO}$ on the smallest molecules (H₂ and LiH), which may be due to infinitesimal changes in $\lambda_Q$ when changing $U^{OO}$. This is why we employed the SLSQP optimizer here.

To identify the origin of the significant reduction of the 1-norm and its scaling, we focus on the transformation of the two-electron integrals defined in Eq. (3) with the following permutational symmetries (for real-valued integrals):

$$ \begin{aligned} (pq|rs) &= (pq|sr) = (qp|rs) = (qp|sr) \\ &= (rs|pq) = (rs|qp) = (sr|pq) = (sr|qp), \end{aligned} \quad (35) $$

and we divide the indices of this tensor in seven separate classes: (1) pppp, (2) pqqq, (3) pqpq, (4) ppqq, (5) pqrq, (6) pprs, (7) pqrs, where the indices p, q, r, and s are all different. The above notations are used in Fig. 2 and correspond to the sum of the absolute values of the electronic integrals associated to these indices [and the ones related by symmetries in Eq. (35)], e.g.,

$$ \begin{aligned} \text{pqrq} &\equiv \sum_{p \neq q \neq r} |g_{pqrq}| + |g_{pqqr}| + |g_{qprq}| + |g_{qpqr}| \\ &= \sum_{p \neq q \neq r} 4|g_{pqrq}| \\ \text{pqrs} &\equiv \sum_{p \neq q \neq r \neq s} |g_{pqrs}| \end{aligned} \quad (36) $$

for the cases (5) and (7), respectively. These contributions can give us some insight as to how the localized orbitals affect the norm of the diagonal and off-diagonal parts of the tensor and

<table><caption>TABLE III. Values of $\lambda_Q$ for relatively small test molecules in the cc-pVDZ basis-set. The active space size is indicated as ($n_e, N$), where $n_e$ is the number of electrons in $N$ the number of spatial orbitals. The lowest 1-norm obtained from standard localization schemes for each system are in bold. Superscripts (a) and (b) refer to the use of the SLSQP and the L-BFGS-B optimizer used for in our 1-norm orbital-optimization scheme (denoted as optimizer in the table), respectively. The rightmost column shows the percentage of reduction obtained from the 1-norm-optimized orbitals compared to the CMOs.</caption><thead><tr><th>Molecule</th><th>Active space</th><th>CMO</th><th>PM</th><th>FB</th><th>ER</th><th>OAO</th><th>Optimizer</th><th>% reduction</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>Full space (2,10)</td><td>101</td><td>135</td><td>116</td><td><b>93</b></td><td>103</td><td>90<sup>(a)</sup></td><td>10.9%</td></tr><tr><td>LiH</td><td>Full space (4,19)</td><td>185</td><td>177</td><td>190</td><td>178</td><td><b>153</b></td><td>134<sup>(a)</sup></td><td>27.6%</td></tr><tr><td>H<sub>2</sub>O</td><td>Full space (10,24)</td><td>717</td><td>678</td><td>710</td><td>662</td><td><b>616</b></td><td>576<sup>(b)</sup></td><td>19.7%</td></tr><tr><td>HLiO</td><td>Full space (12,33)</td><td>993</td><td>787</td><td>831</td><td><b>768</b></td><td>788</td><td>668<sup>(b)</sup></td><td>32.7%</td></tr><tr><td>H<sub>2</sub>CO</td><td>Full space (16,38)</td><td>1792</td><td>1419</td><td>1417</td><td>1343</td><td><b>1327</b></td><td>1101<sup>(b)</sup></td><td>38.6%</td></tr><tr><td>HNCH<sub>2</sub></td><td>Full space (16,43)</td><td>3096</td><td>1813</td><td>1670</td><td><b>1588</b></td><td>1645</td><td>1240<sup>(b)</sup></td><td>60.0%</td></tr><tr><td>C<sub>3</sub>H<sub>6</sub></td><td>(24,45)</td><td>2316</td><td>1223</td><td>1145</td><td><b>1137</b></td><td>N/A</td><td>995<sup>(b)</sup></td><td>57.0%</td></tr><tr><td>C<sub>4</sub>H<sub>6</sub></td><td>(30,45)</td><td>1760</td><td><b>974</b></td><td>1054</td><td>1048</td><td>N/A</td><td>812<sup>(b)</sup></td><td>53.9%</td></tr><tr><td>C<sub>5</sub>H<sub>8</sub></td><td>(38,45)</td><td>1821</td><td>818</td><td>801</td><td><b>799</b></td><td>N/A</td><td>698<sup>(b)</sup></td><td>61.7%</td></tr><tr><td>HNC<sub>3</sub>H<sub>6</sub></td><td>(32,50)</td><td>3796</td><td>1493</td><td>1244</td><td><b>1232</b></td><td>N/A</td><td>1085<sup>(b)</sup></td><td>71.4%</td></tr><tr><td>HNC<sub>5</sub>H<sub>10</sub></td><td>(48,50)</td><td>2632</td><td>1098</td><td><b>989</b></td><td>1002</td><td>N/A</td><td>842<sup>(b)</sup></td><td>68.0%</td></tr><tr><td>HNC<sub>7</sub>H<sub>14</sub></td><td>(50,50)</td><td>2610</td><td>751</td><td><b>705</b></td><td>707</td><td>N/A</td><td>616<sup>(b)</sup></td><td>76.4%</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Sum of absolute values of different parts of the two-electron integral tensor on a logarithmic scale (see text for further details on the notations). It shows the second stable Ruthenium intermediate in a minimal ANO-RCC basis for an active space of 100 electrons in 100 spatial orbitals.</figcaption></figure>

are represented in Fig. 2 for the second stable intermediate Ruthenium complex in different orbital bases. As readily seen in this diagram, localization schemes do not lead to a uniform reduction of the norm of every term, but tend to maximize some of them (the pppp term), minimize others (the pqrs term), and leave others relatively unchanged (the ppqq term). However, the pqrs contribution clearly dominates in the CMO basis, being more than one order of magnitude larger than any of the other terms. Especially, the orbital-optimizer can give an even better reduction of this term. This explains why the localization schemes can be used to reduce the 1-norm, as reducing the pqrs norm by one order of magnitude (as seen in Fig. 2) will play a much more important role than increasing the pppp one by the same order of magnitude. As the magnitude of these integrals depends on the overlap densities $\psi_p(\mathbf{r})\psi_q(\mathbf{r})$ and $\psi_r(\mathbf{r})\psi_s(\mathbf{r})$, it is clear that localization will reduce the number of numerically significant contributors. For an extended system, integrals in which p and q as well as r and s are localized on the same atoms are expected to dominate. This explains the observation that the pqrs norm becomes comparable to the ppqq and pprs terms when employing localized orbitals.

### C. Effect of increasing the size of the active space on $\lambda_Q$

Let us now investigate how $\lambda_Q$ scales with respect to the number of active orbitals for large and complex systems. We selected test molecules based on two important recent papers in the field of quantum computing for quantum chemistry: FeMoco—an iron-molybdenum complex that constitutes the active site of a MoFe protein which plays an important role in nitrogen fixation [1]—and all the stable intermediates and transition states of the catalytic cycle presented in Ref. [4]. This catalytic cycle describes the binding and transformations of carbon dioxide $CO_2$, a molecule with infrared absorption properties that makes it a potent and the most important greenhouse gas.

For the calculations to remain doable in a reasonable amount of time, the ANO-RCC minimal basis (with scalar relativistic corrections) is used. This should be sufficient for the current purpose of studying 1-norm reductions within an active space, but we note that larger basis sets will be required to reach chemical accuracy [111]. The scaling of $\lambda_Q$ with

<figure><img src="image_6.png" /><figcaption>FIG. 3. An illustration of how to decrease the 1-norm on the formaldimine molecule: Moving from CMOs to LOs with a Jacobi rotation between occupied $\pi$-bonding and virtual $\pi^*$-antibonding orbitals.</figcaption></figure>

respect to the number of active orbitals $N_{act}$ for different orbital bases are reported in detail in Fig. 4. Comparing the scaling reported in this section and the ones obtained in Sec. IV A, we see that the scaling is larger when increasing the number of active orbitals than when increasing the number of atoms, as expected. However, the impact of localized orbitals on the 1-norm is similar in both cases, as it reduces the scaling by almost one order of magnitude. We also do not see any significant distinction between the 1-norm associated to the different conventional localization schemes. The brute-force orbital-optimizer, however, gives a consistent 10–25% improvement over the conventional localization schemes. In the largest active spaces considered ($N_{act} \ge 85$), the orbital-optimization scheme can converge very slowly due to the need for four-index transformations at every step, combined with the noncontinuous and complicated landscape of $\lambda_Q$. Because of this, the user is usually forced to cut off the optimizer at some point or choose a high threshold of convergence. In spite of this, we got similar good results in the largest active spaces we considered, indicating that the optimizer is not limited to a specific active space size.

The scaling results for all large molecules are summarized in Table IV. The scaling lies in the order of $O(N^{2.6})$ to $O(N^{2.9})$ when using CMOs, while it lies around $O(N^2)$ or slightly higher when using LOs. One gets a bit of a misleading picture looking at just the scaling of the orbital optimizer, as it has the lowest absolute value of $\lambda_Q$ at every point. The cause is that the cost function of the orbital optimizer has more local minima and is harder to converge compared to conventional localization schemes.

For the localization schemes, it is usual that the bigger the active space is, the more $\lambda_Q$ can be reduced compared to CMOs. This can be rationalized as follows. When the active space increases, it may involve orbitals localized on atoms which have negligible overlap with the other orbitals. Another possibility is that the additional orbitals may have a positive effect on the localization procedure. As an illustration, see the example shown in Fig. 3 where we consider an occupied $\pi$-bonding orbital mixing with its associated virtual $\pi^*$-antibonding orbital for the formaldimine molecule. By rotating those two orbitals via a Jacobi rotation, we move from a
<table><caption>TABLE IV. Scaling of $\lambda_Q$ with respect to the active space size $N_{\text{act}}$ for different orbital bases, assuming $\log_{10} \lambda_Q = \alpha \log_{10} N_{\text{act}} + \beta$ with the associated $R^2$ value. The minimal ANO-RCC basis-set is considered.</caption><thead><tr><th>Molecule</th><th>α<sub>CMO</sub></th><th>β<sub>CMO</sub></th><th>R<sup>2</sup><sub>CMO</sub></th><th>α<sub>PM</sub></th><th>β<sub>CMO</sub></th><th>R<sup>2</sup><sub>PM</sub></th><th>α<sub>FB</sub></th><th>β<sub>CMO</sub></th><th>R<sup>2</sup><sub>FB</sub></th><th>α<sub>ER</sub></th><th>β<sub>CMO</sub></th><th>R<sup>2</sup><sub>ER</sub></th><th>α<sub>OO</sub></th><th>β<sub>OO</sub></th><th>R<sup>2</sup><sub>OO</sub></th></tr></thead><tbody><tr><td>FeMoco</td><td>2.83</td><td>-3.37</td><td>0.9996</td><td>2.03</td><td>-1.60</td><td>0.9477</td><td>1.94</td><td>-1.26</td><td>0.9580</td><td>1.98</td><td>-1.49</td><td>0.949</td><td><b>1.89</b></td><td>-1.2</td><td>0.9688</td></tr><tr><td>Ru SI I</td><td>2.63</td><td>-2.84</td><td>0.9957</td><td>2.36</td><td>-3.13</td><td>0.9465</td><td><b>2.17</b></td><td>-2.47</td><td>0.9865</td><td>2.25</td><td>-2.78</td><td>0.9763</td><td>2.21</td><td>-2.80</td><td>0.9927</td></tr><tr><td>Ru SI II</td><td>2.83</td><td>-3.71</td><td>0.9997</td><td>2.27</td><td>-2.81</td><td>0.9530</td><td><b>2.11</b></td><td>-2.26</td><td>0.9867</td><td>2.22</td><td>-2.68</td><td>0.9770</td><td>2.20</td><td>-2.8</td><td>0.9915</td></tr><tr><td>Ru TS II–III</td><td>2.95</td><td>-4.22</td><td>0.9998</td><td>2.16</td><td>-2.29</td><td>0.8919</td><td><b>1.98</b></td><td>-1.66</td><td>0.9659</td><td>2.13</td><td>-2.26</td><td>0.9234</td><td>2.12</td><td>-2.44</td><td>0.9946</td></tr><tr><td>Ru SI V</td><td>2.75</td><td>-3.41</td><td>0.9979</td><td>2.42</td><td>-3.37</td><td>0.9063</td><td><b>2.13</b></td><td>-2.33</td><td>0.9798</td><td>2.33</td><td>-3.09</td><td>0.9526</td><td>2.21</td><td>-2.83</td><td>0.9926</td></tr><tr><td>Ru SI VIII</td><td>2.84</td><td>-3.79</td><td>0.9989</td><td>2.22</td><td>-2.54</td><td>0.8630</td><td><b>2.12</b></td><td>-2.24</td><td>0.9659</td><td>2.18</td><td>-2.48</td><td>0.9323</td><td>2.19</td><td>-2.73</td><td>0.9901</td></tr><tr><td>Ru TS VIII–IX</td><td>2.82</td><td>-3.57</td><td>0.9995</td><td>2.23</td><td>-2.51</td><td>0.8645</td><td><b>2.09</b></td><td>-2.07</td><td>0.9706</td><td>2.19</td><td>-2.44</td><td>0.9039</td><td>2.19</td><td>-2.66</td><td>0.9883</td></tr><tr><td>Ru SI IX</td><td>2.8</td><td>-3.66</td><td>0.9987</td><td>2.42</td><td>-3.36</td><td>0.9240</td><td><b>2.22</b></td><td>-2.69</td><td>0.9876</td><td>2.34</td><td>-3.14</td><td>0.9663</td><td>2.23</td><td>-2.91</td><td>0.9967</td></tr><tr><td>Ru SI XVIII</td><td>2.79</td><td>-3.46</td><td>0.9986</td><td>2.31</td><td>-2.87</td><td>0.9302</td><td><b>2.12</b></td><td>-2.16</td><td>0.9854</td><td>2.26</td><td>-2.74</td><td>0.9557</td><td>2.24</td><td>-2.87</td><td>0.9929</td></tr></tbody></table>

maximal value of $\lambda_Q$ for $\theta = 0$ (corresponding to the original CMOs) to a minimal value for $\theta = \pi/4$ (corresponding to the new LOs). This way, we see that the larger the active space, the more degrees of freedom we have to perform a better orbital localization. For the biggest active space chosen, an

active space of 100 electrons in 100 orbitals, $\lambda_Q$ is reduced by an order of magnitude. This illustrates the potential of using widely available and computationally cheap classical preoptimizations, together with our dedicated orbital-optimizer, to reduce the cost of quantum simulations.

<figure><img src="image_7.png" /><figcaption>FIG. 4. Scaling of $\lambda_Q$ (in Hartree) with respect to the size of the active space for different orbitals for (a) FeMoco and (b) the Ruthenium transition state II-III in a minimal ANO-RCC basis. The constant term in the Hamiltonian is ignored. Assuming polynomial scaling of $\lambda_Q = O(N_{\text{act}}^\alpha)$, we fit $\log_{10} \lambda_Q = \alpha \log_{10} N_{\text{act}} + \beta$ and show the plot with log-log axes. The fitting parameters and regression coefficients are given in Table IV.</figcaption></figure>

V. CONCLUSIONS

Numerically, we have seen that exploiting the locality of the basis gives rise to a lower variance in the Hamiltonian coefficients, reducing $\lambda_Q$. This results in a significant reduction of the absolute values of the integrals. The off-diagonal elements of the two-electron tensor play the biggest role here.

To find an even better basis rotation that minimizes $\lambda_Q$, one can use a brute-force orbital optimizer as we have done in this paper. This method could be realized in practice only because we have an efficient way of computing $\lambda_Q$ directly in terms of the molecular integrals [Eq. (19)], avoiding doing a fermion-to-qubit mapping explicitly. Expressing it as $\lambda_Q = \lambda_T + \lambda_V$, this improves a bit over the equation of $\lambda_V$ in Ref. [52], holding into account the representation of the Hamiltonian in terms of unique Pauli strings. The direct orbital optimizer can be used in quite large active spaces, reliably converging for active spaces with up to $N_{\text{act}} \approx 80$ orbitals, sometimes more. For the largest active spaces considered, it converges very slowly due to the need for four-index transformations at every step, which are costly because of the large dimensionality of the system, combined with the noncontinuous and complicated landscape of $\lambda_Q$. This gives rise to the need to either set a high threshold of convergence or cut off the optimizer at some point, where the user is satisfied with the result. While this is problematic, there is no a priori reason not to try this optimizer on the (potentially large) active spaces of the molecule one is considering. This can be a point of further study.

We benchmarked a range of various molecules and active spaces for which we showed achieving a significant reduction of $\lambda_Q$. Apart from this benchmark, we investigated the scaling of $\lambda_Q$ with the size of the system $N$. There are multiple paths one can take here to increase $N$. A popular way to do this in the literature is considering either a hydrogen chain or hydrogen ring and increase the amount of atoms. We benchmarked the scaling of $\lambda_Q$ for a hydrogen chain and a linear alkane chain by increasing the number of atoms as well, where our results show a scaling of $O(N^{2.3})$ and $O(N^{2.2})$, respectively. Here we
saw localized orbitals can give approximately a factor of N lower scaling of $O(N^{1.3})$ and $O(N^{1.4})$, respectively. As this gives one limited information how big $\lambda_Q$ will be in real eventual applications of quantum algorithms, we decided to investigate the scaling on relevant highly correlated molecules such as FeMoco (important in nitrogen fixation) and ruthenium metal complexes (important in carbon dioxide capture). Even though we used a minimal basis to make the scaling calculations feasible on these molecules, there is no reason to believe $\lambda_Q$ will scale differently with respect to the active space size on a bigger basis set. Here we showed that the scaling is a factor of N larger, and with localized orbitals can be brought down to $\lambda_Q \approx O(N^2)$, depending on the molecule in consideration. Our dedicated orbital optimizer was able to consistently give an even further improvement of 10–25% on these molecules. As a concrete example, consider the largest active space considered of FeMoco: this factor of N difference would result in a reduction by two orders of magnitude in the amount of measurements needed for tomography [Eq. (14)], if one wants chemical accuracy. This indicates that the simple and efficient classical preprocessing by widely available localization techniques, together with our dedicated optimizer, will help to make simulations with large active spaces feasible.

## ACKNOWLEDGMENTS

We sincerely thank Susi Lethola for his helpful comments on localization schemes. S.Y. and B.S. acknowledge support from the Dutch Research Council (NWO). E.K. acknowledges support from Shell Global Solutions BV.

## APPENDIX A: FROZEN CORE HAMILTONIAN

Applying the frozen core approximation to the electronic-structure Hamiltonian consists of assuming the existence of a set of frozen orbitals (always occupied), another set of active orbitals (belonging to an active space), and a set of virtual orbitals (always unoccupied). Based on this partitioning, every Slater determinant $|\Phi\rangle$ used to describe properties of the system take the form

$$|\Phi\rangle = |\Phi_{\text{frozen}}\Phi_{\text{active}}\rangle, \qquad (A1)$$

where the left contribution $\Phi_{\text{frozen}}$ is the part of the determinant encoding the frozen orbitals of the system (always occupied) whereas $\Phi_{\text{active}}$ is a part encoding the occupancy of the remaining electrons in the active orbitals of the system. In this context, if one considers that every correlated electronic wave function is always expanded with Slater determinants following Eq. (A1), one can demonstrate by projections that the system Hamiltonian takes an effective form

$$\langle\Phi|\hat{H}|\Phi\rangle \equiv \langle\Phi_{\text{active}}|\hat{H}^{\text{FC}}|\Phi_{\text{active}}\rangle, \qquad (A2)$$

with $\hat{H}^{\text{FC}}$ the so-called frozen core Hamiltonian defined as follows:

$$\hat{H}^{\text{FC}} = \hat{H}_{\text{active}} + E_{\text{frozen}}^{\text{MF}} + \hat{V}. \qquad (A3)$$

Here, $\hat{H}_{\text{active}}$ is the Hamiltonian encoding the one- and two-body terms only acting in the active space,

$$\hat{H}_{\text{active}} = \sum_{tu}^{\text{active}} h_{tu} \hat{E}_{tu} + \sum_{tuvw}^{\text{active}} g_{tuvw} \hat{e}_{tuvw}, \qquad (A4)$$

where $t, u, v, w$ denote active space orbitals. The second term $E_{\text{frozen}}^{\text{MF}}$ is a scalar representing the mean-field-like energy obtained from the frozen orbitals,

$$E_{\text{frozen}}^{\text{MF}} = 2 \sum_{i}^{\text{frozen}} h_{ii} + \sum_{ij}^{\text{frozen}} (2g_{iijj} - g_{ijji}), \qquad (A5)$$

and the third term

$$\hat{V} = \sum_{tu}^{\text{active}} V_{tu} \hat{E}_{tu}, \quad \text{with} \quad V_{tu} = \sum_{i}^{\text{frozen}} (2g_{tuii} - g_{tiiu}) \qquad (A6)$$

represents an effective one body potential which encodes the interaction of the frozen electrons with the active space electrons. To summarize, the main effect of the frozen core approximation [Eq. (A3)] is first to introduce an energetic shift [Eq. (A5)] and second to augment the one body term of the Hamiltonian operator [Eq. (A4)] (that only lives in the active space) with an additional effective one body operator [Eq. (A6)].

## APPENDIX B: ELECTRONIC INTEGRALS
TRANSFORMATION

Preparing the electronic Hamiltonian in a given orthogonal orbital basis $\{\phi_p(\mathbf{r})\}_{p=1}^N$ is a crucial step for realizing concrete quantum computing applications. In practice, such a process can be realized classically via electronic integral transformations. For this, we assume that the orthogonal orbitals can be expressed as a linear combination of AOs such that

$$\phi_p(\mathbf{r}) = \sum_{\mu} \chi_{\mu}(\mathbf{r}) C_{\mu p}, \qquad (B1)$$

where we assume that the functions $\chi_{\mu}(\mathbf{r})$ as well as the coefficient matrix $C$ are real valued, as is usually done in nonrelativistic quantum chemistry. Based on our knowledge of this matrix $C$, one can transform the one- and two-electron integrals from the AO basis to the orthogonal orbital basis. To do so, one starts from the one- and two-electron integrals in the AO basis, respectively $h_{\mu\nu}$ and $g_{\mu\nu\delta}$, and implements the following two- and four-index transformations:

$$h_{pq} = \sum_{\mu\nu} h_{\mu\nu} C_{\mu p} C_{\nu q} \qquad (B2)$$

and

$$g_{pqrs} = \sum_{\mu\nu\gamma\delta} g_{\mu\nu\gamma\delta} C_{\mu p} C_{\nu q} C_{\gamma r} C_{\delta s}. \qquad (B3)$$

From a computational point of view, the numerical cost of processing electronic integrals is essentially governed by the four-index transformation. This transformation is known to scale as $O(N^5)$ (see Ref. [112]), when not employing approximations such as density fitting (see Ref. [113]).
APPENDIX C: EXPLICIT FORM OF $\lambda_Q$ IN TERMS OF MOLECULAR INTEGRALS

To derive a formula for $\lambda_Q$ in terms of the molecular integrals, one needs to keep track of what happens to the coefficients in the Hamiltonian [Eq. (A2)] after a fermion-to-qubit transformation of the fermionic operators. This is nontrivial, since, e.g., the Jordan-Wigner transformation transforms a single arbitrary product of fermionic operators $\hat{a}_p^\dagger \hat{a}_s^\dagger \hat{a}_q$ into 16 different Pauli strings because of the products of $\frac{\hat{x}-i\hat{z}}{2}$ and $\frac{\hat{x}+i\hat{z}}{2}$. The Hamiltonian in the form of Eq. (A2) is actually not best suited for this derivation. Instead, it is helpful to express the Hamiltonian in terms of Majorana operators that have the convenient property of being Hermitian operators that square to identity:

$$ \hat{\gamma}_{p\sigma,0} = \hat{a}_{p\sigma} + \hat{a}_{p\sigma}^{\dagger}, \quad \hat{\gamma}_{p\sigma,1} = -i(\hat{a}_{p\sigma} - \hat{a}_{p\sigma}^{\dagger}), \quad (C1) $$

$$ \{\hat{\gamma}_i, \hat{\gamma}_j\} = 2\delta_{ij}\mathcal{I}, \quad \hat{\gamma}_i^\dagger = \hat{\gamma}_i, \quad \hat{\gamma}_i^2 = \mathcal{I}. \quad (C2) $$

To find a representation of the Hamiltonian in terms of Majorana operators, one could directly replace $\hat{a}_p^\dagger = \frac{\hat{\gamma}_{p,0}-i\hat{\gamma}_{p,1}}{2}$ and $\hat{a}_q = \frac{\hat{\gamma}_{q,0}+i\hat{\gamma}_{q,1}}{2}$ in the Hamiltonian or follow the procedure in Ref. [4]. We chose the latter approach, detailed below.

It is well known that the electron repulsion integral tensor, when written as a $N^2 \times N^2$ matrix $g_{pq),(rs)}$ with the composite indices $pq$ and $rs$, is positive semidefinite. This makes it possible to define a Cholesky decomposition $\mathbf{g} = \mathbf{LL}^\dagger$. We write

$$ g_{pqrs} = \sum_{\ell} \sum_{pqrs} L_{pq}^{\ell} L_{rs}^{\ell}, \quad (C3) $$

with $\mathbf{L}$ a lower triangular matrix. Since $g_{pqrs}$ is symmetric in $p, q$ and in $r, s$, for a given $\ell$, $L_{pq}^\ell$ is an $N \times N$ symmetric matrix (also called a Cholesky vector in the Cholesky decomposition). This leads to the following expression for the Hamiltonian:

$$ \begin{aligned} \hat{\mathcal{H}} &= \sum_{pq} \sum_{\sigma} h_{pq} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} + \frac{1}{2} \sum_{pqrs} \sum_{\sigma\tau} g_{pqrs} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{r\tau}^{\dagger} \hat{a}_{s\tau} \hat{a}_{q\sigma} \\ &= \sum_{pq\sigma} \left[ h_{pq} - \frac{1}{2} \sum_{r} g_{prrq} \right] \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} \\ &\quad + \frac{1}{2} \sum_{\ell} \sum_{pqrs} \sum_{\sigma\tau} L_{pq}^{\ell} L_{rs}^{\ell} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} \hat{a}_{r\tau}^{\dagger} \hat{a}_{s\tau} \\ &= \sum_{pq\sigma} h'_{pq} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} + \frac{1}{2} \sum_{\ell} \left( \sum_{pq\sigma} L_{pq}^{\ell} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} \right)^2, \end{aligned} \quad (C4) $$

where $h'_{pq} = h_{pq} - \frac{1}{2} \sum_r g_{prrq}$. Using the relation

$$ \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} + \hat{a}_{q\sigma}^{\dagger} \hat{a}_{p\sigma} = \begin{cases} \mathcal{I} + i(\hat{\gamma}_{p\sigma,0} \hat{\gamma}_{p\sigma,1}), & p=q \\ \frac{i}{2}(\hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} + \hat{\gamma}_{q\sigma,0} \hat{\gamma}_{p\sigma,1}), & p \neq q \end{cases} \quad (C5) $$

one can show that

$$ \sum_{pq\sigma} M_{pq} \hat{a}_{p\sigma}^{\dagger} \hat{a}_{q\sigma} = \sum_{p} M_{pp} \mathcal{I} + \frac{i}{2} \sum_{pq\sigma} M_{pq} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1}, \quad (C6) $$

where $M_{pq}$ can be replaced by any symmetric matrix (like $L_{pq}^\ell$ or $h'_{pq}$). Employing Eq. (C6) in Eq. (C4) leads to the expression of the Hamiltonian in terms of Majorana operators:

$$ \hat{\mathcal{H}} = \sum_{p} h'_{pp} \mathcal{I} + \frac{i}{2} \sum_{pq\sigma} h'_{pq} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} + \frac{1}{2} \sum_{\ell} \left( \sum_{p} L_{pp}^{\ell} \mathcal{I} + \frac{i}{2} \sum_{pq\sigma} L_{pq}^{\ell} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \right)^2. \quad (C7) $$

After working out the square and substituting back $\sum_{\ell} \sum_{pqrs} L_{pq}^{\ell} L_{rs}^{\ell} = g_{pqrs}$, we obtain

$$ \begin{aligned} \hat{\mathcal{H}} = & \left( \sum_{p} h_{pp} + \frac{1}{2} \sum_{pr} g_{pprr} - \frac{1}{2} \sum_{pr} g_{prrp} \right) \mathcal{I} \\ & + \frac{i}{2} \sum_{pq\sigma} \left( h_{pq} + \sum_{r} g_{pqrr} - \frac{1}{2} \sum_{r} g_{prrq} \right) \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \\ & - \frac{1}{8} \sum_{pqrs\sigma\tau} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{r\tau,0} \hat{\gamma}_{s\tau,1}. \end{aligned} \quad (C8) $$

We now want to determine the value of $\lambda_Q$ after a fermion-to-qubit mapping of this Hamiltonian in which all products of Majorana operators should be distinct. We use a fermion-to-qubit mapping such that single Majorana operators are mapped to single unique Pauli strings, such as the Jordan-Wigner transformation [51],

$$ \hat{\gamma}_{i,0} \rightarrow \hat{X}_i \hat{Z}_{i-1} \dots \hat{Z}_0, \quad \text{and} \quad \hat{\gamma}_{i,1} \rightarrow \hat{Z}_i \hat{Z}_{i-1} \dots \hat{Z}_0, \quad (C9) $$

where we used the composite index $i = p\sigma$. Note this is not only true for the Jordan-Wigner transformation but also for example the Bravyi-Kitaev transformation [114]. The only thing left to do is to make sure that the sum over four indices used to evaluate $\lambda_Q$ is indeed truly quartic. Looking at Eq. (C8), this is not yet the case. If $p\sigma = r\tau$ or $q\sigma = s\tau$, the corresponding Majorana operators will square to identity and will reduce to quadratic and identity terms. Below we employ the permutational symmetry of the real-valued integrals $g_{pqrs}$ to obtain the most compact representation.

We first reorder the Majorana's and distinguish between the same and opposite spin cases:

$$ \begin{aligned} & -\frac{1}{8} \sum_{pqrs\sigma\tau} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{r\tau,0} \hat{\gamma}_{s\tau,1} \\ &= \frac{1}{8} \sum_{pqrs} \sum_{\sigma\tau} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\tau,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\tau,1} \\ &= \frac{1}{8} \sum_{pqrs} \sum_{\sigma} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} \\ &\quad + \frac{1}{8} \sum_{pqrs} \sum_{\sigma\neq\tau} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\tau,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\tau,1}. \end{aligned} \quad (C10) $$

For the opposite spin case, all operator products are unique and quartic, so no further work is needed. For the same spin case, it is useful to identify cases in which two or more indices
of the same type of Majorana operators are equal:

$$ 
\begin{aligned}
\frac{1}{8} \sum_{pqrs} \sum_{\sigma} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} &= \frac{1}{8} \sum_{pq} \sum_{\sigma} g_{pqpq} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{q\sigma,1} + \frac{1}{8} \sum_{p \neq r, q} \sum_{\sigma} g_{pqrq} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{q\sigma,1} \\
&\quad + \frac{1}{8} \sum_{q \neq s, p} \sum_{\sigma} g_{pqqs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} + \frac{1}{8} \sum_{p \neq r, q \neq s} \sum_{\sigma} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} \\
&= \frac{1}{4} \sum_{pq} g_{pqpq} \mathcal{I} + \frac{1}{8} \sum_{p \neq r, q \neq s} \sum_{\sigma} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1},
\end{aligned}
\quad (C11) $$

where we applied the identity relation $\hat{\gamma}_i^2 = \mathcal{I}$ and for the second term used that integrals are symmetric in exchanging $p$ and $r$ while the product of Majorana operators is antisymmetric under this exchange. This nullifies this and the third term. The first term can be absorbed in the scalar term. For the truly quartic term, we may make further use of permutational symmetries to get the most compact form:

$$ 
\begin{aligned}
\frac{1}{8} \sum_{p \neq r, q \neq s} \sum_{\sigma} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} \\
&= \frac{1}{8} \sum_{p>r,s>q} \sum_{\sigma} (g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} \\
&\quad + g_{rqps} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} + g_{psrq} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{s\sigma,1} \hat{\gamma}_{q\sigma,1} \\
&\quad + g_{rspq} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{s\sigma,1} \hat{\gamma}_{q\sigma,1}) \\
&= \frac{1}{4} \sum_{p>r,s>q} \sum_{\sigma} (g_{pqrs} - g_{psrq}) \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1}. \quad (C12)
\end{aligned}
$$

Combining all terms, we end up with the expression

$$ 
\begin{aligned}
\hat{H} = & \left( \sum_{p}^{N} h_{pp} + \frac{1}{2} \sum_{pr}^{N} g_{pprr} - \frac{1}{4} \sum_{pr}^{N} g_{prrp} \right) \mathcal{I} \\
& + \frac{i}{2} \sum_{pq\sigma}^{N} \left( h_{pq} + \sum_{r}^{N} g_{pqrr} - \frac{1}{2} \sum_{r}^{N} g_{prrr} \right) \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1} \\
& + \frac{1}{4} \sum_{p>r,s>q} \sum_{\sigma} (g_{pqrs} - g_{psrq}) \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\sigma,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\sigma,1} \\
& + \frac{1}{8} \sum_{pqrs}^{N} \sum_{\sigma \neq \tau} g_{pqrs} \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{r\tau,0} \hat{\gamma}_{q\sigma,1} \hat{\gamma}_{s\tau,1}.
\end{aligned}
\quad (C13) $$

We then take the sum of absolute values of the coefficients and perform the sums over spin in the quartic term explicitly (amounting to a factor of 2) to get the following form of $\lambda_Q$:

$$ \lambda_Q = \lambda_C + \lambda_T + \lambda'_V, \quad (C14) $$

where $\lambda_C$ corresponds to the constant term in Eq. (C13), $\lambda_T$ to the quadratic, and $\lambda_V$ to the quartic term in Majorana operators. They have the form

$$ \lambda_C = \left| \sum_{p}^{N} h_{pp} + \frac{1}{2} \sum_{pr}^{N} g_{pprr} - \frac{1}{4} \sum_{pr}^{N} g_{prrp} \right|, \quad (C15) $$

$$ \lambda_T = \sum_{pq}^{N} \left| h_{pq} + \sum_{r}^{N} g_{pqrr} - \frac{1}{2} \sum_{r}^{N} g_{prrr} \right|, \quad (C16) $$

$$ \lambda'_V = \frac{1}{2} \sum_{p>r,s>q} |g_{pqrs} - g_{psrq}| + \frac{1}{4} \sum_{pqrs} |g_{pqrs}|. \quad (C17) $$

This $\lambda'_V$ is slightly different to the $\lambda_V$ of Ref. [52] in calculating the $\lambda$ for the sparse algorithm of Berry et al. [50], which should actually be the same as $\lambda_Q$. This is due that we took into account the swapping of majorana operators in Eq. (C12). In Ref. [52], this corresponds to the term $\hat{V}' = \frac{1}{8} \sum_{\alpha, \beta \in \{\uparrow, \downarrow\}} \sum_{p,q,r,s} V_{pqrs} \hat{Q}_{p\alpha} \hat{Q}_{r\beta}$. It is hard to see here that the product $\hat{Q}_{p\alpha} \hat{Q}_{r\beta}$ is antisymmetric in swapping $p, r$ and $q, s$, but it becomes clear when one realizes that $\hat{Q}_{p\alpha} = i \hat{\gamma}_{p\sigma,0} \hat{\gamma}_{q\sigma,1}$, indicating the usefulness of Majorana operators. As the absolute values give that

$$ |g_{pqrs} - g_{psrq}| \le |g_{pqrs}| + |g_{psrq}| \quad \forall p, q, r, s \quad (C18) $$

such that $\lambda'_V < \lambda_V$ (except when $g_{pqrs}$ and $g_{psrq}$ always have opposite signs, in which case they would be equal).

[1] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Elucidating reaction mechanisms on quantum computers, Proc. Natl. Acad. Sci. USA 114, 7555 (2017).

[2] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Quantum chemistry in the age of quantum computing, Chem. Rev. 119, 10856 (2019).

[3] B. Bauer, S. Bravyi, M. Motta, and G. Kin-Lic Chan, Quantum algorithms for quantum chemistry and quantum materials science, Chem. Rev. 120, 12685 (2020).

[4] V. von Burg, G. H. Low, T. Häner, D. S. Steiger, M. Reiher, M. Roetteler, and M. Troyer, Quantum computing enhanced computational catalysis, Phys. Rev. Research 3, 033055 (2021).

[5] A. Yu. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, arXiv:quant-ph/9511026.

[6] D. S. Abrams and S. Lloyd, Simulation of Many-Body Fermi Systems on a Universal Quantum Computer, Phys. Rev. Lett. 79, 2586 (1997).

[7] D. S. Abrams and S. Lloyd, Quantum Algorithm Providing Exponential Speed Increase for Finding Eigenvalues and Eigenvectors, Phys. Rev. Lett. 83, 5162 (1999).
[8] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-
Gordon, Simulated quantum computation of molecular ener-
gies, *Science* **309**, 1704 (2005).

[9] T. E. O'Brien, B. Tarasinski, and B. M. Terhal, Quantum phase estimation of multiple eigenvalues for small-scale (noisy) experiments, *New J. Phys.* **21**, 023022 (2019).

[10] J. Preskill, *Quantum computing in the NISQ era and beyond*, *Quantum* **2**, 79 (2018).

[11] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'Brien, A variational eigenvalue solver on a photonic quantum processor, *Nat. Commun.* **5**, 4213 (2014).

[12] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, *New J. Phys.* **18**, 023023 (2016).

[13] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, arXiv:2012.09265.

[14] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, *J. Phys. Soc. Jpn.* **90**, 032001 (2021).

[15] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Noisy intermediate-scale quantum (NISQ) algorithms, arXiv:2101.08448.

[16] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. de Jong, Hybrid quantum-classical hierarchy for mitigation of decoherence and determination of excited states, *Phys. Rev. A* **95**, 042308 (2017).

[17] J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Generalized unitary coupled cluster wave functions for quantum computation, *J. Chem. Theory Comput.* **15**, 311 (2019).

[18] P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. Gambetta, and I. Tavernelli, Quantum equation of motion for computing molecular excitation energies on a noisy quantum processor, *Phys. Rev. Research* **2**, 043140 (2020).

[19] K. M. Nakanishi, K. Mitarai, and K. Fujii, Subspace-search variational quantum eigensolver for excited states, *Phys. Rev. Research* **1**, 033062 (2019).

[20] R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martinez, Hybrid quantum/classical derivative theory: Analytical gradients and excited-state dynamics for the multistate contracted variational quantum eigensolver, arXiv:1906.08728.

[21] R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, Quantum Computation of Electronic Transitions using a Variational Quantum Eigensolver, *Phys. Rev. Lett.* **122**, 230401 (2019).

[22] Y. Ibe, Y. O. Nakagawa, T. Yamamoto, K. Mitarai, Q. Gao, and T. Kobayashi, Calculating transition amplitudes by variational quantum eigensolvers, arXiv:2002.11724.

[23] O. Higgott, D. Wang, and S. Brierley, Variational quantum computation of excited states, *Quantum* **3**, 156 (2019).

[24] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, Variational quantum algorithms for discovering hamiltonian spectra, *Phys. Rev. A* **99**, 062304 (2019).

[25] P. Jouzdani, S. Bringuier, and M. Kostuk, A method of determining excited-states for quantum computation, arXiv:1908.05238.

[26] S. Yalouz, B. Senjean, J. Günther, F. Buda, T. E. O'Brien, and L. Visscher, A state-averaged orbital-optimized hybrid quantum-classical algorithm for a democratic description of ground and excited states, *Quantum Sci. Technol.* **6**, 024004 (2021).

[27] I. Kassal and A. Aspuru-Guzik, Quantum algorithm for molecular properties and geometry optimization, *J. Chem. Phys.* **131**, 224102 (2009).

[28] T. E. O'Brien, B. Senjean, R. Sagastizabal, X. Bonet-Monroig, A. Dutkiewicz, F. Buda, L. DiCarlo, and L. Visscher, Calculating energy derivatives for quantum chemistry on a quantum computer, npj Quantum Inf. **5**, 1 (2019).

[29] K. Mitarai, Y. O. Nakagawa, and W. Mizukami, Theory of analytical energy derivatives for the variational quantum eigensolver, *Phys. Rev. Research* **2**, 013129 (2020).

[30] I. O. Sokolov, P. K. Barkoutsos, L. Moeller, P. Suchsland, G. Mazzola, and I. Tavernelli, Microcanonical and finite-temperature *ab initio* molecular dynamics simulations on quantum computers, *Phys. Rev. Research* **3**, 013125 (2021).

[31] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard-even for logarithmically many qubits and free fermionic systems, arXiv:2101.07267.

[32] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, *Phys. Rev. A* **92**, 042303 (2015).

[33] N. C. Rubin, R. Babbush, and J. McClean, Application of fermionic marginal constraints to hybrid quantum algorithms, *New J. Phys.* **20**, 053020 (2018).

[34] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, *Nature (London)* **549**, 242 (2017).

[35] A. F. Izmaylov, T.-C. Yen, and I. G. Ryabinkin, Revising the measurement process in the variational quantum eigensolver: Is it possible to reduce the number of separately measured operators? *Chem. Sci.* **10**, 3746 (2019).

[36] V. Verteletskyi, T.-C. Yen, and A. F. Izmaylov, Measurement optimization in the variational quantum eigensolver using a minimum clique cover, *J. Chem. Phys.* **152**, 124114 (2020).

[37] A. Jena, S. Genin, and M. Mosca, Pauli partitioning with respect to gate sets, arXiv:1907.07859.

[38] A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, Unitary partitioning approach to the measurement problem in the variational quantum eigensolver method, *J. Chem. Theory Comput.* **16**, 190 (2019).

[39] T.-C. Yen, V. Verteletskyi, and A. F. Izmaylov, Measuring all compatible operators in one series of single-qubit measurements using unitary transformations, *J. Chem. Theory Comput.* **16**, 2400 (2020).

[40] P. Gokhale, O. Angiuli, Y. Ding, K. Gui, T. Tomesh, M. Suchara, M. Martonosi, and F. T. Chong, Minimizing state preparations in variational quantum eigensolver by partitioning into commuting families, arXiv:1907.13623.

[41] W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, Efficient and noise resilient measurements for quantum chemistry on near-term quantum computers, npj Quantum Inf. **7**, 23 (2021).
[42] X. Bonet-Monroig, R. Babbush, and T. E. O'Brien, Nearly Optimal Measurement Scheduling for Partial Tomography of Quantum States, Phys. Rev. X 10, 031064 (2020).

[43] A. Zhao, A. Tranter, W. M. Kirby, S. F. Ung, A. Miyake, and P. J. Love, Measurement reduction in variational quantum algorithms, Phys. Rev. A 101, 062322 (2020).

[44] G. Wang, D. E. Koh, P. D. Johnson, and Y. Cao, Minimizing estimation runtime on noisy quantum computers, PRX Quantum 2, 010346 (2021).

[45] G. García-Pérez, M. A. C. Rossi, B. Sokolov, F. Tacchino, P. K. Barkoutsos, G. Mazzola, I. Tavernelli, and S. Maniscalco, Learning to measure: Adaptive informationally complete povms for near-term quantum algorithms, arXiv:2104.00569.

[46] E. Campbell, Random Compiler for Fast Hamiltonian Simulation, Phys. Rev. Lett. 123, 070503 (2019).

[47] A. M. Childs and N. Wiebe, Hamiltonian simulation using linear combinations of unitary operations, Quantum Inf. Comput. 12, 901 (2012).

[48] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Exponential Improvement in Precision for Simulating Sparse Hamiltonians (Association for Computing Machinery, New York, 2014), p. 283.

[49] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).

[50] D. W. Berry, C. Gidney, M. Motta, J. R. McClean, and R. Babbush, Qubitization of arbitrary basis quantum chemistry leveraging sparsity and low rank factorization, Quantum 3, 208 (2019).

[51] P. Jordan and E. P. Wigner, About the Pauli exclusion principle, Z. Phys. 47, 631 (1928).

[52] J. Lee, D. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, Even more efficient quantum computations of chemistry through tensor hypercontraction, PRX Quantum 2, 030305 (2021).

[53] T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (John Wiley & Sons, Ltd., Chichester, 2000).

[54] J. R. McClean, R. Babbush, P. J. Love, and A. Aspuru-Guzik, Exploiting locality in quantum computation for quantum chemistry, J. Phys. Chem. Lett. 5, 4368 (2014).

[55] R. Babbush, D. W. Berry, I. D. Kivlichan, A. Y. Wei, P. J. Love, and A. Aspuru-Guzik, Exponentially more precise quantum simulation of fermions in second quantization, New J. Phys. 18, 033032 (2016).

[56] M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, and G. K. Chan, Low rank representations for quantum simulation of electronic structure, arXiv:1808.02625.

[57] J. D. Whitfield, J. Biamonte, and A. Aspuru-Guzik, Simulation of electronic structure Hamiltonians using quantum computers, Mol. Phys. 109, 735 (2011).

[58] Y. Ouyang, D. R. White, and E. T. Campbell, Compilation by stochastic Hamiltonian sparsification, Quantum 4, 235 (2020).

[59] K. Kuroiwa and Y. O. Nakagawa, Penalty methods for a variational quantum eigensolver, Phys. Rev. Research 3, 013197 (2021).

[60] S. Saebo and P. Pulay, Fourth-order Mo/ller-Plesset perturbation theory in the local correlation treatment. I. Method, J. Chem. Phys. 86, 914 (1987).

[61] S. Saebo and P. Pulay, Local treatment of electron correlation, Ann. Rev. Phys. Chem. 44, 213 (1993).

[62] M. Schütz, G. Hetzer, and H.-J. Werner, Low-order scaling local electron correlation methods. I. Linear scaling local MP2, J. Chem. Phys. 111, 5691 (1999).

[63] V. Weijo, P. Manninen, P. Jørgensen, O. Christiansen, and J. Olsen, General biorthogonal projected bases as applied to second-order Møller-Plesset perturbation theory, J. Chem. Phys. 127, 074106 (2007).

[64] C. Hampel and H.-J. Werner, Local treatment of electron correlation in coupled cluster theory, J. Chem. Phys. 104, 6286 (1996).

[65] M. Schütz and H.-J. Werner, Low-order scaling local electron correlation methods. IV. Linear scaling local coupled-cluster (LCCSD), J. Chem. Phys. 114, 661 (2001).

[66] T. Korona and H.-J. Werner, Local treatment of electron excitations in the EOM-CCSD method, J. Chem. Phys. 118, 3006 (2003).

[67] O. Christiansen, P. Manninen, P. Jørgensen, and J. Olsen, Coupled-cluster theory in a projected atomic orbital basis, J. Chem. Phys. 124, 084103 (2006).

[68] D. Maynau, S. Evangelisti, N. Guihéry, C. J. Calzado, and J.-P. Malrieu, Direct generation of local orbitals for multireference treatment and subsequent uses for the calculation of the correlation energy, J. Chem. Phys. 116, 10060 (2002).

[69] C. Angeli, C. J. Calzado, R. Cimiraglia, S. Evangelisti, N. Guihéry, T. Leininger, J.-P. Malrieu, D. Maynau, J. V. P. Ruiz, and M. Sparta, The use of local orbitals in multireference calculations, Mol. Phys. 101, 1389 (2003).

[70] N. Ben Amor, F. Bessac, S. Hoyau, and D. Maynau, Direct selected multireference configuration interaction calculations for large systems using localized orbitals, J. Chem. Phys. 135, 014101 (2011).

[71] C. Chang, C. J. Calzado, N. B. Amor, J. S. Marin, and D. Maynau, Multi-scale multireference configuration interaction calculations for large systems using localized orbitals: Partition in zones, J. Chem. Phys. 137, 104102 (2012).

[72] A. S. P. Gomes and C. R. Jacob, Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 108, 222 (2012).

[73] S. J. Bennie, M. Stella, T. F. Miller III, and F. R. Manby, Accelerating wavefunction in density-functional-theory embedding by truncating the active basis set, J. Chem. Phys. 143, 024105 (2015).

[74] B. Hégely, P. R. Nagy, G. G. Ferenczy, and M. Kállay, Exact density functional and wave function embedding schemes based on orbital localization, J. Chem. Phys. 145, 064107 (2016).

[75] S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, A practical guide to density matrix embedding theory in quantum chemistry, J. Chem. Theory Comput. 12, 2706 (2016).

[76] F. Libisch, M. Marsman, J. Burgdörfer, and G. Kresse, Embedding for bulk systems using localized atomic orbitals, J. Chem. Phys. 147, 034110 (2017).

[77] D. V. Chulhai and J. D. Goodpaster, Improved accuracy and efficiency in quantum embedding through absolute localization, J. Chem. Theory Comput. 13, 1503 (2017).

[78] E. R. Sayfutyarova, Q. Sun, G. K.-L. Chan, and G. Knizia, Automated construction of molecular active spaces from atomic valence orbitals, J. Chem. Theory Comput. 13, 4063 (2017).
[79] S. J. R. Lee, M. Welborn, F. R. Manby, and T. F. Miller III, Projection-based wavefunction-in-DFT embedding, *Acc. Chem. Res.* **52**, 1359 (2019).

[80] X. Wen, D. S. Graham, D. V. Chulhai, and J. D. Goodpaster, Absolutely localized projection-based embedding for excited states, *J. Chem. Theory Comput.* **16**, 385 (2020).

[81] D. Claudino and N. J. Mayhall, Automatic partition of orbital spaces based on singular value decomposition in the context of embedding theories, *J. Chem. Theory Comput.* **15**, 1053 (2019).

[82] D. Claudino and N. J. Mayhall, Simple and efficient truncation of virtual spaces in embedded wave functions via concentric localization, *J. Chem. Theory Comput.* **15**, 6085 (2019).

[83] M. R. Hermes and L. Gagliardi, Multiconfigurational self-consistent field theory with density matrix embedding: The localized active space self-consistent field method, *J. Chem. Theory Comput.* **15**, 972 (2019).

[84] M. R. Hermes, R. Pandharkar, and L. Gagliardi, The variational localized active space self-consistent field method, *J. Chem. Theory Comput.* **16**, 4923 (2020).

[85] I. Mayer, On Löwdin's method of symmetric orthogonalization, *Int. J. Quantum Chem.* **90**, 63 (2002).

[86] D. Szczepanik and J. Mrozek, On several alternatives for löwdin orthogonalization, *Comput. Theor. Chem.* **1008**, 15 (2013).

[87] P. Löwdin, On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, *J. Chem. Phys.* **18**, 365 (1950).

[88] B. C. Carlson and J. M. Keller, Orthogonalization procedures and the localization of Wannier functions, *Phys. Rev.* **105**, 102 (1957).

[89] J. Pipek and P. G. Mezey, A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions, *J. Chem. Phys.* **90**, 4916 (1989).

[90] J. M. Foster and S. F. Boys, Canonical configurational interaction procedure, *Rev. Mod. Phys.* **32**, 300 (1960).

[91] C. Edmiston and K. Ruedenberg, Localized atomic and molecular orbitals, *Rev. Mod. Phys.* **35**, 457 (1963).

[92] R. C. Raffenetti, K. Ruedenberg, C. L. Janssen, and H. F. Schaefer, Efficient use of Jacobi rotations for orbital optimization and localization, *Theor. Chim. Acta* **86**, 149 (1993).

[93] R. Barr and H. Basch, Improved convergence in orbital localization methods, *Chem. Phys. Lett.* **32**, 537 (1975).

[94] S. Lehtola and H. Jónsson, Unitary optimization of localized molecular orbitals, *J. Chem. Theory Comput.* **9**, 5365 (2013).

[95] J. M. Leonard and W. L. Luken, Quadratically convergent calculation of localized molecular orbitals, *Theoretica chimica acta* **62**, 107 (1982).

[96] I.-M. Høyvik, B. Jansik, and P. Jørgensen, Trust region minimization of orbital localization functions, *J. Chem. Theory Comput.* **8**, 3137 (2012).

[97] Q. Sun, Co-iterative augmented hessian method for orbital optimization, *arXiv:1610.08423*.

[98] F. Jensen, *Introduction to Computational Chemistry* (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2006).

[99] B. Jansík, S. Høst, K. Kristensen, and P. Jørgensen, Local orbitals by minimizing powers of the orbital variance, *J. Chem. Phys.* **134**, 194104 (2011).

[100] I.-M. Høyvik, B. Jansik, and P. Jørgensen, Orbital localization using fourth central moment minimization, *J. Chem. Phys.* **137**, 224114 (2012).

[101] S. Lehtola and H. Jónsson, Pipek–mezey orbital localization using various partial charge estimates, *J. Chem. Theory Comput.* **10**, 642 (2014).

[102] G. Knizia, Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts, *J. Chem. Theory Comput.* **9**, 4834 (2013).

[103] B. Senjean, S. Sen, M. Repisky, G. Knizia, and L. Visscher, Generalization of intrinsic orbitals to Kramers-paired quater-nion spinors, molecular fragments, and valence virtual spinors, *J. Chem. Theory Comput.* **17**, 1337 (2021).

[104] G. te Velde, F. M. Bickelhaupt, E. J. Baerends, C. Fonseca Guerra, S. J. A. van Gisbergen, J. G. Snijders, and T. Ziegler, Chemistry with ADF, *J. Comput. Chem.* **22**, 931 (2001).

[105] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.3.033127 for the geometries of molecules studied in this work.

[106] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, Pyscf: The python-based simulations of chemistry framework, *WIREs Comput. Mole. Sci.* **8**, e1340 (2018).

[107] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey et al., SciPy 1.0: Fundamental algorithms for scientific computing in Python, *Nat. Methods* **17**, 261 (2020).

[108] F. Gerber and R. Furrer, *optimParallel: An R package providing a parallel version of the L-BFGS-B optimization method*, *R J.* **11**, 352 (2019).

[109] J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby et al., *OpenFermion: The electronic structure package for quantum computers*, *Quantum Sci. Technol.* **5**, 034014 (2020).

[110] F. A. Evangelista, Garnet Kin-Lic Chan, and G. E. Scuseria, Exact parameterization of fermionic wave functions via unitary coupled cluster theory, *J. Chem. Phys.* **151**, 244112 (2019).

[111] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, and A. K. Wilson, Basis-set convergence in correlated calculations on Ne, N₂, and H₂O, *Chem. Phys. Lett.* **286**, 243 (1998).

[112] M. Yoshimine, Construction of the Hamiltonian matrix in large configuration interaction calculations, *J. Comput. Phys.* **11**, 449 (1973).

[113] M. Feyereisen, G. Fitzgerald, and A. Komornicki, Use of approximate integrals in ab initio theory. an application in MP2 energy calculations, *Chem. Phys. Lett.* **208**, 359 (1993).

[114] V. Havlíček, M. Troyer, and J. D. Whitfield, Operator locality in the quantum simulation of fermionic models, *Phys. Rev. A* **95**, 032332 (2017).