Bootstrap Embedding on a Quantum Computer

To provide a foundation for a more concrete exposition of the bootstrap embedding method, we first establish some rigorous notation for discussing molecular Hamiltonians and their associated Hilbert spaces. We will work with the molecular Hamiltonian under the second quantization formalism. Specifically, given a particular molecule of interest, define $O=\{\phi_{\mu}\ |\ \mu=1,\ldots,N\}$ to be an orthonormal set of single-particle local orbitals (LOs), where $N$ is the total number of orbitals; in this work, these LOs are generated through Löwdin’s symmetric orthogonalization method ⁵¹. The full Hilbert space $\mathcal{H}$ for the entire molecular system is thus given by $\mathcal{H}=\mathcal{F}(O)$, where $\mathcal{F}(O)$ denotes the Fock space determined by the LOs in the set $O$. Further define the creation (annihilation) operator $c^{\dagger}_{\mu}$ ($c_{\mu}$) which creates (annihilates) an electron in the LO $\phi_{\mu}$, the molecular Hamiltonian is written in the second-quantized notation

where $h_{\mu\nu}$ and $V_{\mu\nu\lambda\sigma}$ are the standard one- and two-electron integrals.

Note that the number of terms in the full molecular Hamiltonian $\hat{H}$ scales polynomially with the total number of orbitals $N$, but the dimension of $\mathcal{H}$ scales exponentially with $N$. Clearly, for large $N$, it will become prohibitively expensive to directly compute the exact full ground state. To circumvent this issue, we divide the full molecule into multiple smaller fragments, each equipped with its own “embedding Hamiltonian” which contains a number of terms that only scales polynomially with the number of orbitals in the fragment. Given that there are potentially far fewer orbitals in each fragment than in the whole molecular system, computing the ground state of each fragment’s embedding Hamiltonian can be significantly less expensive than computing the ground state of the full system. Furthermore, using the bootstrap embedding procedure to be described later, the ground states of individual fragments can, to a high degree of accuracy, be algorithmically combined to recover the desired electron densities prescribed by the exact ground state of the full system. Thus, this combination of fragmentation and bootstrap embedding can be used to reconstruct the full molecular ground state more efficiently than by direct computation alone.

We now briefly review the construction of embedding Hamiltonians for each fragment. Consider a single fragment associated with a label $A$, without loss of generality, define $O^{(A)}=\{\phi_{\mu}\ |\ \mu=1,\ldots,N_{A}\}$ with $N_{A}\leq N$ to be the set of LOs contained in fragment $A$; we will refer to $O^{(A)}$ as the set of fragment orbitals. Note that $O^{(A)}\subseteq O$, the set of LOs for the entire molecular system. The construction of the embedding Hamiltonian $\hat{H}^{(A)}$ for fragment $A$ begins with any solution of the ground state of the full system $\hat{H}$. For simplicity, the Hartree-Fock (HF) solution $\ket{\Phi_{\textrm{HF}}}$ is often used because it is easy to obtain on a classical computer. By invoking a Schmidt decomposition, we can write $\ket{\Phi_{\textrm{HF}}}$ with the following tensor product structure for $\forall~{}A$

In the above decomposition, the $\ket{f_{i}^{(A)}}$ represent single-particle fragment states contained in the Fock space $\mathcal{F}(O^{(A)})$ of fragment orbitals. On the other hand, the $\ket{b_{i}^{(A)}}$ and $\ket{\Psi_{\textrm{env}}^{(A)}}$ represent Slater determinants contained in the “environment” Fock space $\mathcal{F}(O\setminus O^{(A)})$ of the $N-N_{A}$ orbitals not included in the fragment. The key difference between the single environment state $\ket{\Psi_{\textrm{env}}^{(A)}}$ and the various “bath” states $\ket{b^{(A)}_{i}}$ is that the bath states $\ket{b^{(A)}_{i}}$ are entangled with the fragment states $\ket{f_{i}^{(A)}}$ while $\ket{\Psi_{\textrm{env}}^{(A)}}$ is not; this entanglement is quantified by the Schmidt coefficients $\lambda_{i}^{(A)}$. Crucially, since the HF solution is used, the sum in Eq. (2) only has $N_{A}$ terms (as opposed to $2^{N_{A}}$ for a general many-body wave function). Denote the collection of the $N_{A}$ entangled bath orbitals as $O^{(A)}_{\rm bath}=\{\beta_{\mu}\ |\mu=1,\ldots,N_{A}\}$, where each of the LOs $\beta_{\mu}$ are linear combinations of the original LOs not included in the fragment, $\beta_{\mu}\in\operatorname{Span}\{O\setminus O^{(A)}\}$. Furthermore, we denote the Fock space that corresponds to this set of entangled bath orbitals as $\mathcal{F}(O^{(A)}_{\rm bath})$.

This tensor product structure of $\ket{\Phi_{\rm HF}}$ allows us to naturally decompose the Hilbert space $\mathcal{H}$ for the full molecular system into the direct product of two smaller Hilbert spaces, namely

where

