The Hartree–Fock problem provides the conceptual and mathematical underpinning of a large portion of quantum chemistry. As efforts in quantum technology aim to enhance computational chemistry algorithms, the Hartree–Fock method, central to many other numerical approaches, is a natural target for quantum enhanced algorithms. While quantum computers and quantum simulation offer many prospects for the future of modern chemistry, the non-deterministic polynomial-complete Hartree–Fock problem is not a likely candidate. We highlight this fact from a number of perspectives including computational complexity, practical examples, and the full characterization of energy landscapes for simple systems.

## I. INTRODUCTION

The present study aims to highlight the difficulty in conducting optimization problems in the context of electronic structure, with and without quantum resources. Specifically, we will focus on the Hartree–Fock problem, an optimization problem using the mean-field approximation. The Hartree–Fock problem^{1,2} provides the mathematical setting for molecular orbitals widely used in chemistry and beyond. The ubiquitous self-consistent field (SCF) methodology used to solve the Hartree–Fock problem is also applied to most implementations of density functional theory based on Kohn–Sham theory.^{3} While the solution to an instance of the Hartree–Fock problem is often insufficient for many applications, it often serves as the reference state for post-Hartree–Fock methods. Widely used post-Hartree–Fock methods include the coupled cluster ansatz,^{4–9} Møller–Plesset perturbation theory,^{10–13} equation of motion,^{14–20} multi-reference configuration interaction,^{21–23} and many more. Here, we forgo improving the Hartree–Fock ansatz and instead ask how difficult it is to find the true Hartree–Fock global minimum and its importance.

Many instances of the Hartree–Fock problem can be solved quickly using heuristic approaches. In practice, conventional algorithms for Hartree–Fock scale cubically with the number of basis functions. However, this cost only reflects the correct scaling if the number of iterations is bounded by a constant or if local minima are acceptable in place of a global solution. Linear scaling methods^{24,25} avoid the diagonalization of Hessian entirely and rely on localized properties of the system. This assumption may not be true generally for every case. Iterative procedures, regardless of cost per step, are prone to convergence issues. This article highlights the above-mentioned properties in the context of standard numerical and hybrid-quantum approaches to the Hartree–Fock problem.

Typical approaches to improve SCF convergence use direct inversion of the iterative subspace^{26} (DIIS), level shifting,^{27} quadratically convergent Newton–Raphson techniques,^{28} variational coupled cluster methods,^{29} or varying fractional occupation numbers,^{30} among many other approaches. The success of any above-mentioned methods varies at each instance and depends on the initial parameters chosen, e.g., the size of the iterative subspace, and often work well in combination, leading to attempts to build black-box SCF procedures.^{31,32}

Unfortunately, these methods cannot work in all cases without violating fundamental assumptions on the theory of computation. In previous work, arbitrary spin-glasses have been mapped to instances of the Hartree–Fock problem.^{33} In addition to other works,^{34,35} this shows that the Hartree–Fock problem is difficult in the worst case setting and non-deterministic polynomial (NP)-complete.^{36} The complexity class of non-deterministic polynomial (NP) time problems is the set of problems^{37} that can be solved efficiently with a hint.^{38} It is possible that every problem that can be solved efficiently with a hint can also be efficiently solved without the hint. However, all signs of practical experience suggest that problems in the NP class do not admit efficient black-box solutions. Thus, due to the NP-completeness of the Hartree–Fock problem, it is unlikely that any classical algorithm can solve all instances in a time proportional to a polynomial of the input size. Note that because it is in the NP complexity class, the Hartree–Fock problem can be solved efficiently in polynomial time with a sufficiently strong hint, e.g., as gleamed from experience or luck. It has been known since the late 1990s that quantum computers can promise no more than a quadratic speedup over their classical counterparts on such NP-complete problems.^{39}

It has long been known that any classical algorithm can be simulated using quantum hardware. With the heavy reliance of the variational quantum eigensolvers on classical optimization routines, it is marginal how the quantum approach differs from the conventional approaches. In the recent variational quantum eigensolver (VQE)^{40,41} study, the optimization strategy used was an augmented Hessian approach^{42} rather than the typical DIIS.^{26} However, both of these optimization strategies can be employed by conventional computers to optimize the orbitals.

