I_{2}\otimes\bra{x_{1},x_{1}}\right)\ket{x_{0}}\ket{x_{0},x_{1},x_{0},x_{1}}\ket{x_{1}}$
(69) $\displaystyle=2^{-2}\ket{x_{1},x_{0}},$ (70)

where we used the fact that

$\displaystyle G(H,H)\ket{0,0}$
$\displaystyle=2^{-1/2}\sum_{x\in\\{0,1\\}}\ket{x,x}.$ (71)

Overall, this implies

$\displaystyle\braket{0^{n}}{C}{0^{n}}$
$\displaystyle=\braket{0^{n}}{C_{\mathrm{NN}}}{0^{n}}$ (72)
$\displaystyle=\braket{0^{n}}{C_{\mathrm{swp}}}{0^{n}}$ (73)
$\displaystyle=4^{s}\braket{0^{n+4s}}{C_{\mathrm{mgc}}}{0^{n},\mathrm{M}^{s}}.$
(74)

Defining

$\displaystyle\ket{\psi_{\mathrm{NG}}}$
$\displaystyle=\ket{0^{n},\mathrm{M}^{s}},$ (75)
$\displaystyle\ket{\psi_{[}\mathrm{G}}$
$\displaystyle=C_{\mathrm{mgc}}^{\dagger}\ket{0^{n+4s}}$ (76)

completes the proof. ∎

Finally, we get to the proof of the main result.

###### Proof of Theorem 1.

We reduce the problem of estimating the output probability of a quantum
circuit. Specifically, we show that polynomial-time algorithm that produces an
approximation to the magnitude of the overlap
$\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|$ with
multiplicative error $<1/\sqrt{2}$ can be used to construct a polynomial-time
algorithm to produce an approximation to the output probability
${\left|\braket{0^{n}}{C}{0^{n}}\right|}^{2}$ with multiplicative error
$<1/2$, which is $\\#\P$-hard [Corollary 9 [89]].

Consider a circuit $C$ on $n$ qubits and containing $m=\poly(n)$ gates. By
Lemma 4, we can efficiently construct normalized states
$\ket{\psi_{\mathrm{G}}}$ and $\ket{\psi_{\mathrm{NG}}}$ such that
$\braket{0^{n}}{C}{0^{n}}=4^{s}\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}$,
where $s=\poly(n)$ is known from the contsruction. Suppose that we can produce
an estimate $\hat{v}>0$ to
$\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|$ with
multiplicative error $<1/\sqrt{2}$, i.e.,

$\displaystyle\frac{1}{\sqrt{2}}\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|$
$\displaystyle<\hat{c}<\sqrt{2}\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|.$
(77)

Then $\hat{v}^{\prime}=16^{s}\hat{v}^{2}$ is an estimate of
${\left|\braket{0^{n}}{C}{0^{n}}\right|}^{2}$ with multiplicative error
$<1/2$:

$\displaystyle\frac{1}{2}{\left|\braket{0^{n}}{C}{0^{n}}\right|}^{2}={\left(\frac{1}{\sqrt{2}}4^{s}\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|\right)}^{2}$
(78) $\displaystyle<16^{s}\hat{v}^{2}=\hat{v}^{\prime}$ (79)
$\displaystyle<\sqrt{2}\cdot
4^{s}{\left|\braket{\psi_{\mathrm{G}}}{\psi_{\mathrm{NG}}}\right|}^{2}={2}{\left|\braket{0^{n}}{C}{0^{n}}\right|}^{2}.$
(80)

In other words, the ability to compute an estimate $\hat{v}$ to the overlap
with multiplicative error $<1/\sqrt{2}$ yields the ability to compute an
estimate $\hat{v}^{\prime}$ to the output probability with multiplicative
error $<1/2$, which is $\\#\P$-hard. ∎

##  IV Results and Discussion

The focus now shifts to investigating practical implementation and performance
evaluations of MPS-AFQMC. We proved the combination of MPS and AFQMC is
#P-hard but have several classical heuristics to implement MPS-AFQMC. How the
formal complexity limits average-case chemistry applications is unclear, and
this is what we hope to investigate numerically.

The one- and two-dimensional hydrogen lattice models display rich many-body
effects [90] and serve as a testbed for electronic structure methods where
strong correlation is controlled by inter-hydrogen distances [91, 57]. For
instance, the ab initio hydrogen lattice in a minimal basis closely resembles
the Hubbard model, augmented with long-range Coulomb interactions. Larger
basis sets employ many orbitals per site, as in typical materials simulations.
We consider these systems below for numerical explorations of our heuristics.

We utilize PySCF [92] to obtain the electron integrals for the Hamiltonian in
Eq.(2) and conduct preliminary electronic structure calculations such as
Hartree-Fock. These integrals are then supplied to subsequent DMRG and AFQMC
calculations. DMRG calculations are performed with the Renormalizer [93, 63]
package. The MPS-AFQMC algorithm is developed in the AFQMC package ipie [94,
*malone2022ipie].

###  IV.1 Comparisons of the SD-to-MPS strategies

Figure 4: Comparing performance of GMPS-to-MPS, bipartite strategies for the
overlap calculation across different sybstems. (A) The relative error, the
walker MPS’s bond dimension and CPU time of performing SD-to-MPS, averaged
over 3200 equilibrated walkers of H20, as a function of the threshold for the
bipartite approach ($\epsilon_{\textrm{trunc}}$) and the GMPS-to-MPS approach
($\epsilon_{\textrm{ad}}$). The maximal bond dimension of the GMPS-to-MPS
approach is set as 1000. The reference overlap values are taken from the
bipartite approach with truncation threshold $10^{-7}$. (B) The averaged
relative error and averaged CPU time as a function of walker MPS’s bond
dimension for H50. The reference values are taken from the bipartite approach
with truncation threshold $10^{-4}$. The bipartite results are obtained with
$\epsilon_{\textrm{trunc}}\in[10^{-1},10^{-2},10^{-3},10^{-4}]$.

As mentioned, in ph-AFQMC calculations, the overlap between trial and walker
wavefunctions is the central quantity to calculate when performing the
phaseless approximation. In this section, we compare the performance of
different strategies for this task in terms of accuracy and efficiency. The
approaches we consider are GMPS-to-MPS, bipartite, and perfect sampling, as
outlined in Sec. II.3. To see this, in Fig. 4, we employ one-dimensional
hydrogen chains with varying system sizes and inter-atomic distances.

Both bipartite and GMPS-to-MPS approaches involve converting the walker to MPS
with some bond dimension $D_{\textrm{w}}$ before computing the overlap with
the MPS trial. We aim to understand the tradeoff between accuracy and cost set
by different $D_{\textrm{w}}$. First, we compare these two approaches for H20
in a minimal basis with $r=3a_{0}$ ($a_{0}$ represents the Bohr radius). In
Fig. 4(A), we investigate the impact of the threshold
$\epsilon_{\textrm{trunc}}$, which denotes the coefficient truncation
threshold in the bipartite method (as defined in Eq.(36)), and the threshold
$\epsilon_{\textrm{ad}}$, which adaptively controls the block size in the
GMPS-to-MPS conversion. Tightening the threshold for the bipartite strategy
enhances the precision and increases the bond dimension of the walker MPS,
albeit with a slightly increased time cost. Unlike the bipartite approach,
which utilizes a single parameter to control the conversion accuracy, the
GMPS-to-MPS approach relies on two quantities: the adaptive threshold, which
determines the block size and hence the number of gates, and the maximum
allowed bond dimension during compression. Here, we cap the maximum bond
dimension at 1000 and tune only the adaptive threshold
$\epsilon_{\textrm{ad}}$. Tightening the threshold initially leads to a
decrease in error but eventually to an increase in error. This arises because
a smaller threshold leads to a larger block size and, hence, more gates to be
applied, resulting in a higher bond dimension for achieving the same level of
accuracy. For $\epsilon_{\textrm{ad}}\leq 10^{-5}$, the GMPS-to-MPS approach
strikes a good balance between the cost and accuracy. The bipartite method
demonstrates greater cost-effectiveness and accuracy for H20 overall.

For larger systems (e.g., H50), the bipartite approach becomes significantly
more expensive, as illustrated in Fig. 4(B). The error and computational cost
of the bipartite and GMPS-to-MPS approaches with different thresholds are
analyzed as a function of the walker bond dimension in H50. We select
$\epsilon_{\textrm{trunc}}$ values ranging from $0.1,0.01,0.001$ to $0.0001$
for the bipartite approach, where each threshold yields an MPS with a
different bond dimension $D_{\textrm{w}}$. For the GMPS-to-MPS approach, we
maintain the maximum allowed $D_{\textrm{w}}$ at 1000 and adjust the adaptive
threshold $\epsilon_{\textrm{ad}}$. The reference values of the overlap are
obtained from the bipartite approach with $\epsilon_{\textrm{trunc}}=0.0001$.
It is observed that for similar walker bond dimensions $D\textrm{w}$, the
bipartite approach still achieves higher accuracy but at the cost of
substantially longer computational times. For example, when compared to the
GMPS-to-MPS approach with $\epsilon_{\textrm{ad}}=0.001$ and
$D_{\textrm{w}}=200$, the accuracy of the bipartite approach with
$\epsilon_{\textrm{trunc}}=0.01$ is comparable, yet the computational time is
approximately 1000 times slower.