is the active fragment embedding space and $\mathcal{H}^{(A)}_{\textrm{env}}$ contains the remaining states, including $\ket{\Psi_{\textrm{env}}^{(A)}}$. Note that since both sets $O^{(A)}$ and $O^{(A)}_{\rm bath}$ have size $N_{A}$, the fragment Hilbert space $\mathcal{H}^{(A)}$ is a Fock space spanned of just $2N_{A}$ single-particle orbitals. The core intuition motivating this decomposition is that, in the exact ground state of the full system, states in $\mathcal{H}^{(A)}_{\textrm{env}}$ are unlikely to be strongly entangled with the many-body fragment states (consider the approximate HF ground state in Eq. (2), where they are perfectly disentangled); therefore, in a mean-field approximation, it is reasonable to entirely disregard the states in $\mathcal{H}^{(A)}_{\textrm{env}}$ when calculating the ground state electron densities on fragment $A$. Following this logic, we can define an embedding Hamiltonian $\hat{H}^{(A)}$ for fragment $A$ only on the $2N_{A}$ LOs in $\mathcal{H}^{(A)}$, which will have the form

given some creation and annihilation operators $a^{(A)\dagger}_{p}$ and $a^{(A)}_{p}$, which respectively create and annihilate electrons in orbitals from the combined set $O^{(A)}\cup O^{(A)}_{\rm bath}$ for $\mathcal{H}^{(A)}$. The new one- and two- electron integrals $h^{(A)}_{pq}$ and $V^{(A)}_{pqrs}$ can be computed by projecting $\hat{H}$ into the smaller Hilbert space $\mathcal{H}^{(A)}$ (consult the Supporting Information (SI) Sec. LABEL:app:frag for details on constructing $h_{pq}^{(A)}$ and $V_{qprs}^{(A)}$). Note that since we can choose $2N_{A}\ll N$, the ground state of this embedding Hamiltonian can be solved at a significantly reduced cost when compared to that of the full system Hamiltonian.

We are hence prepared to generate an embedding Hamiltonian for any arbitrary fragment of the original molecular system. However, the ground state electron densities of the fragment embedding Hamiltonian are unlikely to exactly match those of the full system Hamiltonian because, as mentioned above, the embedding process may neglect some small (but nonzero) entanglement of the fragment orbitals with the environment. Because we can expect interactions in the molecular Hamiltonian to be reasonably local, we anticipate that the electron densities on orbitals near the edge of the fragment (those closest to the “environment”) will deviate most significantly from their true values, while electron densities on orbitals toward the center of the fragment will be most accurate. Note that the uneven distribution of entanglement in molecular systems may likely lead to the potential sensitivity of the BE results to particular choices and partitions of fragments ^{52, 53, 54, 55, 28, 56}, while how quantum computers may help to reduce such dependence is an open problem.

To improve the accuracy of the fragment ground state wave function near the fragment edge, we employ the technique of bootstrap embedding. Broadly speaking, we first divide the full molecule into overlapping fragments such that the edge of each fragment overlaps with the center of another. Fig. 1i illustrates this fragmentation strategy: for example, we see that the edge of fragment $A$ (labeled as orbital 3) coincides with the center of fragment $B$. We then apply additional local potentials to the edge sites of each fragment to match their electron densities to those on overlapping center sites of adjacent fragments. Because we expect the electron densities computed on the center sites to be closer to their true values, these added local potentials should improve the accuracy of each fragment wave function near the edges. In the next section, we will formalize this edge-to-center matching process rigorously and discuss its implementation on a classical computer.

As mentioned in the previous section, we intend to correct the electron density error near a fragment’s edge by applying a local potential to the edge; this local potential serves to match the edge electron density of the fragment to the center electron density of an adjacent overlapping fragment, which we expect to be more accurate. In principle, to achieve an exact density matching, all $k$-electron reduced density matrices ($k$-RDM, for any $k$) on the overlapping region have to be matched. However, in practice, such matching beyond the 2-RDM is difficult on a classical computer due to the mathematical challenge that the number of terms in $k$-RDM in general increases exponentially as $k$. In addition, almost all electronic structure codes available on classical computers are programmed to deal with only 1- and 2-RDMs, despite the importance of $k$-RDMs ($k>2$) for computing observables such as entropy and other multi-point correlation functions ⁵⁷. Due to this reason, the discussion of density matching process in classical BE in this section will be based on 1-RDMs. We note that the matching process applies similarly if $k$-RDMs are matched.

We begin by introducing some rigorous notation. Recall that a fragment $A$ is defined by a set of local orbitals $O^{(A)}$ which constitute the fragment. We partition this set of LOs into a subset of edge sites (or orbitals), denoted $\mathbb{E}^{(A)}$, and a subset of center sites, denoted $\mathbb{C}^{(A)}$, such that $\mathbb{E}^{(A)}\cup\mathbb{C}^{(A)}=O^{(A)}$ and $\mathbb{E}^{(A)}\cap\mathbb{C}^{(A)}=\emptyset$. Given the ground state wave function $\ket{\Psi^{(A)}}$ of the embedding Hamiltonian, we further define the 1-electron reduced density matrix (1-RDM) $\mathbf{P}^{(A)}$ according to

where $p,q=1,\ldots,2N_{A}$ and the operators $a^{(A)\dagger}_{p}$ and $a^{(A)}_{q}$ are defined in the previous section.

Suppose, for example, that the edge of fragment $A$ overlaps with the center of another fragment $B$ so that $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}\neq\emptyset$. On a high level, the goal of bootstrap embedding is to find a ground state wave function $\ket{\Psi^{(A)}}$, perturbed by local potentials on the edge sites of $A$, such that $|P^{(A)}_{pq}-P^{(B)}_{pq}|\rightarrow 0$ for indices $p$ and $q$ that correspond to orbitals in the set of overlapping sites $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}$. More generally, and more rigorously, the goal is to find a wave function which minimizes the fragment Hamiltonian energy

subject to the constraints