In this article, we consider how the use of quantum hardware enhances the ability of chemists to solve instances of the Hartree–Fock problem. Quantum hardware has many prospects for applications to physical and chemical simulations;^{43} however, the Hartree–Fock problem is not likely to admit drastic advances using quantum computers.

The present article is inspired by a recent study of the Google group and collaborators^{40} using a variational quantum–classical hybrid formulation of the Hartree–Fock problem. There, the authors pointed out that the primary purpose of considering the Hartree–Fock problem was to benchmark their quantum device. Here, we highlight obstructions to its use as a general purpose replacement for standard SCF solvers. We do so by first introducing the Hartree–Fock problem and characterizations of the solution landscape for an instance of the Hartree–Fock problem. We then give examples that are (1) simple, (2) small, and (3) well-motivated instances of the Hartree–Fock problem that present convergence problems for black-box approaches.

## II. HARTREE–FOCK THEORY

The Hartree–Fock ansatz is important for conventional quantum chemistry and has been thoroughly developed over the past 100 years.^{44}

In this subsection, we wish to succinctly sketch the key ideas of Hartree–Fock theory. First, we give a standard presentation in terms of orbital rotations, and then, we reformulate the problem using bond density matrices as a useful representation of the problem.

### A. Orbital rotations

The Hartree–Fock optimization problem refers to the collection of instances each specific by a given electronic Hamiltonian characterized by one- and two-electron integrals ({*h*_{ij}}, {*h*_{ijkl}}). Once an instance is specified, the Hartree–Fock problem can be stated succinctly as minimizing the electronic energy in the space of single Fock states,

The space of single Fock states $F1$ is the set of all rank 1 [vs rank $MN$ in the general case] *N*-electron Fock states. Here and throughout, *N* is the number of electrons and *M* is the number of spatial basis functions. Each single Fock state is of the form

where we have $[bi,bj\u2020]+=bibj\u2020+bj\u2020bi=\delta ij$ and $[bj,bk]+=0$. The vacuum state with no particles is denoted by |Ω⟩ and is characterized by *a*_{j}|Ω⟩ = 0 for all *j*. Alternative formulations of the Hartree–Fock problem can be done using Slater determinants.^{45}

After the restriction to this ansatz, the variational space is characterized by rotations of the orbitals from an initial set of orbitals (e.g., the atomic orbitals). If the initial set corresponds to operators ${aj:[ai,aj\u2020]+=\u222bdx\u2009\varphi i(x)\varphi j*(x)=Sij,[ai,aj]+=0}$, then

Here, we require that *W*^{†}*W* = **1** when *S* = **1** and $W\u0303\u2020SW\u0303=1$ in the general case. For applications to molecular physics, single-electron spin-orbitals *ϕ*_{i}(*x*) corresponding to the fermionic operators $aj\u2020$ and *a*_{j} are often contracted Gaussians, e.g., STO-3G.^{46} If the orbitals in Eq. (3) yield the minimal energy, then they are called the molecular orbitals.

### B. Rotation of charge density matrix

The optimization problem stated in Eq. (1) can be solved using any of a variety of techniques. The SCF approach to instances of the Hartree–Fock problem is a direct application of the steepest descent approach, but other methods such as semi-definite programing can be used.^{47} In practice, the SCF method to solve Eq. (1) uses an effective potential term that takes into account the averaged two-body interaction (mean-field). In this approach, the *N*-body problem is reduced to a non-linear single particle problem. At each iteration of the simplest implementation, the Fock matrix, *F*, is formed as a function of a previous bond density matrix (*P*_{prev}), also known as one-particle reduced density matrix or the two-point correlation function,

with *η* as the orbital occupancies written as a diagonal matrix. The new transformation matrix, *W*, is determined using the gradient of Eq. (1) with respect to the bond density matrix *P*. The new coefficient matrix is used to form a new bond density matrix. At each iteration, the (real-valued) bond density matrix satisfies the following three properties:^{48}

Next, we consider transformation between bond density matrices. A convenient parameterization of the set of bond density matrices that avoids redundancy is

with *A*_{block} = *P*_{0}*A*(**1** − *P*_{0}) + (**1** − *P*_{0})*AP*_{0} for arbitrary skew-symmetric *A* = −*A*^{T}. In Eq. (8), *W*_{A} = exp(−*A*_{block}) is a change of basis as in Eq. (3).

The objective for the optimization problem in Eq. (1) is also given by the energy matrix functional^{49}