Similarly to H20, tightening the threshold for the GMPS-to-MPS approach
results in decreased accuracy for a fixed bond dimension, $D_{\textrm{w}}$. In
principle, a tighter threshold should lead to improved accuracy, provided that
no approximation is made during the gate applications. However, when we impose
a restriction on the walker’s bond dimension for compressions during the gate
applications, the GMPS-to-MPS method becomes more approximate for that given
threshold. One needs a larger bond dimension to maintain a similar accuracy.
GMPS-to-MPS approach with $\epsilon_{\textrm{ad}}=10^{-3}$ with $D_{w}\leq
200$ seems to strike the balance between cost and accuracy needed for H50.
These observations indicate that the GMPS-to-MPS approach, while potentially
less accurate than the bipartite approach, will likely offer a more scalable
solution for larger systems.

Figure 5: (A) Overlap error from perfect sampling with different system sizes
and (B) the overlap distribution for different system sizes. We use the
nearest neighbor bond length of $r=3a_{0}$ in all cases. In (A), we show the
relative error of the perfect sampling approach for an equilibrated walker of
H10, H20, and H50 as a function of number of samples. The walker state in H50
is projected with a block size of $5$ (see Appendix. D for more information)
for demonstration. The presented distributions are derived from $6.4\times
10^{5}$ equilibrated walkers. For these analyses, the block-averaged
distributions are calculated by averaging groups of 50 walkers per block.

The third strategy circumvents the expensive conversion of the walker state to
MPS and instead samples the overlap value using the perfect sampler. Following
the formalism outlined in Sec.II.3.4, we investigate the perfect sampling
method by selecting an equilibrated walker from systems of varying sizes and
assessing the relative error against the number of samples, as depicted in
Fig.5(A). We sample exactly from the probability distribution given by the
walker matchgate state, which corresponds to Eq. (40). The results are similar
to a strategy where the sampling is performed with the probability from the
trial, which is not shown here for brevity. This technique demands a
significantly higher number of samples for larger systems due to the
exponential decrease of overlap as system size increases. Specifically, the
magnitude of the overlap we examined here —H${}_{10}(r=a_{0})$,
H${}_{10}(r=5a_{0})$, H${}_{20}(r=3a_{0})$, and H${}_{50}(r=3a_{0})$—are on
the order of $10^{-2}$, $10^{-3}$, $10^{-3}$, and $10^{-6}$, respectively. We
show the overlap distribution for systems of H10, H20 and H50 in Fig. 5(B),
where the overlap values become exponentially smaller as system size
increases.

According to Sec. II.3.2, the accuracy of SD-to-MPS via the GMPS approach is
determined by two factors. Although a tighter $\epsilon_{\textrm{ad}}$ leads
to higher accuracy, it concurrently necessitates a substantially larger
$D_{\textrm{w}}$ for convergence. As illustrated in Fig. 6, a balanced choice
employing a moderately small $\epsilon_{\textrm{ad}}$ and an appropriate
$D_{\textrm{w}}$ appears to be the most effective. Based on this, we chose the
GMPS strategy for the rest of the paper. Although careful optimization of
different methods is beyond the scope of this work, it is worth revisiting
this in the future. Additional numerical details for these strategies are
given in the Appendix. B.

Figure 6: Convergence of MPS-AFQMC overlap as a function of system size. All
systems under study are characterized by an interatomic bond distance of
$r=3a_{0}$ and employ the STO-6G basis set. The convergence analysis is
performed using the mean relative error of
$|\langle\Psi_{\textrm{T}}|\psi\rangle|$ over 100 walkers with respect to
$\epsilon_{\textrm{ad}}$ and $D_{\textrm{w}}$.

###  IV.2 Performance of MPS-AFQMC with the GMPS-to-MPS strategy

Given the comparison between various strategies for implementing the overlap
evaluation subroutine for AFQMC, we adopt the GMPS strategy as our primary
approach. In this section, we hope to benchmark how different convergence
parameters in the GMPS strategy affect the final AFQMC energies.

####  IV.2.1 Case study of 1D H10 chain

We test our MPS-AFQMC on a one-dimensional stretched H10 chain. We used a
minimal basis set, STO-6G, and an interatomic distance of 3$a_{0}$ to get
strong correlation effects. With such a small system, it is feasible to
transform the MPS trial into a linear combination of multiple Slater
determinants (MSD) [96], and the MSD-AFQMC results serve as the reference data
for benchmarking MPS-AFQMC,

$|\Psi_{\textrm{T}}\rangle=\sum_{\begin{subarray}{c}s,\
|c_{s}|>\eta\end{subarray}}c_{s}|s\rangle,$ (81)

where the threshold $\eta$ is used to truncate the MSD expansion [96], which
means that SDs with coefficients smaller than $\eta$ are omitted. Here, we set
$\eta=10^{-14}$ to obtain a numerically exact representation of the original
MPS trial.

We use this numerically exact representation to assess the accuracy of MPS-
AFQMC when the overlap could be evaluated exactly. In other words, in this
test, we evaluate the necessary overlap exactly for a given MPS trial. As
illustrated in Fig. 7, with the bond dimension of MPS increasing, the
variational energies gradually improve towards the exact energy. When using
the corresponding MPS as the trial of AFQMC, the deviation from exact full
configuration interaction (FCI) energy decreases. The convergence of MPS-AFQMC
energy to FCI is much faster than MPS itself, demonstrating the synergy
between MPS and AFQMC, even within a small active space. The encouraging
results presented in Fig. 7 are what is expected when SD-to-MPS conversion can
be done exactly or very accurately. However, in practice, the MPS-AFQMC must
make approximations that might cause a large error in the overlap as we saw in
Fig 4 and Fig. 6.

Figure 7: Deviation of MPS trial and MPS-AFQMC from the FCI Energy of H10. The
MPS-AFQMC calculations are conducted without any approximation in SD-to-MPS.
We employ minimal STO-6G atomic orbitals with Löwdin orthogonalization and the
interatomic distance of $r=3a_{0}$. Figure 8: Impact of approximating SD-to-
MPS on AFQMC energy for H10. The heatmaps detail the deviation of the MPS-
AFQMC energy from the converged MPS-AFQMC value shown in Fig. 7.  Figure 9:
Impact of approximating SD-to-MPS on overlap and local energy distribution of
H${}_{\textrm{10}}$. The trial for (A) is the $D_{\textrm{T}}=4$ MPS as in
Fig. 8. The histograms display the distributions of the overlap ratio and the
local energy ratio. The overlap ratio is defined as the ratio of the overlap
calculated with a specific $D_{\textrm{w}}$ and $\epsilon_{\textrm{ad}}$ to
the converged overlap value. Similarly, the local energy ratio compares the
local energy with its converged value. These distributions are aggregated
across the SD of 70,000 equilibrated walkers where spin-projection is used for
walkers [97].

Now, we investigate the impact of these approximations on the final ph-AFQMC
energies using trials with varying quality. Within the GMPS-to-MPS method, we
vary $\epsilon_{\textrm{ad}}$ and $D_{\textrm{w}}$ and compare the results
using three different trials with different bond dimensions
($D_{\textrm{T}}$). Here, the error is measured with respect to the
corresponding MPS-AFQMC results without any approximations made to the overlap
evaluation. As $\epsilon_{\textrm{ad}}$ approaches 0 and $D_{\text{w}}$
approaches $2^{10}$ (maximum possible bond dimension for this problem), we
expect to recover the exact MPS-AFQMC answer. As shown in Fig. 8, for a fixed
$\epsilon_{\textrm{ad}}$, increasing $D_{\textrm{w}}$ leads to a smaller
deviation from the exact MPS-AFQMC energies. A more stringent
$\epsilon_{\textrm{ad}}$ gives a smaller deviation at its converged
$D_{\textrm{w}}$, but it also requires a higher $D_{\textrm{w}}$ to achieve
energy convergence. This is consistent with our observation of the overlap
convergence study in Section IV.1. Similar to the overlap evaluation,
practical applications of GMPS-to-MPS for MPS-AFQMC must balance
$\epsilon_{\text{ad}}$ and $D_{\text{w}}$ to obtain good accuracy without
increasing the cost too much.

In Fig. 9, we examine the distributions of the overlap and local energy ratios
over 70,000 walkers. Here, we use the ratio between the values from specific
$D_{\textrm{w}}$ and $\epsilon_{\textrm{ad}}$ and their converged exact
values. A narrower distribution centered around $1$ indicates smaller relative
errors. For a fixed adaptive threshold $\epsilon_{\textrm{ad}}$, an increase
in $D_{\textrm{w}}$ leads to a narrower distribution of both overlap and local
energy ratios, as expected. In contrast to the overlap, the local energies are
typically much more robust against variations in both $D_{\textrm{w}}$ and
$\epsilon_{\textrm{ad}}$. This robustness largely explains why, despite
sometimes significant deviations in overlap from the exact values, the final
MPS-AFQMC energy deviation from the exact MPS-AFQMC remains small, with most
values falling within the threshold of chemical accuracy (1 m$E_{h}$). The
deviation is even smaller for trials with better quality (higher
$D_{\textrm{T}}$), as shown in Fig. 9. Similar observations are obtained for
H50, as shown in Fig. B3, where the trial is the UHF wavefunction.

The robustness of local energy calculations can be attributed to the zero-
variance principle [98]. More concretely, let us consider the following error
model:

$\displaystyle|\psi_{\text{true}}\rangle$
$\displaystyle=|\psi\rangle+|\epsilon\rangle$ (82)
$\displaystyle\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle$
$\displaystyle=\langle\Psi_{\textrm{T}}|\psi\rangle+\langle\Psi_{\textrm{T}}|\epsilon\rangle$
(83)

where $|\epsilon\rangle$ quantifies the error in SD-to-MPS and the deviation
from the true walker wavefunction. We see that the relative overlap error,
$\langle\Psi_{\textrm{T}}|\epsilon\rangle/\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle$,
can easily become large because the true value in the denominator can be
small. Now, we consider another error model for the local energy evaluation:

$\displaystyle|\Psi_{\textrm{T}}\rangle$
$\displaystyle=|\Psi_{0}\rangle+|\epsilon_{H}\rangle$ (84)
$\displaystyle\langle\Psi_{\textrm{T}}|\hat{H}|\psi_{\text{true}}\rangle$
$\displaystyle=E_{0}\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle+\langle\epsilon_{H}|(\hat{H}-E_{0})|\psi_{\text{true}}\rangle$
(85)

where $|\epsilon_{H}\rangle$ quantifies the distance between
$|\Psi_{\text{T}}\rangle$ and the exact ground state, $|\Psi_{0}\rangle$. It
is possible that the second term in Eq. 85 is much smaller than the first
term, especially if $|\Psi_{\textrm{T}}\rangle$ is close to
$|\Psi_{0}\rangle$. In that case, dividing Eq. 85 by
$\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle$ will give us an accurate
estimation of the local energy even when
$\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle$ is poorly approximated.
This analysis suggests that MPS-AFQMC may achieve high accuracy even with a
relatively loose adaptive threshold $\epsilon_{\textrm{ad}}$ and a smaller
$D_{\textrm{w}}$. These findings also provide insights for quantum-classical
hybrid quantum Monte Carlo methods [46, 99, 49], where the evaluation of
overlap is noisy and correct only up to some additive error and yet the final
AFQMC energies were found to be accurate.

Figure 10: MPS-AFQMC energy error and walker bond dimension as a function of
imaginary time with various thresholds. The trial used here is the same as in
Fig. 8 with $D_{\textrm{T}}=4$, with the maximum allowed bond dimension set at
2000. Spin-projection was used in the initial walker set-up.  Figure 11:
Walker bond dimension as a function of imaginary time at various system sizes.
The weight-averaged $D_{\textrm{w}}$ for AFQMC walkers employing a
unrestricted Hartree-Fock wavefunction as the trial and initial walker for
H10, H20, and H50 systems, evaluated at various interatomic distances
$r\in\\{2a_{0},3a_{0},4a_{0}\\}$. The adaptive threshold is fixed at
$\epsilon_{\textrm{ad}}=10^{-3}$. Each analysis block encompasses 12 steps,
with the time step set at $0.01E_{h}^{-1}$.

In Fig. 10, we fix the trial at $D_{\textrm{T}}=4$, corresponding to the trial
used in Fig. 8, and present the evolution of MPS-AFQMC energy errors from the
exact MPS-AFQMC at each step and walker’s bond dimension as a function of
imaginary time. We observe that, without a cap on the walker’s bond dimension,
a tighter threshold correlates with reduced errors and increased stationary
bond dimension. The walker bond dimension initially increases and can grow
very fast before reaching a plateau. However, as we showed in Fig. 8 and
analyzed in Fig. 9, one appears to need a much smaller walker bond dimension
to obtain good MPS-AFQMC energies.

####  IV.2.2 System size and bond length dependence

We have previously analyzed the increasing hardness of calculating overlap
values for systems as their size increases, as illustrated in Figs. 4 and 5.
Here, we examine how the walker bond dimension changes during imaginary time
propagation as a function of system size. This result is shown in Fig. 11. The
trend in these imaginary-time trajectories is similar to H10, as we discussed
in Fig. 10. Interestingly, while the interatomic distance does not affect the
final saturated bond dimension, it does influence the rate of increase of the
bond dimension; systems with longer interatomic distances exhibit a slower
growth in bond dimension. Contrary to typical DMRG studies for one-dimensional
systems, where bond dimensions needed for a specified accuracy remain fixed
across molecular chains of varying lengths due to the area law, in AFQMC
calculations, the walker’s bond dimension grows with system size. In DMRG
studies, longer interatomic distances necessitate MPS with smaller bond
dimensions, which also contrasts with AFQMC simulations.

###  IV.3 Applications of MPS-AFQMC

In this section, we apply MPS-AFQMC to a variety of systems, including one-
dimensional long chains with varying interatomic distances, two-dimensional
lattices, and scenarios involving large basis sets that incorporate our MPO-
free algorithm.

####  IV.3.1 H50 in minimal basis

Figure 12: The MPS-AFQMC energy deviation for H50 at different bond distances
and comparison with UHF-AFQMC. (A) The dependence of the MPS-AFQMC energy
deviation on the trial energy deviation. The UHF-AFQMC is also plotted, and
the corresponding UHF trial energy deviation is 0.459$\mathrm{Eh}$. We present
the corresponding overlap/local energy distribution in Fig. B3. (B) The MPS-
AFQMC, MPS, and UHF-AFQMC energy for H50 across different bond distances. The
error bars of MPS-AFQMC calculations are all below $0.1$mEh. The total number
of walkers is set to 640, and the time step used is $0.01E_{\textrm{h}}^{-1}$
with population control performed every five steps.

Using the GMPS approach with $\epsilon_{\textrm{ad}}=10^{-3}$ and
$D_{\text{w}}=50$, we apply MPS-AFQMC to calculate the electronic ground state
energy of H50 at various interatomic distances with STO-6G basis. In Fig. 12,
we calculate the deviation of the ground state energy of H50 from the exact
reference provided in [100].

In Fig. 12(A), we present the MPS-AFQMC energy and the corresponding MPS
variational energy using different MPS trials with varying bond dimensions
$D_{\textrm{T}}$ for H50 at $r=3.2a_{0}$. MPS-AFQMC significantly improves
accuracy compared to the trial’s variational energy. Furthermore, as the trial
quality improves, the accuracy of MPS-AFQMC also improves, which aligns with
our initial results observed on smaller systems, as illustrated in Fig. 7.

Next, we apply MPS-AFQMC to various interatomic distances and compare the
results with the variational MPS energy and UHF-AFQMC taken from Ref. 57.
These results are shown in Fig. 12(B). Instead of fixing the bond dimension of
the MPS trial for all interatomic distances, we pick trials with a similar
energy deviation from the true ground state at each distance. Achieving the
same accuracy for a shorter distance becomes harder for DMRG calculations, as
it requires a larger bond dimension MPS to capture the delocalization of
electrons. Here, the MPS trials have bond dimension
$D_{\textrm{T}}=20,10,10,8,5$ for different interatomic distances
$r=2.4a_{0},2.8a_{0},3.2a_{0},3.6a_{0},4.2a_{0}$, respectively. For
calculating the singlet ground state, MPS-AFQMC’s accuracy could improve by
using spin-adapted MPS [28, *keller2016spin] or spin-projected MPS [71], or by
employing identical spin-up and spin-down walker wavefunctions. For
simplicity, we opted for the latter approach, which is referred to as “spin-
projection” in the AFQMC community [97]. Following Eq. (5), where both spin
matrices use the same propagator, they remain identical throughout the
imaginary time evolution. The errors in MPS-AFQMC energy at relatively long
distances (specifically at $r\in\\{3.2a_{0},3.6a_{0},4.2a_{0}\\}$) fall below
the threshold of chemical accuracy (i.e., 1 m$E_{h}$).

We also explored entanglement projection to potentially speed up the
conversion from SD to MPS and improve accuracy by controlling the entanglement
growth in walker wavefunctions. The efficacy of this technique on AFQMC energy
is preliminarily explored in Fig. D6 in Appendix D, warranting further
investigation for potential enhancements. MPS-AFQMC is as accurate as UHF-
AFQMC at relatively shorter distances and produces much more accurate results
for larger bond lengths when the underlying entanglement is smaller in a
localized basis.

####  IV.3.2 Two dimensional lattice: $4\times 4$ hydrogen lattice

Expanding our exploration beyond one-dimensional systems, we investigate the
effectiveness of our MPS-AFQMC in handling more complex two-dimensional
systems. DMRG requires much larger bond dimensions for two-dimensional systems
than the one-dimensional case, primarily due to entanglement growth beyond the
one-dimensional area law.

Our test system involves a $4\times 4$ square lattice of hydrogen atoms with
equivalent nearest-neighbor distances. The orbitals for both DMRG and MPS-
AFQMC calculations are ordered with the genetic ordering algorithm [19] (as
shown in the inset). The bond dimension of the MPO in the test case is 562,
leading to an enlarged half-rotated trial MPS,
$\hat{H}|\Psi_{\textrm{T}}\rangle$, with a bond dimension of $562\times
D_{\textrm{T}}$ which is impractical for local energy computations. We utilize
variational optimization, starting with approximately compressed MPO and MPS
as initial guesses. This strategy successfully reduces the bond dimension of
the half-rotated trial MPS to 100 or less, enabling significantly more
efficient calculations. For more complex active space models or full space
calculations where constructing the MPO is impractical or compressing the
half-rotated MPS trial proves challenging, we suggest using the MPO-free
approach in Section II.4, for which we later present numerical results in Sec.
IV.3.3. With these details, the ground state energies were computed with
different methods and are shown in Fig. 13.