for all other fragments $B$ with $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}\neq\emptyset$ and for all $p,q$ corresponding to orbitals in $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}$. Here, we explicitly write the expectation $\langle\cdot\rangle_{A}=\bra{\Psi^{(A)}}\cdot\ket{\Psi^{(A)}}$ in terms of $\ket{\Psi^{(A)}}$ to indicate that the optimization is over the wave function of $A$.

We can formulate this constrained optimization problem as finding the stationary solution to a Lagrangian by associating a scalar Lagrange multiplier $(\lambda^{(A)}_{B})_{pq}$ to Eq. (8). Since Eq. (8) has to be satisfied for any $p,q$ and $B$ that overlaps with $A$, these constraint can be rewritten in a more compact vector form $\bm{\lambda^{(A)}_{B}}\cdot\bm{\mathcal{Q}_{\textbf{1-RDM}}}(\Psi^{(A)};\mathbf{P}^{(B)})$ where the dot product conceals the implicit sum over $p,q$, and each component of the vector $\mathcal{Q}_{\textrm{1-RDM}}(\Psi^{(A)};\mathbf{P}^{(B)})_{pq}$ represents the constraint associated with Lagrange multiplier $(\lambda^{(A)}_{B})_{pq}$, given by the left hand side of Eq. (8). With this notation, we arrive at the following Lagrangian with the constraint added as an additional term

where once again the $B$ are fragments adjacent to $A$ with $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}\neq\emptyset$ and $p,q$ are indices of orbitals contained in the overlapping set $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}$. Here, the additional constraint with Lagrange multiplier $\mathcal{E}^{(A)}$ is also included to ensure normalization of the ground state wave function $\ket{\Psi^{(A)}}$. Solving for the stationary solution of the Lagrangian in Eq. (9) will only result in a ground state wave function for fragment $A$ whose 1-RDM elements at the edge sites match those at the center sites of adjacent overlapping fragments. However, we would instead like to solve for such a ground state for all fragments in the molecule simultaneously. Toward this regard, we can combine all individual fragment Lagrangians (of the form of Eq. (9)) into a single composite Lagrangian for the whole molecule, given by

where $N_{\textrm{frag}}$ is the number of fragments in the molecule. Observe that we have added one additional constraint

with Lagrange multiplier $\mu$ to restore the desired total number of electrons in the molecule, $N_{e}$. Note in Eq. (11) that $p^{\prime}$ is summed over indices corresponding to orbitals only in $\mathbb{C}^{(A)}$; this is to ensure that there is no double-counting of electrons in the whole molecule. By self-consistently finding ground states $\ket{\Psi^{(A)}}$ for $A=1,\ldots,N_{\textrm{frag}}$ which make the composite Lagrangian in Eq. (10) stationary, we will have completed the density matching procedure for all fragments, and the process of bootstrap embedding will be complete.

We can gain insight into which wave functions $\ket{\Psi^{(A)}}$ will make the composite Lagrangian $\mathcal{L}$ stationary by differentiating $\mathcal{L}$ with respect to $\ket{\Psi^{(A)}}$ for some fixed fragment $A$ and setting the resulting expression equal to zero. Upon some algebraic manipulation, we can recover the eigenvalue equation

where $V_{\textrm{BE}}$, the local bootstrap embedding potential, is given by

where the $p,q$ are indices of orbitals in the overlapping set $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}$, and the $p^{\prime}$ are indices of orbitals in the fragment center $\mathbb{C}^{(A)}$. We see that, when the composite Lagrangian is made stationary with respect to the fragment wave functions, the bare fragment embedding Hamiltonians become dressed with a potential $V_{\textrm{BE}}$ that contains a component local to the edge sites of each fragment (see the left term of Eq. (13)). This observation confirms our intuition that adding a local potential to the edge of one fragment will allow the edge site electron density to be matched to that of a center site on an overlapping neighbor. Note that $V_{\rm BE}$ also contains an additional potential on the center sites of each fragment (see the right term of Eq. (13)); this is simply to conserve the total electron number in the molecule. Moreover, $V_{\rm BE}$ as in Eq. (13) only contains one-body terms because only 1-RDM is used for density matching. In general, $V_{\rm BE}$ will contain up to $k$-body terms if $k$-RDMs are used for matching.

On a classical computer, the composite Lagrangian in Eq. (10) is made stationary through an iterative optimization algorithm ²² until the edge-to-center matching condition for all fragments is satisfied by some criterion. One possible criterion is to terminate the algorithm when the root-mean-squared 1-RDM mismatch, given by

drops below some predetermined threshold. Note again that $p,q$ are indices corresponding to orbitals in the overlapping set $\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}$; also, $N_{\rm sites}$ denotes the total number of overlapping sites in the whole molecule, equal to $N_{\rm sites}=\sum_{A}^{N_{\rm frag}}\sum_{B}\sum_{p,q}1$. The final set of density-matched fragment wave functions $\{\ket{\Psi^{(A)}}\}$ for $A=1,\ldots,N_{\rm frag}$ which solve the composite Lagrangian can then be used to reconstruct the electron densities and other observables for the full molecular system, as desired.

Given the notation established for classical BE, we now begin presenting new material. We discuss the computational resource requirement and typical behaviors of performing BE on classical computers to set the stage for a quantum BE theory. The details of the classical BE algorithms are omitted for succinctness, and we refer the reader to Ref. ^{22, 23, 50, 24} for details.

The space and time resource requirement to perform the classical BE can be broken down into two parts: a) the number of iteration steps to reach a fixed accuracy for $\epsilon$ (Eq. (14)); b) the runtime of the fragment eigensolver. For a), numerical evidence suggests an exponentially fast convergence on total system energy as the number of bootstrap iteration increases (black trace in Fig. 2 for FCI), while a proof of the convergence rate has yet to be established.