Here, the mean-field, *G*, is a function of the bond density matrix

with the antisymmetrized two-electron integral over spatial degrees of freedom defined by $h\mu \kappa \lambda \nu A=h\mu \kappa \lambda \nu \u2212h\mu \kappa \nu \lambda $. Note that Eq. (9) follows directly from evaluating the energy’s expectation value in Eq. (1) using, e.g., Slater–Condon rules.^{50} The matrix derivative of Eq. (9) gives the Fock operator

By rotating the charge density matrix with all possible rotations, we are able to do brute force exploration of the whole space for some small examples. Below, we expand Eq. (8) about *P*_{0} = *P*_{core} where the *N* lowest eigenmodes of $Hcore=\u2211hijai\u2020aj$ are occupied.

It is important to note that the global minimum of the functional equation (9) has significance. If one wishes to obtain approximations to the excited states of the Hamiltonian, Eq. (9) should not be used; instead, excited state functionals^{51} can be derived using excited state mean-field theory.^{52}

Before turning to examples, we will introduce the quantum ansatz for the Hartree–Fock method.

## III. QUANTUM CIRCUIT ANSATZ FOR HARTREE–FOCK

The quantum circuit for creating the Hartree–Fock ansatz state^{53} was applied in the context of VQE.^{40} The description of the ansatz circuit will be aided through the use of QR decomposition.^{54}

The Givens rotations provide a useful canonical characterization of an arbitrary orthogonal matrix, *W*. The QR decomposition of a real *n* × *n* orthogonal matrix *W* can be done using *T* = *n*(*n* − 1)/2 Givens rotations such that

When *W* has a determinant of 1, *D* is just the identity matrix. Each Givens rotation, *G*_{i}, is of the form *G*_{i} = *g*(*a*, *b*, *θ*) with *g*_{kk} = 1, unless *k* is either *a* or *b* when instead *g*_{kk} = cos(*θ*). All off-diagonal elements are zero, except *g*_{ab} = −*g*_{ba} = − sin(*θ*).

Applications of the Givens decomposition to fermionic orbital rotations has been worked out elsewhere,^{53,55} resulting in a quantum circuit that is able to prepare arbitrary Slater determinants following the parameters of the QR decomposition. By ordering the QR decomposition appropriately, a fermionic swap network can be used to rotate each pair of orbitals using the appropriate Givens rotation parameters. This results in an efficient state preparation circuit of the form depicted in Fig. 1. The full compilation down to gates including hardware optimization is given elsewhere.^{40,41}

Our characterization of the fermionic space in Eq. (8) gives us a set of parameters, Θ, that also characterizes the mixing between pairs of orbitals. The resulting orthogonal transformation *W*(Θ) is then given to the QR decomposition and forwarded to the quantum circuit construction.

## IV. CALCULATIONS AND RESULTS

All calculations of the molecular system are done in the STO-3G basis.^{46} Energies are reported in Hartrees, angles of rotation are in radians, and bond lengths are in angstroms.

The Hartree–Fock energy surfaces (HESs) were computed using PySCF.^{56,57} In this paper, we only consider Restricted Hartree–Fock (RHF) solutions where the alpha and beta spatial orbitals are restricted to be identical. The quantum optimization routines were that of OpenFermion-Cirq,^{41} and we only modify the initial state routines and the input molecular data.^{58} Internally, the optimization of the quantum circuit is performed using a conjugate gradient as implemented in SciPy.^{59}

### A. Landscape analysis

We consider H_{2}, $H3+$ as minimal basis model systems whose Hartree–Fock instances we can completely characterize. We begin with the H_{2} example.

When considering H_{2} in the minimal basis, there is only a single orbital mixing parameter. In Fig. 2, we have plotted the 1D HES surface as a function of bond length for H_{2}. The number of minima in HES(*θ*) changes with bond length. Before a bond length of ∼1.2 Å, there is only a single minimum. However, at larger bond lengths, an additional minimum begins to appear.

We continue with our two electron examples with the iso-electronic $H3+$. Now, instead of a 1D HES, we now have two parameters that mix the one occupied orbital with the two virtual orbitals. We plot the HES in Fig. 3 for a linear configuration with hydrogen atoms separated by 2.5 Å. There are three minima for HES(*θ*_{1}, *θ*_{2}). In Fig. 4, we give the HES of $H3+$ at 4.36 Å where there are several minima with the same globally optimal value.