As shown in Fig. 13, the performance of DMRG in this two-dimensional system
demonstrates a remarkably slow convergence with respect to the bond dimension
toward the exact ground state energy. In comparison, using the MPS trials with
$D_{\textrm{T}}=10$ and $D_{\textrm{T}}=50$, MPS-AFQMC is highly effective for
this system and yields energies within chemical accuracy from the exact answer
in both cases. On the other hand, DMRG reaches similar accuracy at
substantially higher bond dimensions of about $200$. However, employing the
UHF trial leads to slightly worse energy, despite the UHF wavefunction having
slightly better energy than $D_{\textrm{T}}=10$ MPS trial. These calculations
used 128 walkers and a time step of $0.1$E${}_{\textrm{h}}^{-1}$. As
demonstrated in Fig. C4, Population bias and timestep errors were negligible.
This is a promising result because MPS-AFQMC could be used for systems where
MPS alone struggles to achieve high accuracy.

Figure 13: Application of MPS-AFQMC to two-dimensional hydrogen lattice with
$r=4.2a_{0}$. Energy deviation from exact results using MPS-AFQMC, UHF-AFQMC,
and DMRG with different bond dimensions.

####  IV.3.3 Towards larger basis set simulations: H10 with cc-pVDZ

Figure 14:  Application of MPS-AFQMC beyond minimal basis. The interatomic
distance for H10 is set at $3.2a_{0}$. (A) The change in AFQMC energy over
imaginary time evolution, comparing the results from brute-force MPS-AFQMC to
those obtained using finite difference MPS-AFQMC. (B) The energy deviation
from the exact reference. The finite difference step is set to $10^{-5}$ for
the one-body energy term and $10^{-4}$ for the two-body energy term.

Our MPS-AFQMC can use MPS within the active space as a trial for AFQMC but
seamlessly incorporate dynamic correlation outside the active space. To
demonstrate this, we employ MPS-AFQMC to calculate the correlation energy for
H10 with the cc-pVDZ basis set. We generate our MPS trial within an active
space of 10 electrons in 10 orbitals (10e, 10o). For demonstration, we
transform this MPS into a linear combination of Slater determinants, which is
possible because this is done in a small active space. The corresponding
results are illustrated in Fig. 14, where the MPS-AFQMC employing finite
difference refers to the algorithm introduced in Section II.4. We found an
excellent agreement between the finite difference algorithm and the brute
force approach, demonstrating its correctness and future utility.

##  V Conclusions and outlooks

In conclusion, we have, for the first time, realized using MPS trial
wavefunctions within fermionic AFQMC calculations. In MPS-AFQMC, evaluating
the overlap between MPS and a walker state (Slater determinant) is formally
and practically the computational bottleneck. We also proved that the
estimation of this overlap up to a multiplicative error is #P-hard. Despite
the hardness of this problem, we proposed several heuristics to approximate
this overlap evaluation. We found that inaccurate overlap evaluation could
still lead to accurate energy estimations due to the zero variance principle
in ph-AFQMC. The strategy we recommend is the one that applies Givens rotation
gates to a product state to obtain a MPS representation of walker Slater
determinants. This choice was made after examining other approaches, including
brute-force conversion, bipartite decomposition, and perfect bit string
sampling. We also discussed balancing accuracy and computational efficiency
for practical applications with suitable thresholds and maximum walker bond
dimensions.

We then showcased accurate MPS-AFQMC results over prototypical strongly
correlated systems. We tested MPS-AFQMC on minimal and double-zeta basis
models of one-dimensional hydrogen chains and a minimal basis model of a two-
dimensional hydrogen lattice. For these systems, our method significantly
improves over ph-AFQMC with a most commonly used trial wavefunction (i.e.,
single Slater determinant). Furthermore, we also proposed an algorithm to
evaluate force bias and local energy via differentiation of certain overlap
evaluations, which removes the significant overhead introduced by prior works
[46, 49, 82]. This allows for using an MPS trial obtained within the active
space and computing the correlation energy of the entire problem only based on
the overlap within the active space. Importantly, this low-scaling approach is
generalizable to other types of trial functions.

Our findings open new avenues for fermionic simulations where two state-of-
the-art techniques, QMC and tensor network methods, can create synergistic
impacts. Furthermore, integrating GPU acceleration into both the SD-to-MPS
subroutine and the AFQMC algorithm [94] can substantially enhance the
efficiency of MPS-AFQMC. Such an implementation will extend its applicability
to more complex systems. We also noted that our insight partially explains the
success of quantum-classical hybrid QMC methods [46]. We hope that the insight
and findings presented in our work will lead to further developments in state-
of-the-art fermionic simulation methods.

## Acknowledgments

T.J. and J.L. were supported by Harvard University’s startup funds, the DOE
Office of Fusion Energy Sciences “Foundations for quantum simulation of warm
dense matter” project, and Wellcome Leap as part of the Quantum for Bio
Program. T.J. thanks Jiajun Ren for stimulating discussions and Zhigang Shuai
for encouragement and support. Computations were carried out partly on the
FASRC cluster supported by the FAS Division of Science Research Computing
Group at Harvard University. This work also used the Delta system at the
National Center for Supercomputing Applications through allocation CHE230032
and CHE230088 from the Advanced Cyberinfrastructure Coordination Ecosystem:
Services & Support (ACCESS) program, which is supported by National Science
Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

##  Appendix A Additional details on ph-AFQMC

This appendix complements the theory of ph-AFQMC in Section II.1. With the
Cholesky decomposition of the two-electron integrals, the Hamiltonian in Eq.
(2) can be expressed as

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}=\hat{v}_{0}-\frac{1}{2}\sum_{\gamma=1}^{N_{\gamma}}\hat{v}_{\gamma}^{2}$
(86)

where

$\hat{v}_{0}=\sum_{pq}(h_{pq}-\frac{1}{2}\sum_{r}(pr|rq))a_{p}^{\dagger}a_{q}$
(87)
$\hat{v}_{\gamma}=i\sum_{pq}L_{pq}^{\gamma}\hat{a}_{p}^{\dagger}\hat{a}_{q}$
(88)

The short-time propagator is approximated using a Trotter decomposition:

$e^{-\Delta\tau\hat{H}}=e^{-\frac{\Delta\tau}{2}\hat{v}_{0}}e^{\frac{\Delta\tau}{2}\sum\hat{v}_{\gamma}^{2}}e^{-\frac{\Delta\tau}{2}\hat{v}_{0}}+\mathcal{O}\left(\Delta\tau^{2}\right)$
(89)

The Hubbard-Stratonovich transformation enables the rewriting of the two-body
propagator as an integration over the auxiliary fields
$\textbf{x}=(x_{1},x_{2},\cdots,x_{N_{\gamma}})$ as follows,

$e^{\Delta\frac{\tau}{2}\sum\hat{v}_{\gamma}^{2}}=\int\mathrm{d}\mathbf{x}p(\mathbf{x})e^{-\sqrt{\Delta\tau}\mathbf{x}\cdot\mathbf{\hat{v}}}$
(90)

where $p(\mathbf{x})$ is the standard normal distribution of the auxiliary
fields.

$\hat{B}(\mathbf{x})=e^{-\frac{\Delta\tau}{2}\hat{v}_{0}}e^{-\sqrt{\Delta\tau}\mathbf{x}\cdot\mathbf{\hat{v}}}e^{-\frac{\Delta\tau}{2}\hat{v}_{0}}$
(91)

In practice, a mean-field shift
$\left\langle\mathbf{\hat{v}}\right\rangle_{\textrm{T}}\equiv\left\langle\Psi_{\textrm{T}}\left|\mathbf{\hat{v}}\right|\Psi_{\textrm{T}}\right\rangle$
subtraction trick is often used to reduce the statistical fluctuations [39,
102],

$\hat{H}_{1}^{\prime}=\hat{H}_{1}-\sum_{\gamma=1}^{N_{\gamma}}\hat{v}_{\gamma}\left\langle\hat{v}_{\gamma}\right\rangle_{\textrm{T}}+\frac{1}{2}\sum_{\gamma=1}^{N_{\gamma}}\left\langle\hat{v}_{\gamma}\right\rangle_{\textrm{T}}^{2}$
(92)
$\hat{H}_{2}^{\prime}=-\frac{1}{2}\sum_{\gamma=1}^{N_{\gamma}}\left(\hat{v}_{\gamma}-\left\langle\hat{v}_{\gamma}\right\rangle_{\textrm{T}}\right)^{2}$
(93)

which changes the expression of the propagator Eq. (91) accordingly.

The force bias in Eq. (5) is defined as

$\overline{\mathbf{x}}_{i}(\Delta\tau,\tau)=-\sqrt{\Delta\tau}\frac{\left\langle\Psi_{\textrm{T}}\left|\hat{\mathbf{v}}\right|\psi_{i}(\tau)\right\rangle}{\left\langle\Psi_{\textrm{T}}|\psi_{i}(\tau)\right\rangle}.$
(94)

In traditional AFQMC, the periodic orthogonalization of walkers is employed to
eliminate round-off numerical errors, thus maintaining numerical stability, in
our MPS-AFQMC algorithm, the orthogonalization is done at every step as
required by the conversion of SD to MPS (SD-to-MPS). Population control is
periodically utilized to remove walkers with small weights and replicate
walkers with large weights.

##  Appendix B Additional methods and results for SD-to-MPS