We focus on resource requirement in b) in the following. Admittedly, an exact classical eigensolver such as full configuration interaction (FCI) can be used to solve the embedding Hamiltonian in Eq. (5). However, both the storage space and time requirement scales exponentially as the the number of orbitals (see blue symbols and dashed line in Fig. 3 for the runtime scaling of FCI). Even with the state-of-the-art classical computational resources, exact solutions using FCI are only tractable for systems up to 20 electrons in 20 orbitals ⁵⁸.

As a result, classical computation of BE resorts to approximate eigensolvers with only polynomial cost in practice, by properly truncating or sampling from the fragment Hilbert space. One example for truncation is the coupled-cluster singles and doubles (CCSD) ⁵⁹, which scales with $N^{6}$ with $N$ being the number of orbitals. Alternately, different flavors of stochastic electronic structure solvers can be employed as fragment solvers in BE. Depending on implementation, these stochastic solvers can be biased or unbiased (if unbiased, with a cost of introducing the phase problem in general) ^{60, 61, 62, 63}. Collecting each sample on a classical computer usually has similar cost as a mean field theory (roughly $O(N^{3})$), while the overall target accuracy $\epsilon$ on observable estimation can be achieved with a sampling overhead of roughly $O(\frac{1}{\epsilon^{2}})$ with a constant prefactor depending on the severity of the sign problem.

Importantly, the sampling feature of these stochastic electronic structure methods on classical computers are strikingly similar to the nature of quantum computers where measurement necessarily collapses the wave function. As a result, only a classical sample (in terms of measurement results) can be obtained from a quantum computer. This similarity suggests a general strategy that many sampling techniques in stochastic classical algorithms can be deployed to design better quantum algorithms. For example, sophisticated importance sampling techniques ^{64, 65} can be employed to greatly improve the sampling efficiency in both classical ⁶³ and quantum cases ⁶⁶.

Due their shared feature on sampling between classical stochastic algorithm and quantum eigensolvers, we shall use one approximate sign-problem-free flavor of stochastic electronic structure method, the variational Monte Carlo (VMC), to serve as an additional baseline scenario for comparison with quantum BE in later sections. In addition to BE convergence behavior with a FCI solver, Fig. 2 also shows, for a VMC eigensolver with single Slater-Jastrow type wave function with two-body Jastrow factors ^{67, 68}, the density mismatch converges exponentially fast initially as iteration number increases with varying number of samples. However, partially due to the statistical noise on estimating the 1-RDM (thus the gradient for the optimization), the final density mismatch plateaus to a finite biased value. Comparing the VMC results across different numbers of samples, we can see that the bias improves as the number of samples increases (dashed horizontal lines). It is also evident that the orange trace (640k samples) has smaller fluctuation as compared to blue (160k samples) and grey (40k samples). Strictly speaking, an increase of sample size by a factor of 16 would decrease the statistical fluctuations by a factor of 4. However, our numerical data in Fig. 2 only shows qualitative but not quantitative agreement with this statement. We attribute part of the bias in the plateau to the intrinsic truncation of the VMC ansatz in addition to statistical fluctuations.

The increasing accuracy of density mismatch with respect to BE iteration also suggests an increasing number of samples are needed. Thus, an optimal number of samples at each BE iteration must be determined to achieve the desired accuracy in the matching conditions. A careful design of such a sampling schedule can potentially save a large amount of computational resources. We defer a thorough discussion of this point to later sections on quantum BE.

By employing the coherent superposition and entanglement of quantum states, the limitation of an exact classical solver can be overcome by substituting it with an exact quantum eigensolver such as the quantum phase estimation (QPE) algorithm ³¹. This section directly compares the cost between the two exact eigensovlers on quantum and classical computers, the QPE and the FCI solvers, using hydrogen chains where the initial trial state with a non-vanishing overlap with the exact eigenstate for QPE can be efficiently prepared on classical computers.

Fig. 3 compares the runtime (gate depth) of FCI and QPE for finding the ground state of linear hydrogen chain H_n for different system size $n$. Clearly, the QPE runtime scales only polynomially as the system size increases as expected ^{32, 30}, while its classical counterpart (FCI) has an exponentially increasing runtime. Note the runtime is normalized to the case of $n=1$ for each solver separately (see SI Sec. LABEL:app:computational-details). The dramatic advantage in the runtime scaling of quantum over classical eigensolvers demonstrated above suggests formulating BE on a quantum computer can bring significant benefits.

One might think that the eigensolver at the heart of the classical BE algorithm could simply be replaced with a quantum one. However, as mentioned before, there are two outstanding challenges for such a quantum bootstrap embedding (QBE) method. First, just as in classical stochastic methods, the results of a quantum eigensolver need to be measured for later use, but quantum wave functions collapse after measurement. Therefore, sampling from the quantum eigensolver is required, and the optimal sampling strategy is unclear. Secondly, with quantum wave function from quantum eigensolvers, it is not wise to achieve matching between fragments in the same way as classical BE, as many incoherent samples are needed to obtain a good estimation of the 1-RDM elements. Clearly, performing matching in a quantum way is desired.

In the next two sections (Secs. 3 and 4), we present how we address these two challenges by an adaptive quantum sampling scheduling algorithm and a quantum coherent matching algorithm in detail.