In both the cases of $H3+$ and H_{2}, there are single occupied (spatial) orbital and *m* virtual orbitals. For *m* = 2, this leads to an *A*_{block} generator of the form

The eigenvalues of this matrix are $\lambda \xb1={0,\xb1i\theta 12+\theta 22}$. Since the matrix exponential of *A*_{block} merely exponentiates the eigenvalues, when *λ*_{±} = *iπ*, the rotation acts trivially on the density matrix. This underlies periodicity to the plots seen in Figs. 3 and 4.

We can explain the periodicity in terms of this invariant by converting to polar coordinates where *θ*_{1} = *R* cos *ϕ* and *θ*_{2} = *R* sin *ϕ*. Now, the nontrivial eigenvalue is *λ*_{±} = ±*iR*, and we can express the periodicity of the plots as HES(*θ*_{1}, *θ*_{2}) = HES(*R*, *ϕ*) = HES(*R* + *nπ*, *ϕ*), with *n* being an integer.

There is a nice generalization of this fact. For a single spatial orbital that is doubly occupied with two electrons and *m* virtual orbitals, the generalization of Eq. (13) is

It is straightforward to calculate that the eigenvalues of this matrix are zero except $\lambda \xb1=\xb1i\theta 12+\theta 22+\theta 32+\cdots +\theta m2$. Following the same argument as in the *m* = 2 case above, we can say that

where $R=\theta 12+\cdots +\theta m2$.

The range of the search space for each *θ*_{j} is [0, 2*π*). Straightforward combination of parameter results is a parameter space for bond density matrices that is a hyper-cube in *m* dimensions with equal side lengths of 2*π*. However, due to the constraint of Eq. (15), the minimal search space is now confined to the interior of a hyper-sphere of radius *π* (also in *m* dimension). Therefore, the ratio of the minimal required search space to the full search goes to zero as *m* tends to infinity. This is a well-known consequence of the vanishing ratio of the volume of a hyper-sphere to the volume of the corresponding hyper-cube.^{60,61}

### B. Convergence analysis

We used the quantum algorithm outlined in Ref. 40 for obtaining RHF solutions for four examples. Depending on the initial guess, it may converge to local rather global solutions.

The different initial guesses were generated using the Givens rotations corresponding to different minima in Figs. 2 and 3, respectively. The results for H_{2}, converging to two different minima on the quantum simulator, are shown in Fig. 5. Similarly, results for $H3+$, converging to two different minima on the quantum simulator, are shown in Fig. 6. In Fig. 5, the values of *λ* = 1 and *λ* = 1.2 were chosen as states far enough in parameter space to have energy sufficiently large. If we select states with small *λ*, the convergence to the minimum is highly likely so long as the system does not climb uphill in energy since the initial state has energy less than all minima except the global minimum.

As final examples, we choose diatomic carbon and its cation. We also consider C_{2} and $C22+$ as instances that are commonly known to confound solvers due to the appearance of saddle points in the optimization landscape. To illustrate the complications of convergence, we modified the initial parameters following the method used in Ref. 40. Namely, we begin with the solution provided by the classical SCF solver, perturb from those solutions, and observe if the quantum algorithm still converges to the correct minimum.

Since the Hartree–Fock problem is an optimization problem, quantum computing via a modified Grover search^{39,62} can be used to find the global solution quadratically faster than the classical brute force approaches. However, in both conventional and quantum solvers, local searches are employed whereby no guarantees on finding the globally optimal solutions are given.

This example allows us to highlight the importance of using the information of Hessian to avoid saddle points during SCF optimization in either quantum or classical methods. In Figs. 7 and 8, convergence results are shown for C_{2} and $C22+$, respectively. In each of the plots, we have plotted the performance of the quantum circuit optimization routine and the classical optimization routines as implemented in PySCF. At an inter-nuclear separation of 1.5 Å, there is a saddle point appearing in the potential energy landscape.

The use of the orbital Hessian helps solvers avoid saddle points where the gradient may vanish at non-optimal values. The solver uses Hessian, which provides a notion of the curvature of the landscape. This allows the solver to avoid saddle points. The augmented Hessian conjugate-gradient method^{42} was used for the quantum optimization of the Hartree–Fock circuit. This allows the solver to avoid convergence to saddle points, as illustrated in Figs. 7 and 8. Many of the quantum chemistry do not check the RHF solutions by default. This makes them prone to failure by convergence to a saddle point rather than the RHF solution.

