{'fraction_non_alphanumeric': 0.0383459937341745, 'fraction_numerical': 0.012681859147676066, 'mean_word_length': 4.312699498404013, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Doubling the size of quantum simulators by entanglement forging

Quantum computers are promising for simulations of chemical and physical systems, but the limited capabilities of today's quantum processors permit only small, and often approximate, simulations. Here we present a method, classical entanglement forging, that harnesses classical resources to capture quantum correlations and double the size of the system that can be simulated on quantum hardware. Shifting some of the computation to classical post-processing allows us to represent ten spin-orbitals on five qubits of an IBM Quantum processor to compute the ground state energy of the water molecule in the most accurate simulation to date. We discuss conditions for applicability of classical entanglement forging and present a roadmap for scaling to larger problems.

Quantum computers are promising for simulations of chemical and physical systems, but the limited capabilities of today's quantum processors permit only small, and often approximate, simulations. Here we present a method, classical entanglement forging, that harnesses classical resources to capture quantum correlations and double the size of the system that can be simulated on quantum hardware. Shifting some of the computation to classical post-processing allows us to represent ten spin-orbitals on five qubits of an IBM Quantum processor to compute the ground state energy of the water molecule in the most accurate simulation to date. We discuss conditions for applicability of classical entanglement forging and present a roadmap for scaling to larger problems.
Simulating quantum systems is an especially hard task for classical computers, making the realization of quantum computers potentially revolutionary for the study of chemistry, materials science, and fundamental physics. However, techniques like quantum phase estimation, which promises accurate chemical simulations, require hardware well beyond the present state of the art. While hardware capabilities continue to steadily advance, limitations on both quantity and quality of qubits are giving rise to a new family of algorithms that leverage additional classical resources to enable quantum computations requiring more qubits than physically available [1][2][3][4][5][6][7][8][9][10][11]. Such algorithms broadly benefit from a strategy of partitioning a problem into weakly interacting clusters, then correlating the results of each on a classical computer. A range of important systems naturally possess suitable partitions, including low-energy eigenstates of chemical [12] and lattice-model [13][14][15] Hamiltonians, systems embedded in a quantum bath [16,17], and static correlations associated with chemical bond-breaking processes [18][19][20]. Such quantum-classical programming requires a profound understanding of the tradeoff between quantum and classical computational resources [21] but is of theoretical and practical importance as it can extend the possibilities of classical and quantum computers alike.
Here we introduce a scheme we call classically forged entanglement, which represents a 2N -qubit wavefunction as multiple N -qubit states embedded in a classical computation. Beyond the reduction in requisite qubit number, offloading entanglement synthesis to classical processing permits the constituent N -qubit quantum circuits to be shallower, relaxing requirements on gate error and connectivity, at the cost of an increased number of circuit executions. In this report, we first theoretically describe entanglement forging, including scalable application to any state that can be partitioned into two weakly entangled halves. We then demonstrate entanglement forging using the variational quantum eigensolver (VQE) algorithm [22] to simulate a model of the water molecule [23], representing 10 spin-orbitals with only five qubits. Con-ducted entirely via the IBM Quantum cloud-computing service, the experiment benefits from state-of-the-art capabilities [24] such as active qubit reset, zero-noise extrapolation [25], chip-level parallel execution, and a twoqubit "hop" gate tailored for chemical simulation.
We begin with Schmidt decomposition (Fig. 1A), a standard application of singular value decomposition (SVD) that allows one to write any state |ψ of a bipartite N + N qubit system as Here |b n are the N -qubit bitstring states, also known as computational basis states, U, V are unitary operators respective to the two subsystems, and the Schmidt coefficients λ n may be taken to be non-negative. The flatter the distribution of Schmidt coefficients, the stronger the entanglement; a uniform distribution λ n = 1/ √ 2 N indicates the two halves of the system are maximally entangled, while only one nonzero coefficient corresponds to no entanglement.
Mixing classical and quantum information, entanglement-forging is more naturally expressed using density operators rather than wavefunctions, even for pure states such as |ψ . As shown in the supplementary material ( SM.1), one can write the density operator as where we have used the definition |φ p xy = |x + i p |y / √ 2 with p ∈ {0, 1, 2, 3} = Z 4 . For example, in the minimal case of two qubits (Fig. 1B), |φ p 01 correspond to four equatorial points on Bloch sphere, rewriting a quantum superposition of product states in terms of classical products of superposition states. Eq. (2) generalizes arXiv:2104.10220v1 [quant-ph] 20 Apr 2021 methods proposed in [1,26,27], and is connected to tensor network representations of quantum circuits [8,28], variational simulation of open quantum systems [29], and the encoding of open-shell singlet and triplet states [30]. The expectation of a 2N -qubit operator , and each constituent requires only N qubits to evaluate. The resulting summation for O is not obviously scalable, involving as many as 2 N +1 (2 N +1 − 1) distinct Nqubit quantum circuits. Nonetheless, if one restricts to simulations of sufficiently weak entanglement, O can be efficiently estimated ( SM.5) by sampling each circuit in proportion to the associated coefficients λ n λ m in (3), with a total number of samples for target precision scaling as Executing a quantum circuit once provides one sample of the corresponding expectation value, such that total runtime scales linearly with S. Since the one-norm decreases toward 1 in the limit of weak entanglement, the overhead cost of entanglement forging is smaller for simulations of states divisible into weakly-entangled halves, such as the spin-up and spin-down components of some molecular ground states (Fig. 1D), and scales efficiently when the one-norm is at most polynomial in the problem size. For example, in some statically correlated ground states, S can be independent of the number of basis orbitals. Outside of the domain of scalability, entanglement forging still enables useful heuristic simulations beyond the standard capacity of given quantum hardware, which may be realized with precision by truncating the list of bitstring states retained in the Schmidt decomposition.
Alternatively, this overhead may be reduced to a constant factor independent of qubit number via a complementary scheme ( SM.2) simulating quantum correlations between subsystems using those within a subsystem, rather than using classical correlations as above. This method can be seen as an application of forging in the Heisenberg picture, reinterpreting an observable acting on N + N qubits as a classical mixture of operators describing the forward and backward time evolution of N -qubits, at a cost of deeper circuits. Provided certain sampling assumptions, this method is not limited to weakly entangled states, so it may be applicable to a wider range of systems.
We use entanglement forging for a VQE simulation of the water molecule, as schematized in Figure 2. Starting from the minimal STO-6G basis, we freeze the core oxy- Experimental ansatz structure. A, The five active molecular orbitals are encoded in a line of qubits via the Jordan-Wigner mapping. The qubits are prepared in either a bitstring state, or a superposition of two bitstrings (shown), and then acted upon by parameterized hop gates (B). In order to evaluate the Pauli strings of the Hamiltonian, each state is prepared many times, and rotated appropriately prior to measurement in the computational basis. C, Parameterized molecular geometry from which molecular orbitals are determined. D, Device map of ibmq dublin, with a highlighted example of how two geometries may be solved in parallel.
gen 1s and the out-of-plane oxygen 2p orbitals, leaving an active space of 10 spin-orbitals. The Jordan Wigner mapping encodes either the spin-up or spin-down orbitals onto each five-qubit quantum circuit, and we simplify the problem structure by asserting V = U per the known symmetry between spin-polarizations in the closed-shell singlet ground state. How best to construct an ansatz circuit for VQE remains an open and active research topic [31]. To facilitate our demonstration of entanglement forging, we first ran classical simulations of VQE without entanglement forging to obtain a 10-qubit circuit that performs well for this problem. Based on this circuit, we selected the quantum gates in U , and truncated the Schmidt decomposition to k = 3 of the 10 possible bitstrings, namely |b 1 = |11100 , |b 2 = |01110 , and |b 3 = |01101 . Such truncation was convenient due to a technical limitation on the number of unique circuits, distinct from the more fundamental limitation of overhead time scaling with the number of samples S discussed above; future implementations may variationally explore larger subsets of possible bitstrings as discussed in SM.5. The resulting ansatz is thus tailored to our particular problem and parameterized both by two qubit gate rotation angles and by the Schmidt coefficients, discussed further below. To realize each circuit required for equation (3), we prepare either one bitstring, or a superposition of two. For example, Fig. 2A shows how to prepare |φ 0 b1b3 , and a general construction is given in SM. 3.
Each initialized state is then acted upon by the unitary U , here comprised of five two-qubit "hop gates" ( Fig.  2A,B), where each acts according to the matrix In the simulated molecule, the rotation by ϕ moves particles between orbitals, while the entry -1 provides a CPHASE-like two-particle interaction. A product of hop gates is universal with respect to real-valued wavefunc-tions of fixed particle number, making the gate appealing for simulating eigenstates of time-reversal invariant Hamiltonians. Details of compilation on quantum hardware appear in SM.4. The runtime cost of entanglement forging makes fast circuit execution critical. We implemented the VQE routine using the Qiskit quantum computing package, and executed the experiment via the IBM Quantum cloud computing service on the ibmq dublin processor [32,33]. Active qubit reset enabled execution at a 10 kHz repetition rate. Noting that distinct molecular geometries present embarrassingly parallel problems, we mapped two independent problems onto separate five-qubit chains (Fig. 2D), doubling throughput to 20,000 five-qubit circuits per second ( SM.8). While freedom to independently specify the number of samples for each unique circuit would best economize runtime, here we approximated this behavior by submitting multiple copies of circuits in proportion to the desired weighting ( SM.6).
Properties of the particular problem and ansatz permitted many circuits to be omitted, further accelerating execution. Generally, superposition states |φ p bnbm may be skipped when evaluating Hamiltonian terms acting on only one of the two subsystems. Likewise, since here the operators O 1 and O 2 in the Hamiltonian are real valued in the bitstring basis, contributions from |φ p bnbm with odd p must vanish. Finally, observing that our chosen hop gates nominally do not modify the Hartree-Fock (HF) state, we omit the corresponding circuits and use instead the classical mean-field result, fixing b 1 |Õ|b 1 = E HF . This simplification cannot capture excitations of unpaired electrons across the HF occupied-virtual threshold, and risks unphysical results as hardware errors impact intermediate calculations unevenly. Nonetheless, unburdening the quantum processor of the dominant, classically-accessible part of the problem reduced S by roughly an order of magnitude for our near-equilibrium experiments, while retaining acceptable overall accuracy. Notably, this acceleration was enabled by the mixture of classical and quantum elements in the entanglement-forging representation.
VQE was repeated to study the behavior of entan-

R R
Symmetric dissociation Single-bond cleavage Angular pro le In the upper plots, VQE results using the entanglement-forging ansatz on ibmq dublin appear as filled blue circles, alongside curves indicating the classically computed Hartree Fock (HF) and full configuration interaction (FCI) values for the active space. Unfilled shapes indicate results from running VQE on a noiseless classical simulator using the same 3-bitstring ansatz used on ibmq dublin (k = 3), and using a larger 6-bitstring ansatz (k = 6). Lower plots indicate absolute errors, |E − EFCI|. Error bars are one standard deviation produced by bootstrapping measurement histograms, representing precision but not systematic error or drift.

A B
Symmetric dissociation Single-bond cleavage Angular pro le glement forging across a variety of molecular geometries. For each geometry, optimization was performed in approximately 100 iterations of an SPSA (simultaneous perturbation stochastic approximation) algorithm [34,35]. Figure  Where not otherwise specified, nominal equilibrium values of R eq = 0.958 A and θ eq = 104.478 • were used [36]. Deviations between the VQE results and exact full configuration interaction (FCI) energies trend from ∼ 1 − 10 mHa to ∼ 10 − 100 mHa as bond lengths are increased, in line with the expectation that entanglement forging should work best for problems with weakly entangled ground states. We repeated all experiments on a noiseless classical simulator, first verifying that this trend persists in the absence of gate errors (orange circles), and second observing that accuracy in the stretched regime improves substantially upon increasing the number of represented bitstring states k from three to six (purple squares), enabling the forging to describe stronger entanglement at the expense of sampling more distinct circuits. Finally, to emphasize that entanglement forging can extend to arbitrary accuracy with no increase in the number of qubits, we ran a noiseless VQE simulation at the nominal equilibrium geometry including all 10 three-occupation bitstrings along with more hop gates, which converged within 1.6 mHa of the FCI energy ( SM.11). The Schmidt coefficients λ n (Fig. 1B) are among the parameters optimized by VQE with entanglement forging ( SM.10). These coefficients are uniquely defined by the ground-state wavefunction and choice of partition, and thus, assuming adequate convergence of VQE to the ground state, do not depend on the selection of gates in U and V up to reordering of the index n. Figure 4 shows the Schmidt coefficients obtained on ibmq dublin, along with those obtained on a noiseless classical simulator with the same 3-bitstring ansatz (k = 3). The plots verify that weaker entanglement correlates with better performance in Fig. 3, and clarify how bond stretching strengthens spin-polarization entanglement. For instance, removing a single hydrogen (Fig. 4B) leads to a wavefunction dominated by two bitstrings, consistent with a picture of dissociation into spin-1/2 fragments with opposite polariza-tions. Further discussion of Schmidt coefficient distributions supported by FCI calculations is provided in SM.9.
We have demonstrated entanglement forging by computing the ground state energy of a water molecule for varying geometries using five qubits to represent ten spin-orbitals, the most accurate VQE simulation of this molecule using quantum hardware to date. Entanglement forging opens a vista of experimental possibilities, such as dividing a lattice system using a spatialrather than spin-partition; introducing an adaptive bisection based on orbital optimization; using entanglement forging in simulation algorithms beyond VQE, including error-corrected techniques such as quantum phase estimation; and more generally for any experiment otherwise inaccessible due to limited qubit number or connectivity where the overhead of forging may be tolerated. Where overhead remains prohibitive, Heisenberg-picture forging ( SM.2) may provide an alternative path to simulation of large systems. Additional exploration of techniques incorporating classical resources with quantum processing may further enhance the power of near term quantum simulations and help realize quantum advantage for practical applications. modeling and analyses. All authors contributed to the manuscript. Competing interests: Elements of this work are included in a patent filed by the International Business Machines Corporation with the US Patent and Trademark Office. Data and materials availability: Experiment data and code available upon reasonable request to the authors.

SM.1. Decomposition to superpositions of N -qubit bitstrings
We would like to forge a 2N -qubit state |Ψ = n λ n |b n ⊗ |b n where λ n are real coefficients. For each distinct pair of bitstrings x, y and each p ∈ {0, 1, 2, 3} ≡ Z 4 define a state One gets Applying the above identity to each off-diagonal term in |Ψ Ψ| we arrive at

SM.2. Entanglement forging in the Heisenberg picture
First let us define a class of quantum states whose entanglement can be efficiently forged. Consider a system of 2N qubits and states of the form where the tensor product separates two N -qubit registers, λ x are real coefficients, and U is a unitary operator with real matrix elements in the standard basis. The coefficients λ x have to be normalized such that x λ 2 x = 1. We assume that U admits an efficient implementation by a quantum circuit. It is natural to consider quantum circuits composed of real gates such as the Hadamard, Y -rotation, CNOT, or a hop gate. This ensures that U has real matrix elements. The vector of 2 N coefficients λ x may be specified either explicitly (for small N ) or implicitly by a classical algorithm that can efficiently perform the following tasks: (i) sample a bit string x from the probability distribution λ 2 x ; (ii) compute the ratio λ y /λ x for a given pair of bit strings x, y. Such implicit description of exponentially large vectors is commonly used in Quantum Monte Carlo simulations based on tensor networks [37] or neural network states [38]. It follows the framework of computationally tractable quantum states introduced in [39].
Given Pauli observables O 1 , O 2 ∈ {I, X, Y, Z} ⊗N , our goal is to estimate the expected value We show how to accomplish this task by a series of N -qubit experiments. A typical experiment prepares a state U † CU |x for a suitable N -qubit Clifford circuit C and a bit string x. Then every qubit of the state U † CU |x is measured in the standard basis. Note that the circuit U † can be obtained from U by inverting the order of gates and replacing each gate by its inverse. The size of the Clifford circuit C depends on the form of the observables O 1 , O 2 . Let |O j | be the Hamming weight of O j , that is, the number of single-qubit terms X, Y , Z that appear in O j . Then C contains at most 2(|O 1 | + |O 2 |) CNOT gates and some single-qubit gates. The number of experiments required to approximate the expected value Eq. (7) with a precision is proportional to 1/ 2 with a constant prefactor. Crucially, the number of experiments does not depend on the number of qubits N or the amount of entanglement in the forged state |ψ .
The key idea behind this method is to convert a Pauli observable O 1 ⊗ O 2 describing a bipartite system of N + N qubits into a classical mixture of observables C * ⊗ C, where C is a self-adjoint N -qubit Clifford operator and C * is the complex conjugate of C. By interpreting C and C * as operators describing the forward and the backward time evolution of the same N -qubit register we express the expected value µ as a mixture of quantum probabilities that can be measured on an N -qubit device. Since the analysis of the method is performed in the Heisenberg picture at the level of observables, we shall refer to it as the Heisenberg forging. Meanwhile the forging method described in the main text works primarily in the Schrödinger picture by decomposing an entangled state of N + N qubits into a classical mixture of product states. We shall refer to the latter as the Schrödinger forging, to avoid a confusion between the two methods.
The Heisenberg forging roughly doubles the size of quantum circuits that need to be executed compared with the Schrödinger forging (the former needs to implement both U and U † while the latter only needs to implement U ). However, the Heisenberg forging requires shorter circuits compared with the preparation of the full state |ψ on a hardware with 2N qubits. Indeed, a quantum circuit preparing the state |ψ includes two copies of U as well as some extra gates generating the initial entanglement between the N -qubit registers. Thus we expect that the Heisenberg forging reduces both the number of qubits and the number of gates per experiment while increasing the number of experiments only by a constant factor.
Using states Eq. (6) as a variational ansatz in quantum simulations is justified if the exact ground state has real amplitudes in the standard basis and is invariant under exchanging the two N -qubit registers. As discussed in the main text, this is the case for molecular electronic structure Hamiltonians assuming that the N -qubit registers represent spin-up and spin-down orbitals. Specializing Eq. (6) to the case U = I gives variational states in which every spatial orbital is either empty or occupied by a pair of electrons in the singlet state. Such states play the central role in the Restricted Hartree-Fock space formalism [40]. As shown in Ref. [19], minimizing the energy of a molecular Hamiltonian over variational states of the form Eq. (6) with U = I is equivalent to minimizing the energy of a quantum spin Hamiltonian with two-spin Heisenberg-type interactions and an external magnetic field over states |λ = x λ x |x , where λ x are the Schmidt coefficients that appear in Eq. (6). We anticipate that the ansatz Eq. (6) with U = I may improve upon the Restricted Hartree-Fock formulation since it is based on more general variational states. We leave applications of the Heisenberg forging for a future work.
In the rest of this section we assume that O 1 and O 2 are non-identity Pauli operators, which is the most interesting case. Otherwise we have O 1 = I or O 2 = I. For concreteness, assume O 2 = I. Then the expected value Eq. (7) becomes µ = x λ 2 x x|U † O 1 U |x . One can estimate µ with a precision by generating M ∼ −2 samples x 1 , . . . , x M ∈ {0, 1} N from the distribution λ 2 x , measuring the eigenvalue of O 1 on each state U |x j , and computing the sample mean of the measured eigenvalues. Let where a j are real coefficients such that |a j | ≤ 1, and C j are N -qubit Clifford operators.
We shall now convert the expected value ψ|C * j ⊗ C j |ψ into a mixture of quantum probabilities that can be measured on an N -qubit device. From Eq. (6) one gets Here we used the identity which follows from the assumption that U has real matrix elements. Combining Eqs. (9,10) one gets where The first term in Eq. (11) can be estimated on an N -qubit device by sampling a bit string x from the distribution λ 2 x , preparing a state U |x , and measuring the eigenvalue of O 1 O 2 (note that the first term in Eq. (11) . Below we focus on the second term in Eq. (11) and show how to estimate the quantities µ j . Define a function By assumption, one can efficiently compute R(x, y) for a given pair x, y. For each bit string x define a conditional probability distribution One can sample a bit string y from P j (y|x) on a quantum device with N qubits by preparing the state U † C j U |x and measuring every qubit in the standard basis. Furthermore, Eq. (12) implies that µ j is the expected value of R(x, y) over the probability distribution λ 2 x P j (y|x), that is, The random variable R(x, y) has the variance at most one since Thus one can estimate µ j with an additive error by generating roughly 1/ 2 samples (x, y) from the distribution λ 2 x P j (y|x) and computing the sample mean of R(x, y). By estimating each term µ j in Eq. (11) with a precision /4 and estimating the expected value Re( x|U † O 1 O 2 U |x ) with a precision /2 gives the desired -estimate of µ (use the triangle inequality and the assumption |a j | ≤ 1). The total number of N -qubit experiments required to obtain these estimates is proportional to 1/ 2 with a constant prefactor.
In the rest of this section we explicitly construct the decomposition Eq. (8) for non-identity Pauli observables O 1 , O 2 ∈ {I, X, Y, Z} ⊗N . Let w j be the number of single-qubit Y terms that appear in O j . Note that O * j = (−1) wj O j , that is, the matrix of O j in the standard basis is real if w j is even and imaginary if w j is odd. Since the state |ψ is real, one has µ = ψ|O 1 ⊗ O 2 |ψ = 0 unless w 1 and w 2 have the same parity. Below we assume that this is the case, that is, O * j = σO j for σ = (−1) w1 = (−1) w2 . Recall that Pauli observables either commute or anti-commute.
The circuit V contains at most |O 1 | + |O 2 | − 2 CNOT gates and some single-qubit Clifford gates. Here |O j | is the Hamming weight O j .
Let V be the Clifford circuit from Lemma 1. Define Clifford operators Taking the partial transpose of Eq. (19) over the first N -qubit register, using the fact that O j and C j are self-adjoint, and the identity O * j = σO j one gets This is the desired decomposition Eq. (8) with a 1 = σ, a 2 = −σ, and a 0 = a 3 = a 4 = 0. Lemma 1 and Eq. (18) imply that C 1 and C 2 can be implemented by circuits with at most 2(|O 1 | + |O 2 |) CNOTs.
The circuit V contains at most |O 1 | + |O 2 | − 2 CNOT gates and some single-qubit Clifford gates.
Let V be the Clifford circuit from Lemma 2. Define N -qubit operators where α, β ∈ {0, 1} and CZ p,q denotes the controlled-Z gate acting on qubits p, q. Using the identity A simple algebra shows that Taking the partial transpose over the first N -qubit register, using the identity O * j = σO j , and the fact that C j and O j are self-adjoint one gets This is the desired decomposition Eq. (8) with C 1 = C 0,0 , C 2 = C 0,1 , C 3 = C 1,0 , C 4 = C 1,1 , a 0 = a 1 = a 4 = σ, and a 2 = a 3 = −σ. It remains to prove Lemmas 1,2.
Proof of Lemma 1. We shall convert O 1 and O 2 to single-qubit Pauli operators X q and Z q respectively by a sequence of steps where W is a Clifford circuit composed of CNOT, CZ, and single-qubit gates. We shall choose W such that at each step the combined weight |O 1 | + |O 2 | is reduced at least by number of two-qubit gates in W . The desired circuit V is then obtained as the composition of the circuits W applied at each step. The first step converts O 1 and O 2 into the standard form such that their action on any qubit falls into one of five cases shown below.
The corresponding circuit W can be easily constructed as a product of single-qubit Clifford gates. This gives rise to a partition of N qubits into five disjoint subsets, [N ] = ABCDE. For example, A contains all qubits j such that O 1 and O 2 act on j by Pauli X and Z respectively. Note that A has an odd size since otherwise O 1 and O 2 would commute. In particular, A = ∅.
Suppose B = ∅. Apply W = CZ a,b for some a ∈ A and b ∈ B. This reduces |O 1 | by one without changing O 2 . Suppose CD = ∅. Apply W = CNOT b,a for some a ∈ A and b ∈ CD. This reduces |O 2 | by one without changing O 1 .
One can easily check that W † O 1 W = X 1 and W † O 2 W = Z 1 . Thus W reduces the combined weight |O 1 | + |O 2 | by 4k. Furthermore, W contains 3k ≤ 4k two-qubit gates.
Proof of Lemma 2. We shall use the notations introduced in the proof of Lemma 1. Consider the standard form of O 1 and O 2 . Suppose A = ∅. Note that |A| is even since O 1 and O 2 commute. Assume wlog that A = {1, 2, . . . , 2k} for some integer k. Choose One can can easily check that W maps X 1 X 2 . . . X 2k and Z 1 Z 2 . . . Z 2k to Z 1 Z 3 . . . Z 2k−1 and Z 2 Z 4 · · · Z 2k respectively. The combined weight |O 1 | + |O 2 | is reduced by 2k and W contains 2k two-qubit gates. From now on we can assume A = ∅.
If B = ∅ and D = ∅ then apply W = CNOT d,b for some b ∈ B and d ∈ D. This reduces |O 1 | by one without changing O 2 .
If C = ∅ and D = ∅ then apply W = CNOT d,c for some c ∈ C and d ∈ D. This reduces |O 2 | by one without changing O 1 .
After a sequence of steps as above we have A, D = ∅. Note that B = ∅ and C = ∅ since we assumed that O 1 = I and O 2 = I.
If |B| ≥ 2 then apply W = CNOT b,b for some qubits b, b ∈ B. This reduces |O 1 | by one without changing O 2 . If |C| ≥ 2 then apply W = CNOT c,c for some qubits c, c ∈ C. This reduces |O 2 | by one without changing O 1 .
After a sequence of steps as above we have A, D = ∅, |B| = 1, and |C| = 1, that is, O 1 = Z p and O 2 = Z q for some pair of qubits p = q.

SM.3. State initialization routines
For the entanglement-forging scheme used in the experiment, we need to initialize N qubits in a superposition state |φ p xy = (|x + i p |y )/ √ 2, where x and y are N qubit bitstrings, and x = y. A general construction of |φ p xy proceeds as follows: 1. Find an index k where x k = y k .
2. If x k = 1 (y k = 0), swap the definitions of x and y, and substitute p → (−p mod 4), using the fact that |φ p xy is the same as |φ −p yx up to an inconsequential global phase. 4. Given a set of N qubits prepared in |0 ⊗N , apply the single-qubit gates ⊗ i∈T X i .
The above procedure may become difficult to execute with high fidelity in simulations of large systems as limited processor connectivity may lead to large numbers of costly swap operations, or simply as the Hamming distance between pairs of included bitstrings grows and necessitates more CNOTs in step 6. Moreover, circuits with different initialization gates will experience different errors, and these uneven error rates can lead to inference of non-physical average behaviors, e.g. energies below the correct ground state energy. To circumvent this issue, one may replace the circuits initialized as N -qubit superposition states by a larger number of circuits initialized as products of single-qubit states, which are readily prepared with high fidelity. Suppose x = y are N -bit strings that differ on d bits. For each qubit j and an integer p define a single-qubit state Define an N -qubit tensor product state |ψ p xy = |ψ p1 xy ⊗ |ψ p2 xy ⊗ · · · ⊗ |ψ pN xy .
Such state can be easily prepared starting from the basis vector |0 N by applying Hadamard gates, Pauli X gates, and single-qubit Z-rotations by the angle ±πp/2d. After simple algebra one gets a decomposition Note that the righthand side is a linear combination of 2N -qubit tensor product states. We can use this decomposition to classically forge entanglement for an arbitrary state Indeed, write the density matrix |ψ ψ| as a linear combination of diagonal terms |b n b n | ⊗2 and off-diagonal terms |b n b m | ⊗2 + |b m b n | ⊗2 with n = m. Applying Eq. (24) with x = b n and y = b m to each off-diagonal term one finally arrives at Thus the method is practical only if d n,m is sufficiently small for all bit strings b n and b m that contribute to |ψ . gate denoted by green circles. Panel b likewise defines a "modified hop gate," which does not involve a swap, and acts like the hop gate except that it leaves the |11 state unchanged. Following panel c, we start from the set of hop gates in Fig. 2, use the fact that in the experiment we fix ϕ 3 = 0 to reduce that hop gate to a CPHASE gate, and unpack gates 1, 4, and 5 according to the hop gate definition. (The CPHASE gate is compiled using Hadamards and a CNOT). We now compile hop gate 2 differently depending on the initialization subcircuits (purple rectangles) of the given circuit. Unpacking gate 2 per its definition and untwisting all crossed wires produces the circuit in the top right of panel d. The wires have been permuted to reveal the linear connectivity; in terms of the letters labeling the qubits in the figure, the required connectivity is C-B-A-D-E. To this set of gates we will need to prepend initialization subcircuits preparing superpositions of bitstrings. Each such subcircuit in this experiment contains one additional CNOT (red); to avoid any gate-based swaps, the two involved qubits must also be connected. The two purple subcircuits in d, |φ 0 b1b2 and |φ 0 b2b3 , respectively require connections A-D and D-E. Since this is compatible with the connectivity C-B-A-D-E, we compile hop gate 2 as per panel d whenever these states are prepared.

SM.4. Connectivity and swaps
In contrast, superposing the third pair of bitstrings by the |φ 0 b1b3 subcircuit in e requires A-E connectivity. To realize this, we note that this state preparation ideally never provides |11 as an input to hop gate 2, and thus hop gate 2 may be replaced by a modified hop gate with no nominal change in behavior. With the swap removed from hop gate 2, the desired linear connectivity of the hop-gate circuit becomes C-B-A-E-D, providing the A-E connection facilitating |φ 0 b1b3 . Nonphysical results may be inferred in cases where hardware noise acts unevenly across the set of executed circuits, as will occur to some extent when varying gate compilations depending on the state preparation. However, the total error in the circuit is reduced by the avoidance of swaps, and we accept the tradeoff for this experiment. An intriguing question is how generally entanglement forging and other circuit decomposition techniques may be used to circumvent limitations of device connectivity [9]. Future entanglement-forging experiments may reduce connectivity requirements using the decomposition to single-qubit stabilizer states presented in SM.3.

SM.5. Weighted sampling of circuits: Theory
Suppose O 1 , O 2 ∈ {I, X, Y, Z} ⊗N are N -qubit Pauli observables and |ψ is a 2N -qubit state. Entanglement forging provides a decomposition where ρ a are N -qubit pure states and µ a are real coefficients simply related to the Schmidt coefficients of |ψ . Let V i be a Clifford circuit that maps each Pauli term X or Y that appears in O i to the Pauli Z. Note that V i includes only single-qubit Clifford gates. Then where O i (x) takes value +1 or −1 if the parity of x restricted to the support of O i is even or odd respectively. We would like to estimate the quantity ψ|O 1 ⊗ O 2 |ψ on a quantum device with only N qubits by a series of experiments such that each experiment prepares a state V i ρ a V † i for some pair (i, a) and measures each qubit in the computational basis. We claim that estimating ψ|O 1 ⊗ O 2 |ψ with an additive error and a confidence level 99% requires at most S experiments, where For large problems, Monte Carlo methods may be used to generate the necessary M = S/2 state preparations ρ a . Indeed, define a normalized probability distribution π a = |µ a |/ µ 1 , where a = 1, 2, . . . , , and let a 1 , a 2 , . . . , a M ∈ {1, 2, . . . , } be M independent samples from the distribution. For each j = 1, 2, . . . , M perform two experiments: (1) prepare the state V 1 ρ aj V † 1 , measure every qubit in the Z-basis, and record the measured bit string x j ∈ {0, 1} N ; (2) prepare the state V 2 ρ aj V † 2 , measure every qubit in the Z-basis, and record the measured bit string y j ∈ {0, 1} N . Define a random variable One can easily check that the mean and the variance of f are Here the expectation values are taken over the random choice of a 1 , . . . , a M and the random measurement outcomes.
We note that estimating ψ|O 1 ⊗ O 2 |ψ on a 2N -qubit device with a precision and confidence level 99% would require at most 100 −2 experiments (each experiment prepares the state V 1 ⊗ V 2 |ψ and measures every qubit in the Z-basis). From Eq. (28) one infers that entanglement forging increases the required number of experiments roughly by the factor 2( µ 1 ) 2 . However, each experiment requires N instead of 2N qubits. In addition, as argued in the main text, quantum circuits preparing the states ρ a are much simpler compared with the circuit preparing the full state ψ.
Let λ = (λ 1 , λ 2 , . . . , λ 2 N ) be the Schmidt coefficients of |ψ . Eq. (2) from the main text gives Thus entanglement forging is mostly useful for weakly entangled states such that λ 1 = 2 N i=1 |λ i | is a constant or a slowly growing function of N (note that λ 1 = √ 2 N in the worst case when |ψ is a maximally entangled state of N + N qubits). In practice, one may wish to limit the maximum number of experiments S by some specified cutoff S max . Combining Eqs. (28,31) that determine the number of experiments S = S( , λ) one obtains a non-linear constraint S( , λ) ≤ S max on the vector of Schmidt coefficients λ. This constraint should be incorporated into the classical optimizer that minimizes the variational energy. SM.6. Weighted sampling of circuits: Implementation As the number of bitstrings included in the experimental demonstration (3 bitstrings) was computationally tractable, Monte Carlo methods were not required. To minimize the variance of our estimate of the sum in equation 3, we wish to allocate the total experimental samples S to the various circuits b n |Õ k |b m in proportion to the respective coefficients λ n λ m . However, at the time of the experiment Qiskit did not yet support non-uniform sampling of a list of distinct quantum circuits, instead sampling all J circuits in the submitted list, or job, an equal number of times s = S/J. We thus wrote a small "copysampling" function to approximate non-uniform sampling by populating the job with proportionally more copies of more important circuits. The function takes as inputs the desired normalized statistical weight w c for each circuit c, and the final size J of the job to be submitted. Thus a target expression for the number of copies desired for circuit c is w c J, though this is not necessarily integer-valued. For simplicity, the copysampling function first ensures the job includes at least one copy of every possible circuit. Next, it deterministically appends floor(w c J −1) additional copies of each circuit c. There are then J − c 1+floor(w c J −1) remaining spots in the job. These are allocated to circuits randomly according to the remaining residuals, without replacement, to approximate the desired weighting distribution. A similar partially-deterministic shot allocation scheme was studied in [41]. The job is then executed, and the results ("counts") for all copies of c are merged before analysis.
The above is limited in that the weighted-sampling is coarse-grained by the number of samples s per copy of each circuit. The desired weighted-sampling distribution could be better approximated by setting s = 1 and J = S. However, at present there is appreciable overhead in execution time scaling with the job size J. In our experiment, we found J ≈ 800 and s ≈ 2000 gave an acceptable compromise between job-size overhead time and efficiency of weighted sampling.
Our sampling weights accounted for both the coefficients λ n λ m and the zero-noise extrapolation stretch factors discussed below. We note that possible extensions could further optimize sampling weights for each quantum circuit by accounting for weights of the Pauli strings in the Hamiltonian, optionally further informed by the expectations of Pauli strings obtained in the previous VQE iteration. Optimally assigning samples accounting for Pauli strings is nontrivial as Hamiltonian terms contain pairwise products of N -qubit Pauli-string expectations, and moreover groups of compatible Paulis are in our experiment measured jointly via the tensor-product basis grouping method [42] standard in Qiskit. Nonetheless, preliminary simulations suggest such weighting schemes can further reduce the required number of samples S by more than ∼ 25%.

SM.7. Repeated-gate zero-noise extrapolation
The accuracy of energies computed on the quantum hardware was significantly improved by mitigating gate errors via zero-noise extrapolation, or ZNE [25,43]. To avoid the overhead of calibrating stretched gates, we implemented a repeated-gate ZNE routine similar to that in [44]. For each quantum circuit, a copy is made by applying the transformation G → GG −1 G to each primitive gate, approximating a noise-amplification factor of 3. We use firstorder extrapolation of the results from the original and copied circuits to infer the approximate zero-noise result for each distinct circuit (i.e. for each distinct [state-preparation, Pauli-string] combination) used to compute the desired observable. Additionally, because the extrapolated result depends more sensitively on the original circuit than on the noise-amplified circuit, to improve overall precision we assign a proportionally higher weight to the original circuit in the weighted sampling routine described above.

SM.8. Parallel circuit execution
To reduce the total data acquisition time, we executed pairs of VQE problems simultaneously on different subsets of the Dublin processor. Physical qubits were selected based on automated device benchmarking performed shortly before to the start of each VQE run. To reduce the risk of crosstalk between problems, we required that no direct connectivity exist between the two sets of physical qubits, with a buffer of at least one idle physical qubit. A similar parallelization scheme was recently characterized in [45]. Future parallelized experiments might suppress readout crosstalk as needed by choosing sets of qubits disjoint with respect to readout-line multiplexing, or by applying efficient mitigation methods such as in [46].
This parallelization routine was implemented in software to be a modular wrapper of the experiment code. Using the standard Python multiprocessing package, the wrapper first launches a Python process executing one copy of the experiment code for each independent VQE problem. In each iteration of the VQE optimization, each of these processes generates a list of circuits, or job, it needs executed. Rather than passing the job directly to the backend as usual, each process passes its job to a common "multiplexer" Python process. The multiplexer respectively merges FIG. 6. Accuracy of truncated Schmidt decompositions. The three panels correspond to the three sweeps of molecular geometries as in Figs. 3,4. For each choice of bitstring truncation k, we plot one minus the projection of the FCI wavefunction into the space defined by the leading k bitstrings. In other words, we plot the sum of squares of the FCI Schmidt coefficients excluded by the truncation. Values near 1 (0) indicate that the FCI wavefunction cannot (can) be accurately represented in the truncated space. Unlike in the rest of the text, these FCI calculations do not freeze the oxygen 2p orbital.
the circuits in the separate jobs into a single new job, sends that to the backend for execution, and then separates the component results and returns them to the corresponding VQE processes. This framework is convenient as the multiplexer can be swapped in as a virtual backend with little modification of the original experiment code.

SM.9. Schmidt coefficients and truncation
For technical convenience in the forging demonstration, rather than sampling all possible bitstrings, we truncated the space to include only the leading k = 3 bitstrings, that is, set λ n>k = 0. Noiseless simulations with k = 3 and k = 6 were performed for comparison (Fig. 3). Figure SM.9 illustrates how the choice of k limits accuracy of describing the various ground state wavefunctions, here obtained via FCI calculations. For the case of symmetric stretching, a qualitative change occurs for k ≥ 6 such that accuracy improves as the molecule dissociates further. Single-bond cleavage exhibits a similar, though less pronounced, change for k ≥ 2, as per the discussion in the main text. In contrast, sweeping the bond angle does not reveal any critical value of k. Note that unlike the experiments in the main text, the FCI calculations here do not freeze the oxygen 2p orbital; including this orbital increases the number of allowed bitstrings from 10 to 15, such that a nonzero residual remains for k = 10.

SM.10. Updating Schmidt coefficients
The Schmidt coefficients λ n must be updated each iteration of the VQE optimization. One way to do this is by treating them on the same footing as the gate parameters (ϕ in Fig. 2), such that the optimizer explores a space with dimension equal to the number of Schmidt coefficients plus the number of gate parameters. Here however, since the problems are small, we instead can and do remove the Schmidt coefficients from the optimizer search space by having the classical processor exactly minimize the energy as a function of the Schmidt coefficients based on the latest measurement results. This process is detailed below.
In each iteration of the VQE, quantum circuits are sampled to estimate the expectation values b n |Õ i |b n and φ p bnbm |Õ i |φ p bnbm in Eq. 3, where O i here represents an N -qubit Pauli string. Computing the energy involves plugging these values into a product over i (indicating spin-up or -down), and summations over p, n, m, and the weighted list of 2N -qubit Pauli strings defining the Hamiltonian. Deferring the operations over Schmidt-coefficient indices n, m and performing all others leaves where h is a symmetric matrix of known values. By the variational principle, the eigenvector of this matrix with the lowest eigenvalue minimizes the energy. Thus we update λ by computing this eigenvector, such that the VQE optimizer need only minimize E(ϕ) rather than E(ϕ, λ). Although Schmidt coefficients are typically taken to be non-negative for definiteness, in this computation we allow them to take any real values, facilitating the exact-minimization update routine above and increasing ansatz flexibility for fixed U, V . In the main text, empirical λ values (Fig. 4) are represented by their absolute values for simplicity of discussion.