When mapping electronic structure problem to qubits on quantum computers, it is well-known that the global anti-symmetric property of fermionic wave functions necessarily leads to an overhead in operator lengths or qubit counts ⁶⁹. On the other hand, chemical information is usually local if represented using localized single-particle orbitals ^{70, 71}. In the case of performing bootstrap embedding, this tension between locality of chemical information and global fermionic anti-symmetry is more subtle. Because bootstrap embedding intrinsically uses the fermionic occupation number in the local orbitals (LOs) to perform matching, it is therefore convenient to preserve such locality when constructing the mapping. Throughout the discussion, without loss of generality, we assume a mapping that preserves fermionic local occupation number, such as the Jordan-Wigner mapping where each spin-orbital is mapped to one qubit. Our discussion equally applies to cases where a non-local mapping is used (such as parity mapping). In that case, a unitary transformation from the non-local mapping to a local mapping will be required before actually computing the matching conditions.

It is possible to formulate QBE using matching conditions on either qubit reduced density matrices (RDMs) ⁷² or $k$-electron RDMs ⁷³ for all $k$, both with an exponential number of matrix elements. For simplicity, in the present work we use qubit RDMs in our QBE and leave an efficient formulation in terms of fermionic $k$-electron RDMs for future work. The full density matrix of fragment $A$ is thus provided by $\rho^{(A)}=\ket{\Psi_{A}}\bra{\Psi_{A}}$. Given an orbital set $R\subset O^{(A)}$ for $O^{(A)}$ being set of orbitals in fragment $A$. Let $\rho_{R}^{(A)}$ signify the RDM obtained from $\rho^{(A)}$ by tracing out the set of qubits not in $R$. Specially, if $R$ only contains orbitals on the edge (center) of fragment $A$, then $\rho_{R}^{(A)}$ represents information about the density information (for example the occupation number) on the edge (center) of $A$.

These RDMs can be expanded under an arbitrary set of orthonormal basis $\{\Sigma_{\alpha}\}$ as follows

where $\langle\Sigma_{\alpha}\rangle_{A}=\bra{\Psi_{A}}\Sigma_{\alpha}\ket{\Psi_{A}}=\Tr[\rho^{(A)}~{}\Sigma_{\alpha}],~{}\forall\alpha\in[1,4^{m}-1]$, and $m=|R|$ is the number of orbitals in the set $R$. One convenient orthonormal basis set is the generalized Gell-Mann basis ⁷⁴. In the special case of a 1-qubit RDM, $\{\Sigma_{\alpha}\}$ ($\alpha=x,y,z$) is the familiar Pauli matrices.

A naive implementation of BE on a quantum computer is to simply replace 1-RDM in Eq. (6) with the qubit RDM in Eq. (15) on the fragment overlapping regions. Such an extension imposes matching constraints on each elements of the RDMs, resulting the following constraint vector in analogous to Eq. (8)

It is obvious that $\rho_{R}^{(A)}-\rho_{R}^{(B)}=0$, if and only if all the $(4^{m}-1)$ components in the above constraint are satisfied.

Similarly, we can associate a scalar Lagrange multiplier to each constraint in Eq. (16) and use this linear RDM constraint in place of the 1-RDM constraint $\bm{\mathcal{Q}_{\textbf{1-RDM}}}(\Psi^{(A)};\mathbf{P}^{(B)})$ in Eq. (9). Finding the stationary point of this new Lagrangian gives the same eigenvalue equation as Eq. (12) with a new BE potential given by

where $\bm{\Sigma_{r}}=\begin{bmatrix}\Sigma_{1},\cdots,\Sigma_{\alpha},\cdots,\Sigma_{4^{m}-1}\end{bmatrix}$ is a $(4^{m}-1)$-dimensional vector of the orthonormal basis in Eq. (15), and $\bm{\lambda_{B}^{(A)}}$ is the Lagrange multipliers now modulating the local potentials on each qubit basis, and $n$ is the number of overlapping sites between $A$ and $B$.

To perform the optimization, the eigenvalue equation Eq. (12) with the above new BE potential in (17) can be solved on a quantum computer to obtain an updated wave function for fragment $A$. By iteratively solving the eigenvalue equation and updating the Lagrange multipliers $\{\bm{\lambda},\mu\}$ using either gradient-based or gradient-free methods ⁷⁵, an algorithm can be formulated to solve the optimization problem. For completeness, we document the algorithm from the naive linear matching of RDMs in Sec. LABEL:app:linear-alg of the SI.

The above is a convenient way to impose the constraint on quantum computers, but it is computationally costly as the number of constraints in (16) increases exponentially as the number of overlapping sites $n$ on neighboring fragments. For each constraint equation, the expectation values $\langle\Sigma_{\alpha}\rangle$ has to be measured on the quantum computer, which therefore introduces an exponential overhead on the sampling complexity.

In the next section, we introduce a simple alternative to evaluate the mismatch between two RDMs on a quantum computer much faster based on a $\mathtt{SWAP}$ test.

The wave functions of two overlapping fragments are stored coherently as many amplitudes that suppose with each other. The beauty of quantum computers and algorithms lies at the ability to coherently manipulating such amplitudes simultaneously. We may naturally ask: are there quantum algorithms or circuits that can coherently achieve matching between an exponentially large number of amplitudes, without explicitly measuring each amplitude?

In quantum information, there is a class of quantum protocols to perform the task of estimating the overlap between two wave functions or RDMs under various assumptions ⁷⁶. Among these protocols, the $\mathtt{SWAP}$ test is widely used ^{77, 47}. Such a $\mathtt{SWAP}$ test on a quantum computer can also be naturally implemented by simple controlled-$\mathtt{SWAP}$ operations as in Fig. 4, showing a $\mathtt{SWAP}$ test between two qubits. The essence of a SWAP test is to entangle the symmetric and anti-symmetric subspaces of the two quantum states ($\ket{\phi}$ and $\ket{\psi}$) to a single ancillary qubit, such that the quantum state of the system before the final measurement is