## V. DISCUSSION

In this study, we showed that the convergence to a local minimum or global minimum is a function of the initial guess. The proposed algorithm on the quantum simulator carries this feature from the classical algorithm. There is no *a priori* guarantee that it will find the global solution.

Since the Hartree–Fock problem is an optimization problem, quantum computing via a modified Grover search^{39,62} can be used to find the global solution quadratically faster than the classical brute force approaches. However, in both conventional and quantum solvers, local searches are employed whereby no guarantees on finding the globally optimal solutions are given.

In analyzing the Hartree–Fock functional for H_{2}, $H3+$, C_{2}, and $C22+$, we note that the number of critical points in the solution space^{63} changes as the nuclear separation changes. This also implies that there are additional difficulties in applying Hartree–Fock for nuclear dynamics or other non-equilibrium configurations. Most solvers for the Hartree–Fock problem do not explore the entire space and usually only choose a single starting point (rather than multiple starting points). While this is not often an issue, it can cause serious pathology when used in post-Hartree–Fock methods, e.g., the coupled cluster method. This will be true for both on conventional computers and in its unitary formulation for quantum computers.

## VI. CONCLUSIONS

The recent application of quantum technology to the Hartree–Fock problem may serve as a hardware benchmark but is unlikely to have a dramatic impact on the practical approaches to this problem.

While the application to the Hartree–Fock problem is not likely to change the workhorse routines used on conventional computers, there are still interesting use cases for the results from Ref. 40. The projection methods and purity extrapolation presented in Ref. 40 will still be useful.

We should point out that the obstructions to efficiently solving the Hartree–Fock problem using quantum resources is not directly related to the concentration of measure phenomena, leading to barren plateaus discussed in the context of VQE ansatz initialization.^{64} It would be interesting to combine the insights of computational complexity with similar arguments for the Hartree–Fock ansatz.

When considering application areas of quantum computers, it is far more likely to make major breakthroughs when considering time-dependent phenomena. For example, the quantum–classical hybrid algorithm for obtaining the Kohn–Sham potential of the time-dependent density functional^{65–67} also requires measuring the bond density matrix.

This article highlights the wide body of knowledge on the SCF method and its difficulty and shows that the lack of verification of the solution is material in both the quantum and conventional computing domains.

## ACKNOWLEDGMENTS

The authors thank Z. Zimboras and N. Anand for many useful discussions. J.D.W. acknowledges support from the Department of Energy under Grant No. DE-SC0019374 “Quantum Chemistry for Quantum Computers.” Additional support for J.D.W. comes from the NSF (Grant Nos. PHYS-1820747 and EPSCoR-1921199) and from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under programs Quantum Computing Application Teams and Accelerated Research for Quantum Computing. S.G. would like to thank USC for the Dornsife Graduate School Fellowship.

## DATA AVAILABILITY

The data that support the findings of this study are available from the authors upon reasonable request.

## REFERENCES

In this article, we use the term problem to refer to a specific collection of problem instances. That is to say, the Hartree–Fock problem is the collection of all instances of Hartree–Fock, i.e., Eq. (1) with a fixed choice {*h*_{ij}, *h*_{ijkl}}.

Each Fock state corresponds to a Slater determinant $\Psi J(x1,x2,\u2026xN)=\u2009\Omega |\psi ^(x1)\psi ^(x2)\u2026\psi ^(xN)|\Psi J\u2009/N!$ once orbitals {*ϕ*_{j}(*x*)} have been assigned to operators {*a*_{j}} via field modes $\psi ^(x)=\u2211j\varphi j(x)aj$.

If *S* ≠ **1**, then a bond density matrix in this space $P\u0303$ satisfies Tr$(P\u0303S)=N/2$ and $P\u0303=P\u0303SP\u0303$. We can convert to an orthogonal basis using an orthogonalizing transform *X*, where *X*^{†}*SX* = **1**, e.g., canonical orthogonalization $Xcan=US\u2020s\u22121/2$ with $S=USsUS\u2020$. Then, we orthogonalize as $P=X\u2020P\u0303X$.

Note that only the global minima are constructed to have any connection to the energies of the Hamiltonian.