Figure B1: (a) Conversion of SD-to-MPS by using brute force long ranged MPO-
MPS application (Eq. (20)). (b) Comparing the brute force MPO-MPS approach and
the Gaussian MPS to MPS approach in terms of relative error of overlap and the
time cost as a function of walker’s bond dimension.

##### Brute-force MPO-MPS application.

The brute-force MPO-MPS application is a straightforward, though not
computationally efficient method [72]. The conversion of SD-to-MPS is
accomplished by sequentially applying a series of creation operators to the
initial vacuum state. These creation operators inherently involve long-range
MPO, which contributes to a rapid increase in the entanglement levels of the
MPS and encounters challenges in efficiently compressing the resultant MPS.
The high entanglement levels make it difficult to represent the state with a
low-dimensional MPS without losing significant information. In Fig. B1, we
show its demanding computational time by comparing with the GMPS-to-MPS
approach.

##### The GMPS-to-MPS approach.

We add some details on the GMPS-to-MPS approach we mentioned in Sec. II.3.2.
The Abelian symmetry of MPS [103, 28] is conserved by employing quantum number
conserving SVD to preserve the numerical stability in Eq. (31). When
converting the SD to an MPS, we observed that a random phase factor, denoted
as $e^{i\theta}$, emerges due to the U(1) gauge freedom. In some applications,
this random phase factor is inconsequential, as it does not affect the
resulting state remaining an eigenstate of the one-body operator [74, 77].
However, in our current context, it affects the phaseless constraint, which
relies on the phase change in the overlap. To eliminate this gauge freedom and
fix the gauge, we multiply the MPS by the following scaling factor
$e^{i\theta}$, where

$\theta=\arg{(\langle\psi_{\textrm{SD}}|\psi_{\textrm{max}}\rangle})-\arg{(\langle\psi_{\textrm{MPS}}|\psi_{\textrm{max}}\rangle)}$
(95)

where $|\psi_{\textrm{max}}\rangle$ is the product state with the largest
overlap with the MPS, which can be deterministically obtained easily by the
sampling algorithm introduced in Section II.3.4.

Another problem is that the compression applied during the contraction of the
gates makes the resulting MPS unnormalized. Therefore, we apply the
normalization after all the gates are contracted, which mitigates the error.
We analyze the impact of normalization on the error mitigation of GMPS-to-MPS.
We assume the GMPS-to-MPS approach produces an error of $|\epsilon\rangle$.
Without doing normalization of the resulting MPS, the relative error of the
overlap value is equal to

$\left|\frac{\langle\Psi_{\textrm{T}}|\epsilon\rangle}{\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle}\right|$
(96)

On the other hand, by doing normalization of the resulting MPS, the overlap
value is expressed as

$\frac{\langle\Psi_{\textrm{T}}|(|\psi_{\text{true}}\rangle+|\epsilon\rangle)}{\sqrt{1+\langle\psi_{\text{true}}|\epsilon\rangle+\langle\epsilon|\psi_{\text{true}}\rangle+\langle\epsilon|\epsilon\rangle}}$
(97)

and by performing first-order Taylor expansion of the denominator, the overlap
is approximated to be

$\displaystyle\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle+\langle\Psi_{\textrm{T}}|\epsilon\rangle-\frac{1}{2}\left(\langle\psi_{\text{true}}|\epsilon\rangle+\langle\epsilon|\psi_{\text{true}}\rangle+\langle\epsilon|\epsilon\rangle\right)$
$\displaystyle\times(\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle+\langle\Psi_{\textrm{T}}|\epsilon\rangle)$
(98)

Assuming a small error vector, which implies negligible second-order
contributions, the relative error can be approximated (up to the first order)
as

$\left|\frac{\langle\Psi_{\textrm{T}}|\epsilon\rangle}{\langle\Psi_{\textrm{T}}|\psi_{\text{true}}\rangle}-\frac{1}{2}\left(\langle\psi_{\text{true}}|\epsilon\rangle+\langle\epsilon|\psi_{\text{true}}\rangle\right)\right|$
(99)

While Eq. 99 appears to have a larger error than Eq. 96, the extra terms in
Eq. 99 are expected to be orders of magnitude smaller since they do not have a
prefactor that is inverse of a small quantity. Given that there is only a
small additional error suggested by our analysis, we chose to normalize the
resulting MPS representation of walkers. Based on several numerical tests, we
found that this choice does not affect the final MPS-AFQMC energies in a
statistically significant way.

The conversion to MPS can be performed with various orderings of the site
basis. However, the final ordering must match the trial MPS ordering to
evaluate their overlap. Fig. (2)(A) and Fig. B2 illustrate two distinct
schemes for SD-to-MPS, both finally maintaining the same ordering. The second
conversion is achieved using SWAP gates and associated parity gates to recover
the alternating ordering [104]. In the specific example depicted in Fig.
B2(B), the parity gate is denoted as:

$\hat{\mathcal{P}}=(-1)^{\sum_{k=1}^{N_{\downarrow}}\hat{\sigma}_{k}^{z\downarrow}\Pi_{l\leq
k}^{N_{\uparrow}}\hat{\sigma}_{l}^{z\uparrow}}(-1)^{\sum_{i=1}^{N_{\uparrow}}\hat{\sigma}_{i}^{z\uparrow}\Pi_{j<i}^{N_{\downarrow}}\hat{\sigma}_{j}^{z\downarrow}}$
(100)

where $\sigma_{i}^{z\uparrow(\downarrow)}$ is the Pauli Z operator for the
$i$th spin-up(down) orbitals. As shown in Figure B2(B), the SD-to-MPS
conversion using SWAP gates is both less accurate and more time-consuming
compared to the alternating ordering method used in the main text of the
manuscript.

Figure B2: (A) Performing the GMPS-to-MPS for spin up and spin down block
separately then reordering; (B) Comparing two ordering schemes with H10 with
interatomic distance $r=5.0a_{0}$; The left panel is the mean absolute errors
compared to for the energy of selected 5 equilibrium AFQMC blocks and the
right panel shows the time cost.

##### The bipartite approach.

For the bipartite approach we mentioned in Sec. II.3.3, a phase emerges for
Eq. (33), $e^{i\theta}V^{\dagger}$, which should be removed by dividing the
first row of $V^{\dagger}$ with the scaling factor det($V^{\dagger}$).

##### The sampling approach.

The lower part in Fig. 2(C) shows using perfect sampler from the MPS to obtain
the probability of occupation for the $i-$th bit of the occupancy bitstring,
given the partial occupation bitstring
$\tilde{s}_{1}\tilde{s}_{2}\cdots\tilde{s}_{i-1}$, the probability of
occupation ($s_{i}=1$) and unoccupation ($s_{i}=0$) is obtained by

$P(s_{i}|\tilde{s}_{1}\tilde{s}_{2}\cdots\tilde{s}_{i-1})=\left(\sum_{a_{i}}(\tilde{A}_{a_{i}}^{s_{i}})^{*}\tilde{A}_{a_{i}}^{s_{i}}\right)/P(\tilde{s}_{1}\tilde{s}_{2}\cdots\tilde{s}_{i-1})$
(101)
$\displaystyle\tilde{A}_{a_{i}}^{s_{i}}=\sum_{a_{j},j<i}A_{a_{1}}^{\tilde{s}_{1}}A_{a_{1}a_{2}}^{\tilde{s}_{2}}\cdots
A_{a_{i-2}a_{i-1}}^{\tilde{s}_{i-1}}A_{a_{i-1}a_{i}}^{s_{i}}$

as digramically expressed as follows,