By measuring the top single ancillary qubit in the usual computational $Z$-basis (collapsing it to either the $\ket{0}$ or $\ket{1}$ state), the overlap of the two qubit wave function, $|\langle\phi|\psi\rangle|$, can be directly obtained from the measurement outcome probability:

without requiring explicit estimation of the density matrix elements of each individual qubit.

Can we recast the linear matching conditions as linear combination of several $\mathtt{SWAP}$ tests? Observe that an equivalent condition alternative to Eq. (16) is the following quadratic matching condition

Interestingly, the above quadratic constraint can be rewritten as a linear combination of three different multi-qubit generalization of the $\mathtt{SWAP}$ tests (with each repeated multiple times), regardless of the number of overlapping sites (Fig. 1iiiq). Two of the $\mathtt{SWAP}$ tests are to estimate the purity of $\rho^{(A)}_{R}$ and $\rho^{(B)}_{R}$ each, while the third one is to estimate the overlap between $\rho^{(A)}_{R}$ and $\rho^{(B)}_{R}$. See SI Sec. LABEL:app:constraint-equivalent-proof for a proof of the equivalence between the two quantum matching conditions) and Sec. LABEL:app:swap-test on how to generalize the $\mathtt{SWAP}$ test on two qubits to a multi-qubit setting and how to relate the $\mathtt{SWAP}$ test results to the quadratic constraint.

The reformulation of the quadratic constraint allows us to estimate the mismatch between two fragments by measuring only a single ancilla qubit (estimating three different amplitudes). As compared to the linear constraint case where an exponentially large number of constraints have to be estimated individually ($4^{m}-1$ where $m=|R|$ is the number of overlapping sites again), the quadratic matching based on $\mathtt{SWAP}$ tests achieves an exponential saving in the types of measurements required.

Furthermore, the reduction of the mismatch to the estimation of only a few (three) amplitudes in $\mathtt{SWAP}$ tests allows an additional quadratic speedup by amplifying the amplitude of the ancilla qubit before measure it. We will discuss more details on how to achieve the quadratic speedup in Sec. 4.3. Admittedly, such amplitude amplification algorithm may be applied even to the naive linear RDM matching by boosting individual RDM amplitude, but the resulting quantum circuit will be much more complicated.

With an efficient way to estimate the quadratic penalty constraint established in Eq. (20), it now appears feasible to use this new constraint in Eq. (9) as in the case of linear constraint. However, the nature of the quadratic matching in Eq. (20) makes the same Lagrange multiplier optimization method used in the linear case invalid. We first discuss in more detail why this approach fails, in Sec. 3.4.1; we then describe an alternative way of treating the quadratic constraint as a penalty term to optimize the resulting objective function in Sec. 3.4.2.

We present a high-level framework of the main algorithm in this section. As a comparison, the QBE algorithm with naive linear matching can be found in SI Sec. LABEL:app:linear-alg. Code for the algorithms and data for generating the plots are available as open source on github ⁸².

To quantify the mismatch across all fragments, we define $\Delta\rho$ to be the root mean square density matrix mismatch averaged over all the overlapping sites of all the fragments according to

where $\Tr[\left(\rho^{(B)}_{r}-\rho^{(A)}_{r}\right)^{2}]=\mathcal{Q}_{quad}(\rho^{(A)}_{r};\rho^{(B)}_{r})$ as in Eq. (20), which may also be recognized as the Frobenius norm of $(\rho^{(B)}_{r}-\rho^{(A)}_{r})$. $N_{sites}$ is the total number of terms in the double sum in Eq. (24), $N_{sites}=\sum_{A\neq B}|\mathbb{E}^{(A)}\cap\mathbb{C}^{(B)}|$, with $|\mathcal{S}|$ denoting the number of elements in set $\mathcal{S}$.

The cost function $\mathcal{L}^{(A)}(\lambda)$ being optimized is discussed in Sec. 3.4.1. For clarity, we write it explicitly here

with $\mathcal{Q}_{quad}$ given by Eq. (20). We have omitted the term $\mathcal{E}^{(A)}$ for simplicity since the normalization of the wave function is guaranteed for a fault-tolerant quantum computer. However, this term can be important on a noisy quantum computer where the purity of the wave function can be contaminated. Note the expectation value in Eq. (25) has to be estimated by collecting samples on a quantum computer.

The quantum bootstrap embedding algorithm with quadratic penalty method is presented below in Alg. 1. The algorithm takes as its input the total Hamiltonian of the original system, and then perform the fragmentation and parameter initialization, followed by the main optimization loop to achieve the matching. Finally, it returns the optimized BE potential $V_{\rm BE}^{(A)}$ for any fragment $A$ and the final mismatch $\Delta\rho$. Inside the main loop (line 9 of Alg. 1), the cost function $\mathcal{L}^{(A)}(\lambda)$ for each fragment $A$ is minimized for a fixed penalty parameter $\lambda$ (line 10 and 11). The penalty $\lambda$ is then increased geometrically (line 12) until the mismatch criteria is met, i.e., $\Delta\rho\leq\varepsilon$.

A key step of the algorithm is the minimization of $\mathcal{L}^{(A)}(\lambda)$ at line 11, which consists of repeatedly generating the BE potential $V_{\rm BE}^{(A)}$ and estimate the mismatch using $\mathtt{SWAP}$ test. BE potentials $V_{\rm BE}^{(A)}$ are generated differently for different optimization algorithms. In our implementation, a quasi-Newton method, the L-BFGS-B ⁸³ algorithm, is used at line 11 for minimizing $\mathcal{L}^{(A)}(\lambda)$, where $V_{\rm BE}^{(A)}$ is proposed by the optimizer in order to estimate the inverse Hessian matrix to steer the optimization properly. Alternatively, if derivative-free methods such as Nelder-Mead ⁸⁴ is used, $V_{\rm BE}^{(A)}$ will be generated in a high-dimensional simplex defined by the coefficients $\{v_{\alpha}\}$ in Eq. (22), which is repeatedly refined.

Once $V_{\rm BE}^{(A)}$ is generated, the first term in the cost function in Eq. (25) is estimated by invoking the quantum eigensolver for the Hamiltonian $\left(H^{(A)}+V_{\rm BE}^{(A)}\right)$. The second term, the mismatch in Eq. (25) can be estimated by measurement outcomes of the ancilla qubit in the $\mathtt{SWAP}$ test (Sec. 3). The mismatch estimation at line 13 is performed in the same way as those in line 11. Note that the number of samples $N_{samp}^{\mathtt{SWAP}}$ (Eq. (27)) for the $\mathtt{SWAP}$ test estimation can be changed adaptively in different BE iterations for different accuracy, which we discuss in detail in the next section.

Two major quantum eigensolvers, QPE ⁸⁵ and VQE ³⁴ can be used in line 11 and 14 of Alg. 1 to estimate the cost function. QPE is an exact eigensolver, where the system wave function collapses to the exact ground state regardless of the number of evaluation qubits used. In contrast to QPE, VQE is an approximate eigensolver and the results depends on the choice of ansatz and the optimization algorithm used.

A crucial feature of a quantum eigensolver is its probabilistic nature, in a sense that any measurement collapses the entire quantum state. This perspective allows us to treat a quantum eigensolver as a sign-problem-free sampling oracle for correlated electronic structure problems where Ref. ⁸⁶ provides a concrete example.

The stochastic nature also means a more careful treatment on the number of samples is required to fully quantify any potential quantum speedup. In general, for typical iterative mixed quantum-classical algorithms, some parameters are usually passed from one iteration to the next, where the parameters are estimated by repeatedly sampling from a quantum eigensolver oracle through proper measurement. This means the uncertainty on these parameters estimated from one iteration has to be small enough to avoid a divergence of the algorithm as iteration continues.

In particular in the bootstrap embedding case, the sampling accuracy on the fragment overlap of each iteration has to be good enough such that the uncertainty of the mismatch passed to the next iteration will not spoil the iteration and lead to diverging results as iterations continue. In the following, sampling complexities of classical matching and $\mathtt{SWAP}$-test-based quantum matching are compared.

When estimating the overlap $S$ to an accuracy $\epsilon$ naively by density matrix tomography (TMG) of individual RDM elements, it is shown under mild assumptions that the total number of samples required (Sec. LABEL:app:sample-advantage-lin-quad of SI)

where $n$ is the number of qubits on the overlapping region, and $D$ is a system-dependent constant as a function of the two RDMs. In contrast, the quantum matching based on $\mathtt{SWAP}$ test costs

which is independent of the size $n$ of the overlapping region of two fragments. This demonstrates that our quadratic quantum matching achieves an exponential speedup compared to naive tomography of density matrices. This dramatic speedup is perhaps not that surprising because we only care about one particular observable (the overlap) instead of the full subsystem RDMs. Therefore, if the observable can be mapped to measurement outcome of few qubits by some quantum operations ($\mathtt{SWAP}$ test in this case), advantages are expected in general.

Moreover, the dependence of $N_{samp}^{\mathtt{SWAP}}(S,\epsilon)$ on the overlap $S$ and estimation accuracy $\epsilon$ allows an adaptive sampling schedule to be implemented for line 11 and 14 of Alg. 1. For example, we may use the overlap $S$ estimated from the previous BE iteration to compute the required $N_{samp}^{\mathtt{SWAP}}$ in the current BE iteration. The accuracy $\epsilon$ can also be dynamically tuned according to the error of the first term in Eq. (25), as well as the value of the penalty parameter $\lambda$. For example, at the beginning BE iterations, the mismatch ($\Delta\rho$ or more precisely $\mathcal{Q}_{quad}(\rho^{(A)}_{r};\rho^{(B)}_{r})$) is large so that a moderate $\epsilon$ suffices. As the BE iteration proceeds, the overlap converges exponentially, therefore an exponentially decreasing $\epsilon$ has to be used as well. A numerical value of $\epsilon$ needs be determined from case to case.

In addition, Eq. (27) suggests an interesting behavior. As the QBE algorithm proceeds and the overlap $S$ increases, fewer samples are needed to achieve a target accuracy. If $S$ approaches 1 exponentially fast as $S\sim 1-e^{-\gamma\cdot n_{\text{iter}}}$ for some constant $\gamma$, then the required number of samples for $\mathtt{SWAP}$ will decrease exponentially as BE iteration $n_{\rm iter}$ goes $N_{samp}^{\mathtt{SWAP}}\sim e^{-\gamma\cdot n_{\text{iter}}}/\epsilon^{2}$. In practice, the overlap of two subsystem can never approach 1 but saturates to a constant $0<c<1$ when matching is achieved, and therefore $N_{samp}^{\mathtt{SWAP}}\sim(1-c)/\epsilon^{2}$ still obeys the $1/\epsilon^{2}$ scaling generally. This, on the other hand, suggests that a larger overlapping region is advantageous to reduce $N_{samp}^{\mathtt{SWAP}}$ because the RDM of a larger subsystem of a pure state will have greater purity (hence larger $c$) in general.