$P(s_{i}|\tilde{s}_{1}\tilde{s}_{2}\cdots\tilde{s}_{i-1})\Leftarrow\leavevmode\hbox
to32.04pt{\vbox to65.88pt{\pgfpicture\makeatletter\hbox{\hskip
11.58092pt\lower-32.93864pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}}
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0.74609375,0.6640625,0.8359375}\pgfsys@color@rgb@fill{0.74609375}{0.6640625}{0.8359375}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{11.38092pt}\pgfsys@moveto{5.69046pt}{11.38092pt}\pgfsys@curveto{5.69046pt}{14.52371pt}{3.14279pt}{17.07138pt}{0.0pt}{17.07138pt}\pgfsys@curveto{-3.14279pt}{17.07138pt}{-5.69046pt}{14.52371pt}{-5.69046pt}{11.38092pt}\pgfsys@curveto{-5.69046pt}{8.23813pt}{-3.14279pt}{5.69046pt}{0.0pt}{5.69046pt}\pgfsys@curveto{3.14279pt}{5.69046pt}{5.69046pt}{8.23813pt}{5.69046pt}{11.38092pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{11.38092pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{}}{{}}{}\pgfsys@moveto{0.0pt}{11.38092pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{0.0pt}{11.38092pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \par
\par\pgfsys@beginscope\pgfsys@invoke{ }
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgf@tempcolor}{rgb}{1,1,1}\pgfsys@invoke{\lxSVG@setpatternuncolored@{3}{1}{1}{1}
}{}\pgfsys@moveto{0.0pt}{11.38092pt}\pgfsys@moveto{5.69046pt}{11.38092pt}\pgfsys@curveto{5.69046pt}{14.52371pt}{3.14279pt}{17.07138pt}{0.0pt}{17.07138pt}\pgfsys@curveto{-3.14279pt}{17.07138pt}{-5.69046pt}{14.52371pt}{-5.69046pt}{11.38092pt}\pgfsys@curveto{-5.69046pt}{8.23813pt}{-3.14279pt}{5.69046pt}{0.0pt}{5.69046pt}\pgfsys@curveto{3.14279pt}{5.69046pt}{5.69046pt}{8.23813pt}{5.69046pt}{11.38092pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{11.38092pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \par\IfSubStr{ur}{u}{
{}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{17.07138pt}\pgfsys@lineto{0.0pt}{22.76186pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{ur}{d}{ {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{5.69046pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{ur}{l}{ {}{{}}{}
{}{}{}\pgfsys@moveto{-5.69046pt}{11.38092pt}\pgfsys@lineto{-11.38092pt}{11.38092pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{ur}{r}{ {}{{}}{}
{}{}{}\pgfsys@moveto{5.69046pt}{11.38092pt}\pgfsys@lineto{11.38092pt}{11.38092pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } ; {}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0.74609375,0.6640625,0.8359375}\pgfsys@color@rgb@fill{0.74609375}{0.6640625}{0.8359375}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{-11.38092pt}\pgfsys@moveto{5.69046pt}{-11.38092pt}\pgfsys@curveto{5.69046pt}{-8.23813pt}{3.14279pt}{-5.69046pt}{0.0pt}{-5.69046pt}\pgfsys@curveto{-3.14279pt}{-5.69046pt}{-5.69046pt}{-8.23813pt}{-5.69046pt}{-11.38092pt}\pgfsys@curveto{-5.69046pt}{-14.52371pt}{-3.14279pt}{-17.07138pt}{0.0pt}{-17.07138pt}\pgfsys@curveto{3.14279pt}{-17.07138pt}{5.69046pt}{-14.52371pt}{5.69046pt}{-11.38092pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-11.38092pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{{}}{}{{}}{{}}{}\pgfsys@moveto{0.0pt}{-11.38092pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{0.0pt}{-11.38092pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \par
\par\pgfsys@beginscope\pgfsys@invoke{ }
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgf@tempcolor}{rgb}{1,1,1}\pgfsys@invoke{\lxSVG@setpatternuncolored@{3}{1}{1}{1}
}{}\pgfsys@moveto{0.0pt}{-11.38092pt}\pgfsys@moveto{5.69046pt}{-11.38092pt}\pgfsys@curveto{5.69046pt}{-8.23813pt}{3.14279pt}{-5.69046pt}{0.0pt}{-5.69046pt}\pgfsys@curveto{-3.14279pt}{-5.69046pt}{-5.69046pt}{-8.23813pt}{-5.69046pt}{-11.38092pt}\pgfsys@curveto{-5.69046pt}{-14.52371pt}{-3.14279pt}{-17.07138pt}{0.0pt}{-17.07138pt}\pgfsys@curveto{3.14279pt}{-17.07138pt}{5.69046pt}{-14.52371pt}{5.69046pt}{-11.38092pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-11.38092pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope \par\IfSubStr{dr}{u}{
{}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{-5.69046pt}\pgfsys@lineto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{dr}{d}{ {}{{}}{}
{}{}{}\pgfsys@moveto{0.0pt}{-17.07138pt}\pgfsys@lineto{0.0pt}{-22.76186pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{dr}{l}{ {}{{}}{}
{}{}{}\pgfsys@moveto{-5.69046pt}{-11.38092pt}\pgfsys@lineto{-11.38092pt}{-11.38092pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } \IfSubStr{dr}{r}{ {}{{}}{}
{}{}{}\pgfsys@moveto{5.69046pt}{-11.38092pt}\pgfsys@lineto{11.38092pt}{-11.38092pt}\pgfsys@stroke\pgfsys@invoke{
} }{ } ;
{}{{}}{}{}\pgfsys@moveto{11.38092pt}{11.38092pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{} {{}{}}{{}}
{{{}}{{}}}{{}}{{{}}{{}}}{}{{}}{}{}{}{}\pgfsys@moveto{11.38092pt}{11.38092pt}\pgfsys@curveto{20.25739pt}{11.38092pt}{20.25739pt}{-11.38092pt}{11.38092pt}{-11.38092pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}{}{{}}{{}}{}\pgfsys@moveto{0.0pt}{28.45232pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-2.31584pt}{27.59172pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$s_{i}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{}{{}}{}{{}}{{}}{}\pgfsys@moveto{0.0pt}{-28.45232pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-2.31584pt}{-29.31291pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$s_{i}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$
$\cdots$$\tilde{s}_{1}$$\tilde{s}_{2}$$\tilde{s}_{i-1}$$s_{i}$$s_{i}$$=$ (102)

A random number is proposed to determine the value of $s_{i}\in[0,1]$. The
following algorithm describes the process of using the perfect sampler to
compute the overlap between MPS and SD, which uses the sampling probability
from the trial. The algorithm that uses the probability in Eq. (40) can be
adjusted accordingly.

1 Set initial overlap $f=0$;

2 for  _$sample=1$ to $N_{\textrm{samp}}$_ do

3 for  _$i=1$ to $N$_ do

4 Compute the probability matrix of site $i$, get $P_{0}$ and $P_{1}$;

5 Draw $u\sim\text{Uniform}(0,1)$;

6 if  _$u\leq P_{0}$_ then

7 Assign $s_{i}=0$;

8 else

9 Assign $s_{i}=1$;

10 end if

11 Next site sampling;

12

13 end for

14 Obtain $|s\rangle=|s_{1}s_{2}\cdots s_{N}\rangle$ and
$f_{0}=\langle\Psi_{\textrm{T}}|s\rangle$;

15 Compute $f_{1}=\langle s|\psi\rangle$;

16 $f=f+f_{1}/f_{0}$;

17

18 end for

Get the sampled overlap, $f/N_{\textrm{samp}}$.

Algorithm 1 Perfect sampling of the overlap Figure B3: Distributions of the
overlap ratio and the local energy ratio for H50 with the MPS form of UHF
trial. The overlap ratio is defined as the ratio of the overlap calculated
with a specific $D_{\textrm{w}}$ and $\epsilon_{\textrm{ad}}$ to the exact
overlap value, which can be obtained by AFQMC with UHF trial in the Slater
determinant form. Similarly, the local energy ratio compares the local energy
with its converged value. These distributions are aggregated across the SD of
1,280,000 equilibrated walkers.

##  Appendix C The AFQMC convergence test with time step and number of walkers

In Fig. C4 we test the number of walkers and time step used in AFQMC
calculation for the ${4\times 4}$ hydrogen lattice system. When using
different setups of the number of walkers and the time step, the AFQMC
energies remain the same up to statistical error bars.

Figure C4: The UHF-AFQMC Energy for the ${4\times 4}$ hydrogen lattice. The
energy differences associated with varying the number of walkers and the time
step, compared to using 640 walkers and $dt=0.01$ E${}_{\textrm{h}}^{-1}$, are
plotted.

##  Appendix D Low bandwidth projection of correlation matrix

Figure D5: Low entanglement projection. Projecting the SD after imaginary time
evolution to a subspace with fixed block size for the correlation matrix.
Figure D6: Low entanglement projection. The system is H${}_{\textrm{50}}$ with
STO-6G basis and inter-atomic distance of $r=3a_{0}$. The upper panel of (A)
shows the energy deviation from the exact answer during the imaginary time
evolution and the lower panel of (A) gives the weight averaged bond dimension
of walkers. The first stage (green shaded) corresponds to simulation with the
adaptive threshold $\epsilon_{\textrm{ad}}=10^{-3}$ and maximally allowed
walker’s bond dimension $D_{\textrm{w}}=50$, and the second stage turns on the
low entanglement projection with projected block size of 5. (B) The
statistical energy deviation for choosing different projected block sizes with
different adpative threshold. The inset show the weight-averaged bond
dimension after performing low entanglement projection.

A product state can be represented as an MPS with a bond dimension of 1, and
the correlation matrix is characterized by having exactly $N_{e}$ diagonal
elements set to one. These correspond to the occupied atomic orbitals in the
product state. All other elements of the correlation matrix are zero. As the
randomness (Eq. (91)(5)) or off-diagonal correlation in a system increases,
deviating from a purely diagonal form, and the entanglement of the system also
increases. Consequently, a larger bond dimension is required to accurately
represent the state as a MPS, as shown in Fig. 11. This increase in bond
dimension allows the MPS to encapsulate the more intricate correlation
patterns that arise in the system’s state. Consequently, compressing the
walker wavefunction becomes more challenging, necessitating larger block sizes
and a greater number of Givens gates. This, in turn, makes the compression of
the MPS less effective. To address this challenge, we propose the subspace low
entanglement projection during the time evolution, which can also speed up the
conversion of SD-to-MPS.

Drawing inspiration from the method to approximate the Gaussian MPS of SD, the
correlation matrix can be approximated as a band matrix. In this
approximation, the matrix retains significant elements along the diagonal and
in a finite band around it, effectively capturing the dominant correlations.
This band matrix representation simplifies the structure while preserving the
essential features of the correlation matrix, facilitating a more efficient
and tractable representation in the context of MPS. This approach could be
particularly useful when dealing with systems where the relevant correlations
are localized or decay rapidly with distance.:

$\Lambda^{\prime}_{pq}=\begin{cases}\Lambda_{pq}&\text{if }p\leq
q<p+\mathcal{B}\\\ 0&\text{otherwise}\end{cases}$ (103)

and the coefficient matrix of the projected walker $\psi^{\prime}$ is obtained
by the singular value decomposition of its correlation matrix

$\Lambda^{\prime}=VSV^{\dagger}$ (104)

with $\mathcal{B}$ in Eq. (103) restrict the overlap
$\langle\psi|\psi^{\prime}\rangle\geq 1-\epsilon_{\textrm{proj}}$. With a
sufficiently small $\epsilon_{\textrm{proj}}$, $\Lambda^{\prime}$ will exhibit
$N_{e}$ singular values that are nearly equal to 1, and $N-N_{e}$ eigenvalues
that are nearly equal to 0. We can express the approximated walker
wavefunction, denoted as $U^{\prime}$, with columns corresponding to those in
$V$ with corresponding singular values close to 1. A phase introduced from the
projection should be addressed using $\langle\psi|\psi^{\prime}\rangle$. The
resulting $|\psi^{\prime}\rangle$ should multiply the factor $e^{-i\theta}$
where

$\theta=\arg{(\langle\psi|\psi^{\prime}\rangle)}$ (105)

which is similar to what we did in Eq. (95).

This procedure is referred to as low-entanglement subspace projection or low
bandwidth projection of correlation matrix,

$|\psi^{\prime}\rangle=\hat{P}|\psi\rangle$ (106)

which leads to a modification of Eq. (5),

$\left|\psi(\tau+\Delta\tau)\right\rangle=\hat{P}\hat{B}\left(\Delta\tau,\mathbf{x}-\overline{\mathbf{x}}\right)\hat{P}\left|\psi(\tau)\right\rangle$
(107)

We present the schematics in Fig. D5. We apply the projection for every step
and monitor the overlap of $\langle\psi|\psi^{\prime}\rangle\geq
1-\epsilon_{\textrm{proj}}$, and it turns out the $\mathcal{B}$ will gradually
decrease and saturate as a small block size. The low entanglement subspace
projection, as implied by its name, serves to reduce the entanglement of
walker’s SDs, leading to faster convergence with respect to $D_{\textrm{w}}$
and the adaptive threshold $\epsilon_{\textrm{ad}}$ in the GMPS approach. The
projection makes the SD-to-MPS easier. Subsequently, we can obtain a more
accurate MPS representation of the approximated walker wavefunction,
effectively eliminating the errors associated with the conversion of the
wavefunction to the GMPS format. As an example, in Fig. 4(B), we demonstrate
the convergence of accuracy and the time cost concerning the bond dimension
$D_{\textrm{w}}$ of equilibrated walkers for both the bipartite approach and
the GMPS approach with different projected block size. Smaller block sizes
greatly reduce computational time and improve accuracy.

Fig. D6(a) displays the energy deviation from the exact answer during the
imaginary time evolution for H50 with an interatomic distance of $r=3a_{0}$.
The upper panel demonstrates the energy deviation, and the lower panel
presents the corresponding bond dimension at each time. The first stage,
marked with a green shade, corresponds to the simulation with the adaptive
threshold $\epsilon_{\textrm{ad}}=10^{-3}$ and a fixed walker’s bond dimension
$D_{\textrm{w}}=50$. The second stage, marked with a red shade, corresponds to
a projected block size equal to 5. Without projection, there is a bias below
the exact answer. Turning on the low entanglement projected technique with a
block size of 5, as shown in Fig.D6(a), the energy bias can be removed, and
the final statistical error is below 1 $m\textrm{H}$, as demonstrated in
Fig.D6(b). As the statistical evolution progresses, the walker wavefunction
becomes more entangled, which makes converting to MPS more challenging.
Increasing the projected block size for a given adaptive threshold leads to
more accurate results. However, a larger block size also requires a tighter
threshold and a larger bond dimension. For a given block size, tightening the
adaptive threshold will also improve the accuracy of the results.

## References

  * Helgaker  _et al._ [2000] T. Helgaker, P. Jørgensen, and J. Olsen,  _Molecular Electronic‐Structure Theory_ (John Wiley & Sons, 2000).
  * Friesner [2005] R. A. Friesner, Proc. Natl. Acad. Sci. 102, 6648 (2005).
  * Jones [2015] R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
  * Koch and Holthausen [2015] W. Koch and M. C. Holthausen,  _A chemist’s guide to density functional theory_ (John Wiley & Sons, 2015).
  * Lee  _et al._ [2023] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, M. Kastoryano, R. Babbush, J. Preskill, D. R. Reichman, E. T. Campbell, E. F. Valeev, L. Lin, and G. K.-L. Chan, Nat. Commun. 14, 1952 (2023).
  * Dagotto [2005] E. Dagotto, Science 309, 257 (2005).
  * Li  _et al._ [2019] Z. Li, S. Guo, Q. Sun, and G. K.-L. Chan, Nat. Chem. 11, 1026 (2019).
  * Sayfutyarova and Hammes-Schiffer [2019] E. R. Sayfutyarova and S. Hammes-Schiffer, J. Chem. Theory Comput. 15, 1679 (2019).
  * Goings  _et al._ [2022] J. J. Goings, A. White, J. Lee, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush, and N. C. Rubin, Proc. Natl. Acad. Sci. 119, e2203533119 (2022).
  * Cui  _et al._ [2022] Z.-H. Cui, H. Zhai, X. Zhang, and G. K.-L. Chan, Science 377, 1192 (2022).
  * Verstraete  _et al._ [2023] F. Verstraete, T. Nishino, U. Schollwöck, M. C. Bañuls, G. K. Chan, and M. E. Stoudenmire, Nat. Rev. Phys. 5, 273 (2023).
  * Baiardi and Reiher [2020] A. Baiardi and M. Reiher, J. Chem. Phys. 152, 040903 (2020).
  * Schollwöck [2011] U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Chan and Sharma [2011] G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011).
  * Austin  _et al._ [2012] B. M. Austin, D. Y. Zubarev, and W. A. J. Lester, Chem. Rev. 112, 263 (2012).
  * Becca and Sorella [2017] F. Becca and S. Sorella,  _Quantum Monte Carlo Approaches for Correlated Systems_ (Cambridge University Press, 2017).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * White [1993] S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Olivares-Amaya  _et al._ [2015] R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang, and G. K. Chan, J. Chem. Phys. 142 (2015).
  * Wolf  _et al._ [2008] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Phys. Rev. Lett. 100, 070502 (2008).
  * Vidal  _et al._ [2003] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. 90, 227902 (2003).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
  * Orús [2014] R. Orús, Ann. Phys. 349, 117 (2014).
  * Hachmann  _et al._ [2007] J. Hachmann, J. J. Dorando, M. Avilés, and G. K. Chan, J. Chem. Phys. 127 (2007).
  * Marti  _et al._ [2008] K. H. Marti, I. M. Ondík, G. Moritz, and M. Reiher, J. Chem. Phys. 128 (2008).
  * Ghosh  _et al._ [2008] D. Ghosh, J. Hachmann, T. Yanai, and G. K. Chan, J. Chem. Phys. 128 (2008).
  * Kurashige and Yanai [2009] Y. Kurashige and T. Yanai, J. Chem. Phys. 130 (2009).
  * Sharma and Chan [2012] S. Sharma and G. K.-L. Chan, J. Chem. Phys. 136, 124121 (2012).
  * Ma and Ma [2013] Y. Ma and H. Ma, J. Chem. Phys. 138 (2013).
  * Yanai  _et al._ [2010] T. Yanai, Y. Kurashige, E. Neuscamman, and G. K. Chan, J. Chem. Phys. 132 (2010).
  * Kurashige and Yanai [2011] Y. Kurashige and T. Yanai, J. Chem. Phys. 135 (2011).
  * Wouters  _et al._ [2013] S. Wouters, N. Nakatani, D. Van Neck, and G. K.-L. Chan, Phys. Rev. B 88, 075122 (2013).
  * Sharma  _et al._ [2014] S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan, Nat. Chem. 6, 927 (2014).
  * Sharma and Chan [2014] S. Sharma and G. K. Chan, J. Chem. Phys. 141 (2014).
  * Cheng  _et al._ [2022] Y. Cheng, Z. Xie, and H. Ma, J. Phys. Chem. Lett. 13, 904 (2022).
  * Larsson  _et al._ [2022] H. R. Larsson, H. Zhai, K. Gunst, and G. K.-L. Chan, J. Chem. Theory Comput. 18, 749 (2022).
  * Zhang  _et al._ [1995] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995).
  * Zhang and Krakauer [2003] S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
  * Motta and Zhang [2018] M. Motta and S. Zhang, Wiley Interdiscip Rev Comput Mol Sci 8, e1364 (2018).
  * Lee  _et al._ [2022a] J. Lee, H. Q. Pham, and D. R. Reichman, J. Chem. Theory Comput. 18, 7024 (2022a).
  * Lee and Reichman [2020] J. Lee and D. R. Reichman, J. Chem. Phys. 153 (2020), 10.1063/5.0015077.
  * Lee  _et al._ [2021a] J. Lee, M. A. Morales, and F. D. Malone, J. Chem. Phys. 154 (2021a).
  * Lee  _et al._ [2021b] J. Lee, S. Zhang, and D. R. Reichman, Phys. Rev. B 103, 115123 (2021b).
  * Mahajan  _et al._ [2022] A. Mahajan, J. Lee, and S. Sharma, J. Chem. Phys. 156, 174111 (2022).
  * Mahajan  _et al._ [2023] A. Mahajan, J. S. Kurian, J. Lee, D. R. Reichman, and S. Sharma, J. Chem. Phys. 159 (2023).
  * Huggins  _et al._ [2022] W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee, Nature 603, 416 (2022).
  * Wan  _et al._ [2023] K. Wan, W. J. Huggins, J. Lee, and R. Babbush, Commun. Math. Phys 404, 629 (2023).
  * Amsler  _et al._ [2023] M. Amsler, P. Deglmann, M. Degroote, M. P. Kaicher, M. Kiser, M. Kühn, C. Kumar, A. Maier, G. Samsonidze, A. Schroeder, _et al._ , J. Chem. Phys. 159 (2023).
  * Kiser  _et al._ [2023] M. Kiser, A. Schroeder, G.-L. R. Anselmetti, C. Kumar, N. Moll, M. Streif, and D. Vodola, New J. Phys. (2023).
  * Mahajan and Sharma [2021] A. Mahajan and S. Sharma, J. Chem. Theory Comput. 17, 4786 (2021).
  * Vitali  _et al._ [2019] E. Vitali, P. Rosenberg, and S. Zhang, Phys. Rev. A 100, 023621 (2019).
  * Chang  _et al._ [2016] C.-C. Chang, B. M. Rubenstein, and M. A. Morales, Phys. Rev. B 94, 235144 (2016).
  * Kanno  _et al._ [2023] S. Kanno, H. Nakamura, T. Kobayashi, S. Gocho, M. Hatanaka, N. Yamamoto, and Q. Gao, “Quantum computing quantum Monte Carlo with hybrid tensor network for electronic structure calculations,”  (2023).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Al-Saidi  _et al._ [2007] W. A. Al-Saidi, S. Zhang, and H. Krakauer, J. Chem. Phys. 127, 144101 (2007).
  * Simons Collaboration on the Many-Electron Problem  _et al._ [2020a] Simons Collaboration on the Many-Electron Problem, K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo, R. J. Anderson, J. Li, L. N. Tran, C.-N. Yeh, B. Mussard, S. Sharma, F. Bruneval, M. van Schilfgaarde, G. H. Booth, G. K.-L. Chan, S. Zhang, E. Gull, D. Zgid, A. Millis, C. J. Umrigar, and L. K. Wagner, Phys. Rev. X 10, 011041 (2020a).
  * Simons Collaboration on the Many-Electron Problem  _et al._ [2017] Simons Collaboration on the Many-Electron Problem, M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, A. J. Millis, N. V. Prokof’ev, U. Ray, G. E. Scuseria, S. Sorella, E. M. Stoudenmire, Q. Sun, I. S. Tupitsyn, S. R. White, D. Zgid, and S. Zhang, Phys. Rev. X 7, 031059 (2017).
  * Wouters  _et al._ [2014] S. Wouters, B. Verstichel, D. Van Neck, and G. K.-L. Chan, Phys. Rev. B 90, 045104 (2014).
  * Qin [2020] M. Qin, Phys. Rev. B 102, 125143 (2020).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Ann. Phys. 411, 167998 (2019).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Ren  _et al._ [2022] J. Ren, W. Li, T. Jiang, Y. Wang, and Z. Shuai, Wiley Interdiscip Rev Comput Mol Sci 12, e1614 (2022).
  * Chan  _et al._ [2016] G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White, J. Chem. Phys. 145, 014102 (2016).
  * Ren  _et al._ [2020] J. Ren, W. Li, T. Jiang, and Z. Shuai, J. Chem. Phys. 153, 084118 (2020).
  * Hubbard [1959] J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
  * Stratonovich [1958] R. Stratonovich, Sov. Phys. Dokl. 2, 416 (1958).
  * Thouless [1960] D. J. Thouless, Nuclear Physics 21, 225 (1960).
  * Thouless [1961] D. Thouless, Nuclear Physics 22, 78 (1961).
  * Jordan and Wigner [1928] P. Jordan and E. P. Wigner, Z. Phys. 47, 631 (1928).
  * Li and Chan [2017] Z. Li and G. K.-L. Chan, J. Chem. Theory Comput. 13, 2681 (2017).
  * Wu  _et al._ [2020] Y.-H. Wu, L. Wang, and H.-H. Tu, Phys. Rev. Lett. 124, 246401 (2020).
  * Jin  _et al._ [2020] H.-K. Jin, H.-H. Tu, and Y. Zhou, Phys. Rev. B 101, 165135 (2020).
  * Fishman and White [2015] M. T. Fishman and S. R. White, Phys. Rev. B 92, 075132 (2015).
  * Peschel [2012] I. Peschel, Braz. J. Phys. 42, 267 (2012).
  * Klich [2006] I. Klich, J. Phys. A Math. Gen. 39, L85 (2006).
  * Petrica  _et al._ [2021] G. Petrica, B.-X. Zheng, G. K.-L. Chan, and B. K. Clark, Phys. Rev. B 103, 125161 (2021).
  * Ferris and Vidal [2012] A. J. Ferris and G. Vidal, Phys. Rev. B 85, 165146 (2012).
  * Sandvik and Vidal [2007] A. W. Sandvik and G. Vidal, Phys. Rev. Lett. 99, 220602 (2007).
  * Schuch  _et al._ [2008] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 040501 (2008).
  * Bravyi  _et al._ [2022] S. Bravyi, D. Gosset, and Y. Liu, Phys. Rev. Lett. 128, 220503 (2022).
  * Huang  _et al._ [2024] B. Huang, Y.-T. Chen, B. Gupt, M. Suchara, A. Tran, S. McArdle, and G. Galli, arXiv (2024), 2404.18303 .
  * Gottesman [1998] D. Gottesman,  _The Heisenberg representation of quantum computers_, Tech. Rep. LA-UR-98-2848; CONF-980788- (Los Alamos National Lab. (LANL), Los Alamos, NM (United States), 1998).
  * Cai  _et al._ [2018] J.-Y. Cai, H. Guo, and T. Williams, Theor. Comput. Sci. 745, 163 (2018).
  * Terhal and DiVincenzo [2002] B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A 65, 032325 (2002).
  * Cai and Lu [2007] J.-Y. Cai and P. Lu, in  _Proceedings of the thirty-ninth annual ACM symposium on Theory of computing_, STOC ’07 (Association for Computing Machinery, New York, NY, USA, 2007) pp. 401–410.
  * Markov and Shi [2008] I. L. Markov and Y. Shi, SIAM Journal on Computing 38, 963 (2008).
  * O’Gorman [2019] B. O’Gorman, in  _DROPS-IDN/v2/document/10.4230/LIPIcs.TQC.2019.10_ (Schloss-Dagstuhl - Leibniz Zentrum für Informatik, 2019).
  * Hangleiter and Eisert [2023] D. Hangleiter and J. Eisert, Rev. Mod. Phys. 95 (2023), 10.1103/revmodphys.95.035001.
  * Simons Collaboration on the Many-Electron Problem  _et al._ [2020b] Simons Collaboration on the Many-Electron Problem, M. Motta, C. Genovese, F. Ma, Z.-H. Cui, R. Sawaya, G. K.-L. Chan, N. Chepiga, P. Helms, C. Jiménez-Hoyos, A. J. Millis, U. Ray, E. Ronca, H. Shi, S. Sorella, E. M. Stoudenmire, S. R. White, and S. Zhang, Phys. Rev. X 10, 031058 (2020b).
  * Stella  _et al._ [2011] L. Stella, C. Attaccalite, S. Sorella, and A. Rubio, Phys. Rev. B 84, 245117 (2011).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020).
  * Ren  _et al._ [2021] J. Ren, W. Li, T. Jiang, Y. Wang, and Z. Shuai, https://github.com/shuaigroup/renormalizer (2021).
  * Malone  _et al._ [2022a] F. D. Malone, A. Mahajan, J. S. Spencer, and J. Lee, https://github.com/JoonhoLee-Group/ipie (2022a).
  * Malone  _et al._ [2022b] F. D. Malone, A. Mahajan, J. S. Spencer, and J. Lee, J. Chem. Theory Comput. 19, 109 (2022b).
  * Lee  _et al._ [2021c] S. Lee, H. Zhai, S. Sharma, C. J. Umrigar, and G. K.-L. Chan, J. Chem. Theory Comput. 17, 3414 (2021c).
  * Purwanto  _et al._ [2008] W. Purwanto, W. Al-Saidi, H. Krakauer, and S. Zhang, The Journal of chemical physics 128 (2008).
  * Assaraf and Caffarel [1999] R. Assaraf and M. Caffarel, Phys. Rev. Lett. 83, 4682 (1999).
  * Lee  _et al._ [2022b] J. Lee, D. R. Reichman, R. Babbush, N. C. Rubin, F. D. Malone, B. O’Gorman, and W. J. Huggins, “Response to ”Exponential challenges in unbiasing quantum Monte Carlo algorithms with quantum computers”,”  (2022b).
  * Hachmann  _et al._ [2006] J. Hachmann, W. Cardoen, and G. K.-L. Chan, J. Chem. Phys. 125, 144101 (2006).
  * Keller and Reiher [2016] S. Keller and M. Reiher, J. Chem. Phys. 144, 134101 (2016).
  * Rom  _et al._ [1997] N. Rom, D. M. Charutz, and D. Neuhauser, Chem. Phys. Lett. 270, 382 (1997).
  * Zgid and Nooijen [2008] D. Zgid and M. Nooijen, J. Chem. Phys. 128, 014107 (2008).
  * Li  _et al._ [2022] W. Li, J. Ren, H. Yang, and Z. Shuai, J. Phys. Condens. Matter 34, 254003 (2022).