The core of many quantum speedups over classical algorithms lie at the ability of quantum computers to directly manipulate the probability amplitude instead of probability itself, while classical computers only have access to probability. With this idea, the above perspective of treating a quantum eigensolver as an oracle where some amplitude is estimated through proper measurements allows us to achieve an additional quadratic speedup in our quantum bootstrap embedding algorithm. This section compares two different versions of quantum matching algorithms in QBE, the $\mathtt{SWAP}$ and the $\mathtt{SWAP}+$AE algorithms. However, the same argument of quadratic speedup applies to classical sampling based eigensolvers such as VMC as discussed in detail at the end of this section.

The intuition is that instead of directly measuring a small quantum amplitude to accumulate enough counts to reduce the error bar, we may use quantum algorithms to first amplify the amplitude before the measurement. One way of understanding this is that Eq. (27) contains an overlap-dependent prefactor $(1-S^{2})$ as discussed above. If the overlap $S$ (as a probability amplitude) can be manipulated on the quantum computer easily such that $(1-S^{2})$ is on the order of $\epsilon$, then $N_{samp}^{\mathtt{SWAP}}$ will be proportional to only $1/\epsilon$ instead of $1/\epsilon^{2}$. There are well-established ways of performing such amplitude amplification task via coherent quantum algorithms ⁴⁸. See SI Sec. LABEL:app:aa+binary-search for the construction of the amplitude amplification and binary search quantum algorithm.

In particular, in each iteration of the algorithm, it can be shown that by combining oblivious amplitude amplification and a binary search protocol, estimating the overlap up to precision $\epsilon$ between adjacent fragments takes $N_{samp}^{\mathtt{SWAP}+AE}$ samples (state preparation and $\mathtt{SWAP}$ tests)

regardless of the overlap $S$.

Comparing (28) with (27), the above analysis suggests that our coherent quantum matching algorithm achieves a quadratic speed up (up to a factor of ${\rm polylog}(\frac{1}{\epsilon})$) as compared to the $\mathtt{SWAP}$ test based quantum matching algorithm, which is consistent with typical behavior of a Grover-type of search algorithm. Moreover, in contrast to (26), an exponential advantage is present with respect to the size of the overlapping region, indicating the benefit of using our quadratic QBE algorithm for fragment matching in the presence of large overlapping region.

In the above, we leverage amplitude amplification to achieve a quadratic speedup of a quantum subroutine based on $\mathtt{SWAP}$ test. More generally, such amplitude amplification technique can be utilized to achieve a general quadratic speedup in the required number of samples for any Monte Carlo classical algorithms ^{87, 88, 89}. This can be understood by realizing that classical probability distributions may be encoded in the amplitudes of the quantum state of a quantum computer, where measurements performed after some unitary quantum computation is similar to sample from the quantum computer to extract the probability distributions. When treating the unitary quantum computation part as a quantum blackbox, it is then easier to understand the quadratic speedup in the number of samples as compared to classical Monte Carlo methods. In our case, the quantum blackbox is the quantum eigensolver used to find the ground state for each fragment, while the classical blackbox is the stochastic classical eigensolvers such as VMC.

We focus on demonstrating the convergence of QBE in the infinite sampling limit by using exact deterministic solver with the quadratic constraint in Eq. (20) and linear constraint in Eq. (16). As a standard benchmark system for electronic structure, we perform QBE on a H₈ chain under a minimal STO-3G basis, which is fragmented into six overlapping fragments each with six embedding orbitals. Fig. 5a shows the exponential convergence of the density mismatch for an H₈ molecule in both linear and quadratic constraint cases. This convergence behavior of QBE matches the convergence of classical BE in Fig. 2 with exact classical solver (FCI), demonstrating the correctness of the new constraints. The agreement on the convergence with classical BE in Fig. 2 is expected since at infinite sampling limit, the outer iterations in both classical and quantum BE are the same classical optimization routine.

To quantify how much energy error the final converged result has, Fig. 5b shows the absolute value of the error in energy using the energy in the last (11^th) iteration as a reference. We can see that the energy errors from both the linear and quadratic constraint algorithm exhibit similar exponential convergence as the density mismatch. Moreover, the energy in both cases converge to the same value within 10^-6 in the last iteration (not shown in the figure). We note that the linear constraint case shows a slightly oscillatory convergence, while the quadratic case is free of such oscillatory behavior. The fact that quadratic appears to converge slightly faster than linear may be coincidence for the system investigated, and the convergence rate in general depends on the optimization algorithm chosen. See Sec. LABEL:app:qbe-calculation of the SI for a detailed description on definition of the energy.

In the previous section, we have seen that our quantum bootstrap embedding algorithm convergence as expected in the infinite sampling limit. It is also seen (in the SI) that the approximate VQE leads to biased behavior on the density matching. In practice, only a finite number of samples can be collected on a quantum computer, and we will focus on theis scenario in this section. In particular, we present numerical data demonstrating the sampling advantage of our coherent quantum matching algorithm. Sec. 5.2.1 discusses the sampling advantage of the quantum matching algorithm for an overlapping region of increasing size, echoing the analytical sampling complexity derived in Sec. 4.2. In Sec. 5.2.2, the additional quadratic speedup in estimating the overlap via amplitude amplification and binary search (AE) is presented, which agrees with the theoretical sampling complexity in Sec. 4.3.