# Floating block method for quantum Monte Carlo simulations

Avik Sarkar<EMAIL_ADDRESS>Institut für Kernphysik, Institute for
Advanced Simulation and Jülich Center for Hadron Physics, Forschungszentrum
Jülich, D-52425 Jülich, Germany  Facility for Rare Isotope Beams and
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, USA  Dean Lee<EMAIL_ADDRESS>Facility for Rare Isotope Beams and
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, USA  Ulf-G. Meißner<EMAIL_ADDRESS>Helmholtz-Institut
für Strahlen- und Kernphysik and Bethe Center for Theoretical Physics,
Universität Bonn, D-53115 Bonn, Germany  Institut für Kernphysik, Institute
for Advanced Simulation and Jülich Center for Hadron Physics,
Forschungszentrum Jülich, D-52425 Jülich, Germany  Tbilisi State University,
0186 Tbilisi, Georgia

###### Abstract

Quantum Monte Carlo simulations are powerful and versatile tools for the
quantum many-body problem. In addition to the usual calculations of energies
and eigenstate observables, quantum Monte Carlo simulations can in principle
be used to build fast and accurate many-body emulators using eigenvector
continuation or design time-dependent Hamiltonians for adiabatic quantum
computing. These new applications require something that is missing from the
published literature, an efficient quantum Monte Carlo scheme for computing
the inner product of ground state eigenvectors corresponding to different
Hamiltonians. In this work, we introduce an algorithm called the floating
block method, which solves the problem by performing Euclidean time evolution
with two different Hamiltonians and interleaving the corresponding time
blocks. We use the floating block method and nuclear lattice simulations to
build eigenvector continuation emulators for energies of 4He, 8Be, 12C, and
16O nuclei over a range of local and non-local interaction couplings. From the
emulator data, we identify the quantum phase transition line from a Bose gas
of alpha particles to a nuclear liquid.

##  I Introduction

Quantum Monte Carlo simulations are widely used for first principles
calculations of solid state systems and condensed matter Zhang et al. (1997);
Foulkes and Needs (2001); Drut and Lähde (2009); Kolorenč and Mitas (2011);
Ostmeyer et al. (2020); Luu et al. (2022), quantum chemistry Umrigar et al.
(1993); Hammond et al. (1994); Morales et al. (2012), atomic and molecular
physics Lester Jr and Hammond (1990); Booth and Alavi (2010); Morales et al.
(2012), lattice field theories Creutz et al. (1983); Degrand and DeTar (2006);
Gattringer and Lang (2009), nuclear physics Lee (2009); Carlson et al. (2015);
Gezerlis et al. (2014); Lähde and Meißner (2019), degenerate quantum gases
Ceperley and Alder (1980); Astrakharchik et al. (2004); Capogrosso-Sansone et
al. (2008), and other quantum many-body systems. In cases where sign
oscillations are under control, the computational effort scales favorably as a
low-order polynomial in the number of particles. There are many examples of
quantum many-body systems with strong correlations where quantum Monte Carlo
simulations are the only tools currently available for reliable and accurate
predictions. The standard quantities calculated in quantum Monte Carlo
simulations are low-lying energy levels and their corresponding eigenstate
observables. Since wave functions are not constructed explicitly, the set of
linear algebra operations that can be performed on energy eigenstates are
often more limited than that for other methods which explicitly construct and
store wave functions.

One basic linear algebra operation that would be extremely useful is the
calculation of inner products between energy eigenstates corresponding with
different quantum Hamiltonians. Calculating eigenstate inner products within
the framework of quantum Monte Carlo simulations would allow for the
construction of fast and accurate emulators using eigenvector continuation
Frame et al. (2018) for systems that might otherwise be inaccessible using
other methods. Another potential application is the ability to use quantum
Monte Carlo simulations to design time-dependent Hamiltonians $H(t)$ for
efficient adiabatic quantum computing Farhi et al. (2000, 2001). The initial
Hamiltonian $H(0)$ is any trivial Hamiltonian whose ground state can be
prepared on a quantum computer, and the final Hamiltonian $H(T)$ is the
quantum Hamiltonian of interest. Starting from the ground state of $H(0)$ and
slowly evolving with $H(t)$, one can accurately prepare the ground state of
$H(T)$ when $T$ is sufficiently large. Quantum Monte Carlo simulations on
classical computers can be used to optimally select $H(0)$ and the time-
dependent path $H(t)$ such that the inner products between ground states at
times $t$ and $t+dt$ remain as large as possible for each $t$, thereby
reducing the need for very slow time evolution. These calculations require
calculations of inner products between the ground states corresponding to
$H(t)$ and $H(t+dt)$. After preparing the ground state of the desired
Hamiltonian $H(T)$, the quantum computer could then be used to perform real-
time dynamics that would otherwise be beyond classical computing capabilities.

In this work, we introduce a quantum Monte Carlo algorithm called the floating
block method that computes the inner product between two eigenstates produced
by two different Hamiltonians. This is achieved by performing Euclidean time
evolution for the two Hamiltonians and interleaving the corresponding time
blocks. We demonstrate the floating block method using nuclear lattice
simulations and eigenvector continuation to build emulators for the energies
of 4He, 8Be, 12C, and 16O nuclei over a range of local and non-local
interaction couplings.

##  II Floating block method

Let $H_{i}$ and $H_{j}$ be two different quantum Hamiltonians with ground
state wave functions $\ket{v^{0}_{i}}$ and $\ket{v^{0}_{j}}$, respectively.
The only requirement is that $H_{i}$ and $H_{j}$ can be simulated using
quantum Monte Carlo simulations. The floating block method is based on the
identity,

$\lim_{t\rightarrow\infty}\frac{\braket{\psi_{I}}{e^{-H_{i}t}e^{-H_{j}t}e^{-H_{i}t}e^{-H_{j}t}}{\psi_{I}}}{\braket{\psi_{I}}{e^{-2H_{i}t}e^{-2H_{j}t}}{\psi_{I}}}=|\braket{v^{0}_{i}}{v^{0}_{j}}|^{2},$
(1)

where $\ket{\psi_{I}}$ is any initial state that is not orthogonal to both
$\ket{v^{0}_{i}}$ and $\ket{v^{0}_{j}}$. We note that the ground state energy
values for $H_{i}$ and $H_{j}$ drop out of Eq. (1) in the limit of large $t$.
This is a key feature and advantage of the floating block method. While other
quantum Monte Carlo methods were used to calculate
$\braket{v^{0}_{i}}{v^{0}_{j}}$ in Ref. Frame et al. (2018); Frame (2019), the
signal-to-noise ratio dropped rapidly for large quantum systems. The floating
block method avoids these problems completely by canceling the dependence on
the ground state energies. For cases where the phase of the inner product is
nontrivial, we calculate the phase using

$\lim_{t\rightarrow\infty}\frac{\braket{\psi_{I}}{e^{-2H_{i}t}e^{-2H_{j}t}}{\psi_{I}}}{|\braket{\psi_{I}}{e^{-2H_{i}t}e^{-2H_{j}t}}{\psi_{I}}|}=\frac{\braket{v^{0}_{i}}{v^{0}_{j}}}{|\braket{v^{0}_{i}}{v^{0}_{j}}|}.$
(2)

The floating block method can be used for quantum Monte Carlo simulation on
the lattice or in continuous space and is compatible with worldline Monte
Carlo simulations where fermions are explicitly sampled or auxiliary-field
Monte Carlo simulations where fermions are integrated out.

##  III Applications to eigenvector continuation

Eigenvector continuation (EC) is a variational technique for finding the
extremal eigenvectors and eigenvalues of a parameter-dependent Hamiltonian
matrix $H(c)$ Frame et al. (2018). The method relies on the analyticity of the
eigenvectors as function of parameters of the Hamiltonian. The resulting
smooth parametric dependence allows the eigenvector manifold to be well-
approximated by a low-dimensional subspace. By reducing the dimensionality of
the problem, calculations orders of magnitude faster, and this aspect has been
used to build fast and accurate emulators using EC König et al. (2019);
Ekström and Hagen (2019). Its convergence properties were investigated in
Sarkar and Lee (2021), and a greedy algorithm to optimally select training
points was presented in Ref. Sarkar and Lee (2022). In Ref. Bonilla et al.
(2022); Melendez et al. (2022), EC was identified as a specific example of a
wider class of techniques called reduced basis methods Quarteroni et al.
(2016, 2011); Field et al. (2011). Meanwhile, there have been numerous
theoretical developments and applications of eigenvector continuation and
other reduced basis methods in nuclear physics Demol et al. (2020); Frame
(2019); Furnstahl et al. (2020); Bai and Ren (2021); Wesolowski et al. (2021);
Yoshida and Shimizu (2021); Melendez et al. (2021, 2022).

Let us consider a Hamiltonian of the form $H(c)=H_{0}+cH_{1}$. In order to
perform EC calculation of the ground state, we first choose some set training
points $c=c_{0},\cdots,c_{N}$ and find the exact ground state eigenvectors
$\ket{v^{0}_{0}},\cdots,\ket{v^{0}_{N}}$. If our target system corresponds to
parameter value $c=c_{t}$, we calculate the projected Hamiltonian matrix
$H_{ij}=\braket{v^{0}_{i}}{H(c_{t})}{v^{0}_{j}}$, norm matrix
$N_{ij}=\braket{v^{0}_{i}}{v^{0}_{j}}$, and solve the corresponding
generalized eigenvalue problem. We use the floating block method to compute
the norm matrix $N_{ij}$. In order to determine the projected Hamiltonian
matrix elements, we use the result

$\lim_{t\rightarrow\infty}\frac{\braket{\psi_{I}}{e^{-2H_{i}t}H(c_{t})e^{-2H_{j}t}}{\psi_{I}}}{\braket{\psi_{I}}{e^{-2H_{i}t}e^{-2H_{j}t}}{\psi_{I}}}=\frac{H_{ij}}{N_{ij}}.$
(3)

##  IV Methods

Reviews of auxiliary-field lattice Monte Carlo simulations are given in Ref.
Lee (2009); Lähde and Meißner (2019). We start with an initial state and
propagate it over a large but finite number of time steps, which we denote by
$L_{t}$. In each of these time steps, we multiply the current state by the
transfer matrix $M=\exp(-H(c)\Delta t)$. For notational simplicity, we omit
the normal ordering symbols $::$ that should enclose each transfer matrix.
From the data we obtain from different $L_{t}$, we can extrapolate to infinite
time.

Consider a family of Hamiltonians of the form $H(c)=H_{0}+cH_{1}$. For the
lattice simulations presented here, we consider two independent parameters,
but the formalism is the same. We are interested in performing calculation at
some target value point $c=c_{t}$, using the training data at
$c=\\{c_{1},\cdots,c_{k}\\}$. With a finite $L_{t}$ number of time steps, we
write Eq. (3) as,

$\displaystyle\frac{H_{ij}}{N_{ij}}=\frac{\braket{\psi_{I}}{e^{-H(c_{i})\Delta
t}}{\cdots|H(c_{t})|\cdots|e^{-H(c_{j})\Delta
t}|\psi_{I}}}{\braket{\psi_{I}}{e^{-H(c_{i})\Delta
t}}{\cdots|e^{-H(c_{j})\Delta t}|\psi_{I}}}$ (4)

We simplify the notation above to,

$\displaystyle\frac{H_{ij}}{N_{ij}}=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|H(c_{t})|c_{j}|\cdots|c_{j}|\psi_{I}}}{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{j}|\cdots|c_{j}|\psi_{I}}},$
(5)

where $|c_{i}|$ represents a time step where we multiply by transfer matrix
$e^{-H(c_{i})\Delta t}$, and there $L_{t}/2$ time steps of $|c_{i}|$ and
$|c_{j}|$ each. Note again that since there are an equal number of such time
steps in the numerator and denominator, the ground state energy factors cancel
out. The magnitude of the norm matrix $N_{ij}$ is given by

$\displaystyle|N_{ij}|^{2}=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{j}|\cdots|c_{j}|c_{i}|\cdots|c_{i}|c_{j}|\cdots|c_{j}|\psi_{I}}}{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{i}|\cdots|c_{i}|c_{j}|\cdots|c_{j}|c_{j}|\cdots|c_{j}|\psi_{I}}},$
(6)

in the limit that $L_{t}$ goes to infinity. We prove this expression in the
Supplemental Materials.

For our quantum Monte Carlo simulations, we select configurations according to
the absolute value of the expression in the denominator of Eq. (6) and
evaluate the expressions in the numerator and the denominator. However, the
stochastic noise will be large if the expressions in the numerator and
denominator are very different. We first discuss how to optimize the
performance of the floating block method for auxiliary-field Monte Carlo
simulations, and then we discuss how to optimize the process for worldline
Monte Carlo simulations.

We can reduce the stochastic error in auxiliary-field Monte Carlo simulations
with two techniques. The first is to reorder the auxiliary fields in the
numerator. Usually, the same auxiliary field configurations are used for
calculating the numerator and the denominator in Eq. (6). Here we instead
reorder the auxiliary fields in the numerator so that the sequence of
auxiliary fields for coupling $c_{i}$ are the same for numerator and
denominator and the sequence of auxiliary fields for coupling $c_{j}$ are the
same for numerator and denominator. After performing this reordering, the only
difference between the numerator and denominator then comes from the
commutator of the transfer matrices with different Hamiltonian couplings. We
can illustrate the calculation of $|N_{ij}|^{2}$ for 12 time steps,

(7)

The time steps in the numerator and denominator with the same highlighted
boxes have the same auxiliary field configurations.

When the displacement of time blocks is large, the numerator and denominator
are quite different calculations, and the stochastic noise can become quite
large. This is fixed by the second of our two techniques. We allow the time
blocks to “float” gradually, similar to a block of ice detaching from a large
ice mass and floating into the sea. In the following example, we float three
time blocks one step at a time (floating block shown in gray),

(8) (9) (10)

Each term, $N_{1},N_{2},N_{3}$, is computed separately and multiplied together
to give $N_{1}N_{2}N_{3}=|N_{ij}|^{2}$. This gradual floating of time blocks
is also useful for worldline Monte Carlo simulations. For that case, the
worldlines for the numerator and denominator are chosen to be the same, and
the gradual floating of the time blocks produces correlated samples with
reduced stochastic noise.

##  V Results

We perform nuclear lattice simulations using a simple leading order
interaction that is independent of spin and isospin. Despite the simplicity of
the interaction, previous lattice studies have shown that they produce a good
description of the bulk properties of nuclear structure and thermodynamics Lu
et al. (2019, 2020); Shen et al. (2023); Ren et al. (2023). We use a spatial
lattice spacing of 1.97 fm and a time lattice spacing of 0.197 fm$/c$. The
interactions are composed of two types of interactions, which we write as
$V_{L}({\bf r}^{\prime},{\bf r})$ and $V_{NL}({\bf r}^{\prime},{\bf r})$. Here
${\bf r}$ and ${\bf r}^{\prime}$ are the relative separation of the incoming
and outgoing nucleons respectively. For $V_{NL}$ the interaction is non-local,
which means that ${\bf r}$ and ${\bf r}^{\prime}$ are in general different.
Meanwhile, $V_{L}$ is a local interaction where ${\bf r}^{\prime}$ and ${\bf
r}$ are the same. These local and non-local interactions are smeared in space
with parameters $s_{L}$ and $s_{NL}$ respectively, and the normalizations of
$V_{L}$ and $V_{NL}$ are chosen so that the physical 4He binding energy is
reproduced for interactions $H_{A}=K+V_{L}$ and $H_{B}=K+V_{NL}$, where $K$ is
the kinetic energy operator. We describe the interactions in detail in the
Supplemental Material.

We consider general linear combinations of the local and non-local
interactions, $H(c_{L},c_{NL})=K+c_{L}V_{L}+c_{NL}V_{NL}$, and compute the
binding energies of 4He, 8Be, 12C, and 16O. In Table 1, we show EC predictions
for the ground state energy of 12C using the floating block method in
comparison with full simulation results, for $L_{t}=400$. For the EC emulator,
we use the training points $(c_{L},c_{NL})=(0.5,0.5)$, $(0.2,0.8)$, $(0,1)$
and $(0.2,0.6)$, with the corresponding ground state vectors included in that
order. The reported errors for the full simulation are one standard deviation
stochastic errors. The errors for EC calculations are also one standard
deviation errors determined from propagating stochastic errors using the
trimmed sampling method described in Ref. Hicks and Lee (2023). We see that
the EC emulator very accurately reproduces the full simulations results for
test points $(c_{L},c_{NL})$ within the interpolation region of the training
points. For test points outside of the interpolation region, the quality of
the predictions is less accurate but still quite good for most cases. We
should note that these EC calculations would not be possible without the
floating block method. Calculations of the norm matrix elements without the
floating block method would require at least three or four orders of magnitude
more computational effort.

$(c_{L},c_{NL})$ | Full Simulation | 2nd order EC | 3rd order EC | 4th order EC
---|---|---|---|---
(0.8,0.2) | -338.57 ± 0.03 | -330.63 ± 0.26 | -333.15 ± 1.85 | -333.24 ± 1.14
(0.8,0.1) | -295.33 ± 0.02 | -290.04 ± 0.24 | -292.08 ± 1.51 | -292.17 ± 1.03
(0.9,0.1) | -381.81 ± 0.02 | -369.12 ± 0.28 | -372.63 ± 2.24 | -372.77 ± 1.83
(0.8,0.3) | -382.41 ± 0.03 | -371.22 ± 0.27 | -374.24 ± 2.25 | -374.34 ± 1.42
(0.4,0.6) | -177.15 ± 0.05 | -177.25 ± 0.19 | -177.32 ± 0.46 | -177.33 ± 0.24
(0.4,0.7) | -217.73 ± 0.04 | -217.70 ± 0.19 | -217.71 ± 0.38 | -217.72 ± 0.21
(0.3,0.7) | -141.11 ± 0.05 | -139.63 ± 0.19 | -141.09 ± 0.71 | -141.10 ± 0.35
(0.4,0.8) | -259.35 ± 0.06 | -258.19 ± 0.20 | -258.31 ± 0.30 | -258.32 ± 0.28
(0.2,0.6) | -41.54 ± 0.05 | -31.91 ± 0.15 | -37.75 ± 0.50 | -41.64 ± 0.14
(0.2,0.2) | 31.26 ± 0.02 | 101.48 ± 0.17 | 99.18 ± 1.86 | 69.72 ± 1.77
(0.8,0.8) | -606.63 ± 0.05 | -574.23 ± 0.37 | -579.97 ± 4.02 | -580.12 ± 2.70
Table 1: Comparison of the EC emulator prediction and full simulations results
for 12C. The training points added to the EC subspace progressively are:
$(c_{L},c_{NL})=(0.5,0.5)$, $(0.2,0.8)$, $(0,1)$, and $(0.2,0.6)$.

In Fig. 1, we plot the ground state energy of 8Be relative to the two-alpha
threshold in Panel (a), ground state energy of 12C relative to the three-alpha
threshold in Panel (b), and ground state energy of 16O relative to the four-
alpha threshold in Panel (c). For each case, we perform second-order EC with
training points at $(c_{L},c_{NL})=(0.5,0.5)$ and $(0,1)$ in a periodic box of
length $L=15.76$ fm. We see that for $c_{L}=0$ and $c_{NL}\leq 1$, the finite-
volume ground state energies of 8Be, 12C, and 16O are each higher than their
corresponding multi-alpha thresholds. This is consistent with the finding in
Ref. Elhatisari et al. (2016) that, without sufficiently attractive local
interactions, symmetric nuclear matter forms a Bose gas of alpha particles
rather than a nuclear liquid. The location of the quantum phase transition in
the $(c_{L},c_{NL})$ plane can be seen clearly from the zero contours in Fig.
1 where the ground state energy equals the corresponding multi-alpha
threshold.

(a) $E({}^{8}\text{Be})-2E({}^{4}\text{He})$

(b) $E({}^{12}\text{C})-3E({}^{4}\text{He})$

(c) $E({}^{16}\text{O})-4E({}^{4}\text{He})$

Figure 1: Contour plots for the difference between the EC emulated energy and
the corresponding multi-alpha threshold energy for 8Be, 12C, and 16O. The
dashed lines show the experimental values.

##  VI Summary and Discussion

We have introduced the floating block method, a new computational algorithm
for quantum Monte Carlo simulations that allows for the calculation of inner
products between ground-state wave functions corresponding to different
Hamiltonians. Such calculations were not practical using previously existing
methods, and the floating block method provides a computational advantage of
at least three or four orders of magnitude.

The floating block method allows for the efficient construction of EC
emulators using quantum Monte Carlo simulations. In this work, we have used to
explore the $(c_{L},c_{NL})$ phase diagram corresponding to local and non-
local interactions. We find that the EC emulators are very accurate for
interpolation and also fairly reliable for extrapolation. We are able to
locate the quantum phase transition between a Bose gas of alpha particles and
a nuclear liquid.

While the examples we have considered here are lattice calculations using
auxiliary fields, the application of the floating block method for continuum
calculations with worldline Monte Carlo is straightforward. We hope that this
theoretical development will enable the construction of new emulators using
quantum Monte Carlo methods as well as possible novel applications for quantum
computing.

##  VII Acknowledgement

We are grateful for useful discussions with members of the Nuclear Lattice
Effective Field Theory Collaboration. DL was supported in part by the U.S.
Department of Energy (DE-SC0021152, DE-SC0013365, DE-SC0023658, SciDAC-5
NUCLEI Collaboration). UGM was supported by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 101018170), by Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) (project-ID 196253076 - TRR 110), the
Chinese Academy of Sciences (CAS) President’s International Fellowship
Initiative (PIFI) (Grant No. 2018DM0034), Volkswagen Stiftung (Grant No.
93562) and by the MKW NRW (funding code NW21-024-A). Computational resources
provided by the Oak Ridge Leadership Computing Facility through the INCITE
award “Ab-initio nuclear structure and nuclear reactions”, the Southern
Nuclear Science Computing Center, the Gauss Centre for Supercomputing e.V.
(www.gauss-centre.eu) for computing time on the GCS Supercomputer JUWELS at
the Jülich Supercomputing Centre (JSC), and the Institute for Cyber-Enabled
Research at Michigan State University.

## References

  * Zhang et al. (1997) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. B 55, 7464 (1997), URL https://link.aps.org/doi/10.1103/PhysRevB.55.7464.
  * Foulkes and Needs (2001) W. Foulkes and R. Needs, Reviews of Modern Physics 73, 33 (2001).
  * Drut and Lähde (2009) J. E. Drut and T. A. Lähde, Phys. Rev. B 79, 165425 (2009), eprint 0901.0584.
  * Kolorenč and Mitas (2011) J. Kolorenč and L. Mitas, Reports on Progress in Physics 74, 026502 (2011).
  * Ostmeyer et al. (2020) J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Phys. Rev. B 102, 245105 (2020), eprint 2005.11112.
  * Luu et al. (2022) T. Luu, U.-G. Meißner, and L. Razmadze, Phys. Rev. B 106, 195422 (2022), eprint 2204.02742.
  * Umrigar et al. (1993) C. Umrigar, M. Nightingale, and K. Runge, The Journal of chemical physics 99, 2865 (1993).
  * Hammond et al. (1994) B. L. Hammond, W. A. Lester, and P. J. Reynolds, _Monte Carlo methods in ab initio quantum chemistry_ , vol. 1 (World Scientific, 1994).
  * Morales et al. (2012) M. A. Morales, J. McMinis, B. K. Clark, J. Kim, and G. E. Scuseria, Journal of chemical theory and computation 8, 2181 (2012).
  * Lester Jr and Hammond (1990) W. A. Lester Jr and B. L. Hammond, Annual Review of Physical Chemistry 41, 283 (1990).
  * Booth and Alavi (2010) G. H. Booth and A. Alavi, The Journal of chemical physics 132, 174104 (2010).
  * Creutz et al. (1983) M. Creutz, L. Jacobs, and C. Rebbi, Physics Reports 95, 201 (1983).
  * Degrand and DeTar (2006) T. A. Degrand and C. DeTar, _Lattice methods for quantum chromodynamics_ (World Scientific, 2006).
  * Gattringer and Lang (2009) C. Gattringer and C. Lang, _Quantum chromodynamics on the lattice: an introductory presentation_ , vol. 788 (Springer Science & Business Media, 2009).
  * Lee (2009) D. Lee, Progress in Particle and Nuclear Physics 63, 117 (2009), eprint 0804.3501.
  * Carlson et al. (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015), eprint 1412.3081.
  * Gezerlis et al. (2014) A. Gezerlis, I. Tews, E. Epelbaum, M. Freunek, S. Gandolfi, K. Hebeler, A. Nogga, and A. Schwenk, Phys. Rev. C 90, 054323 (2014), URL https://link.aps.org/doi/10.1103/PhysRevC.90.054323.
  * Lähde and Meißner (2019) T. Lähde and U.-G. Meißner, _Nuclear Lattice Effective Field Theory: An Introduction_ , Lecture notes in physics (Springer, 2019), ISBN 9783030141882.
  * Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980), URL https://link.aps.org/doi/10.1103/PhysRevLett.45.566.
  * Astrakharchik et al. (2004) G. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Physical review letters 93, 200404 (2004).
  * Capogrosso-Sansone et al. (2008) B. Capogrosso-Sansone, Ş. G. Söyler, N. Prokof’ev, and B. Svistunov, Physical Review A 77, 015602 (2008).
  * Frame et al. (2018) D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, Phys. Rev. Lett. 121, 032501 (2018), eprint 1711.07090.
  * Farhi et al. (2000) E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv preprint quant-ph/0001106 (2000).
  * Farhi et al. (2001) E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).
  * Frame (2019) D. K. Frame, Ph.D. thesis, Michigan State University (2019), eprint 1905.02782.
  * König et al. (2019) S. König, A. Ekström, K. Hebeler, D. Lee, and A. Schwenk (2019), eprint 1909.08446.
  * Ekström and Hagen (2019) A. Ekström and G. Hagen, Phys. Rev. Lett. 123, 252501 (2019), eprint 1910.02922.
  * Sarkar and Lee (2021) A. Sarkar and D. Lee, Phys. Rev. Lett. 126, 032501 (2021).
  * Sarkar and Lee (2022) A. Sarkar and D. Lee, Phys. Rev. Res. 4, 023214 (2022).
  * Bonilla et al. (2022) E. Bonilla, P. Giuliani, K. Godbey, and D. Lee (2022), eprint 2203.05284.
  * Melendez et al. (2022) J. A. Melendez, C. Drischler, R. J. Furnstahl, A. J. Garcia, and X. Zhang (2022), eprint 2203.05528.
  * Quarteroni et al. (2016) A. Quarteroni, A. Manzoni, and F. Negri, _Reduced Basis Methods for Partial Differential Equations_ (Springer, 2016).
  * Quarteroni et al. (2011) A. Quarteroni, G. Rozza, and A. Manzoni, J.Math.Industry 1 (2011).
  * Field et al. (2011) S. E. Field, C. R. Galley, F. Herrmann, J. S. Hesthaven, E. Ochsner, and M. Tiglio, Phys. Rev. Lett. 106, 221102 (2011).
  * Demol et al. (2020) P. Demol, T. Duguet, A. Ekström, M. Frosini, K. Hebeler, S. König, D. Lee, A. Schwenk, V. Somà, and A. Tichai, Phys. Rev. C101, 041302(R) (2020), eprint 1911.12578.
  * Furnstahl et al. (2020) R. J. Furnstahl, A. J. Garcia, P. J. Millican, and X. Zhang, Phys. Lett. B 809, 135719 (2020), eprint 2007.03635.
  * Bai and Ren (2021) D. Bai and Z. Ren, Phys. Rev. C 103, 014612 (2021), eprint 2101.06336.
  * Wesolowski et al. (2021) S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips (2021), eprint 2104.04441.
  * Yoshida and Shimizu (2021) S. Yoshida and N. Shimizu (2021), eprint 2105.08256.
  * Melendez et al. (2021) J. A. Melendez, C. Drischler, A. J. Garcia, R. J. Furnstahl, and X. Zhang (2021), eprint 2106.15608.
  * Lu et al. (2019) B.-N. Lu, N. Li, S. Elhatisari, D. Lee, E. Epelbaum, and U.-G. Meißner, Phys. Lett. B797, 134863 (2019), eprint 1812.10928.
  * Lu et al. (2020) B.-N. Lu, N. Li, S. Elhatisari, D. Lee, J. E. Drut, T. A. Lähde, E. Epelbaum, and U.-G. Meißner, Phys. Rev. Lett. 125, 192502 (2020), eprint 1912.05105.
  * Shen et al. (2023) S. Shen, S. Elhatisari, T. A. Lähde, D. Lee, B.-N. Lu, and U.-G. Meißner, Nature Commun. 14, 2777 (2023), eprint 2202.13596.
  * Ren et al. (2023) Z. Ren, S. Elhatisari, T. A. Lähde, D. Lee, and U.-G. Meißner (2023), eprint 2305.15037.
  * Hicks and Lee (2023) C. Hicks and D. Lee, Phys. Rev. Res. 5, L022001 (2023).
  * Elhatisari et al. (2016) S. Elhatisari, N. Li, A. Rokash, J. M. Alarcón, D. Du, N. Klein, B.-n. Lu, U.-G. Meißner, E. Epelbaum, H. Krebs, et al., Phys. Rev. Lett. 117, 132501 (2016).
  * Epelbaum et al. (2010) E. Epelbaum, H. Krebs, D. Lee, and U.-G. Meißner, The European Physical Journal A 45, 335 (2010).

##  VIII Supplemental Material

###  VIII.1 Derivation of the floating block method

Let us consider the following ratio of amplitudes,

$\displaystyle
R_{ij}=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{j}|\cdots|c_{j}|c_{i}|\cdots|c_{i}|c_{j}|\cdots|c_{j}|\psi_{I}}}{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{i}|\cdots|c_{i}|c_{j}|\cdots|c_{j}|c_{j}|\cdots|c_{j}|\psi_{I}}},$
(S1)

where we use the notation $|c_{i}|$ for the transfer matrix with coupling
$c_{i}$ and $|c_{j}|$ for the transfer matrix with coupling $c_{j}$. In the
numerator, we have four blocks of $L_{t}/4$ time steps with couplings
$|c_{i}|$ and $|c_{j}|$ interleaved with each other. We now project onto
eigenstates $\sum_{n}\ket{v^{n}_{i}}\bra{v^{n}_{i}}$ for $H(c_{i})$ and
$\sum_{n}\ket{v^{n}_{j}}\bra{v^{n}_{j}}$ for $H(c_{j})$. We then get

$\displaystyle
R_{ij}=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|\sum_{n}\ket{v^{n}_{i}}\bra{v^{n}_{i}}|c_{j}|\cdots|c_{j}|\sum_{n}\ket{v^{n}_{j}}\bra{v^{n}_{j}}|c_{i}|\cdots|c_{i}|\sum_{n}\ket{v^{n}_{i}}\bra{v^{n}_{i}}c_{j}|\cdots|c_{j}|\sum_{n}\ket{v^{n}_{j}}\bra{v^{n}_{j}}\psi_{I}}}{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|\sum_{n}\ket{v^{n}_{i}}\bra{v^{n}_{i}}c_{j}|\cdots|c_{j}|\sum_{n}\ket{v^{n}_{j}}\bra{v^{n}_{j}}\psi_{I}}}.$
(S2)

If we now take the limit as $L_{t}$ goes to infinity, we then have

$\displaystyle
R_{ij}=\frac{\braket{\psi_{I}}{v^{0}_{i}}\braket{v^{0}_{i}}{v^{0}_{j}}\braket{v^{0}_{j}}{v^{0}_{i}}\braket{v^{0}_{i}}{v^{0}_{j}}\braket{v^{0}_{j}}{\psi_{I}}}{\braket{\psi_{I}}{v^{0}_{i}}\braket{v^{0}_{i}}{v^{0}_{j}}\braket{v^{0}_{j}}{\psi_{I}}}=|\braket{v^{0}_{i}}{v^{0}_{j}}|^{2}=|N_{ij}|^{2},$
(S3)

where $\ket{v^{0}_{i}}$ is the ground state eigenvector for $H(c_{i})$ and
$\ket{v^{0}_{j}}$ is the ground state eigenvector for $H(c_{j})$.

###  VIII.2 Lattice interactions

For our lattice calculations, we use a spatial lattice spacing
$a=(100\;\text{MeV})^{-1}=1.97$ fm, and time lattice spacing
$a_{t}=(1000\;\text{MeV})^{-1}=0.197$ fm. We take the nucleon mass to be
$m=938.92$ MeV. We use $\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}}$ to
denote the sum over nearest-neighbor lattice sites of $\mathbf{n}$. Similarly,
$\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}_{i}}$ denotes sum over nearest-
neighbor lattice sites of $\mathbf{n}$ along the $i^{th}$ axis,
$\sum_{\braket{\braket{\mathbf{n^{\prime}}\mathbf{n}}}_{i}}$ denotes sum over
next-to-nearest-neighbor lattice sites of $\mathbf{n}$ along the $i^{th}$
axis, and
$\sum_{\braket{\braket{\braket{\mathbf{n^{\prime}}\mathbf{n}}}}_{i}}$ denotes
sum over next-to-next-to-nearest-neighbor lattice sites of $\mathbf{n}$ along
the $i^{th}$ axis.

For each lattice site $\mathbf{n}$ and smearing parameter $s_{NL}$, we define
non-local annihilation and creation operators,

$\displaystyle a_{NL}(\mathbf{n})$
$\displaystyle=a(\mathbf{n})+s_{NL}\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}}a(\mathbf{n^{\prime}}),$
(S4) $\displaystyle a_{NL}^{\dagger}(\mathbf{n})$
$\displaystyle=a^{\dagger}(\mathbf{n})+s_{NL}\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}}a^{\dagger}(\mathbf{n^{\prime}}),$
(S5)

Similar to the density operator
$\rho(\mathbf{n})=a^{\dagger}(\mathbf{n})a(\mathbf{n})$, we define the smeared
non-local density,

$\displaystyle\rho_{NL}(\mathbf{n})=a_{NL}^{\dagger}(\mathbf{n})a_{NL}(\mathbf{n}).$
(S6)

The smeared local density is defined as

$\displaystyle\rho_{L}(\mathbf{n})=a^{\dagger}(\mathbf{n})a(\mathbf{n})+s_{L}\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}}a^{\dagger}(\mathbf{n^{\prime}})a(\mathbf{n^{\prime}}),$
(S7)

with local smearing parameter $s_{L}$. The local interaction is

$\displaystyle
V_{L}=\frac{q_{L}}{2}\sum_{\mathbf{n}}:\rho_{L}(\mathbf{n})\rho_{L}(\mathbf{n}):,$
(S8)

and the non-local interaction is

$\displaystyle
V_{NL}=\frac{q_{NL}}{2}\sum_{\mathbf{n}}:\rho_{NL}(\mathbf{n})\rho_{NL}(\mathbf{n}):.$
(S9)

The symbols :: indicate normal ordering, where the annihilation operators are
on the right-hand side and the creation operators are on the left-hand side.
The parameters $q_{L},q_{NL},s_{L},$ and $s_{NL}$ are tuned so that both the
interactions $H_{A}=K+V_{L}$ and $H_{B}=K+V_{NL}$ give the physical binding
energy for 4He. The kinetic term is taken to be Epelbaum et al. (2010),

$\displaystyle
K=\frac{49}{12m}\sum_{\mathbf{n}}a^{\dagger}(\mathbf{n})a(\mathbf{n})-\frac{3}{4m}\sum_{\mathbf{n},i}\sum_{\braket{\mathbf{n^{\prime}}\mathbf{n}}_{i}}a^{\dagger}(\mathbf{n})a(\mathbf{n})+\frac{3}{40m}\sum_{\mathbf{n},i}\sum_{\braket{\braket{\mathbf{n^{\prime}}\mathbf{n}}}_{i}}a^{\dagger}(\mathbf{n})a(\mathbf{n})-\frac{1}{180m}\sum_{\mathbf{n},i}\sum_{\braket{\braket{\braket{\mathbf{n^{\prime}}\mathbf{n}}}}_{i}}a^{\dagger}(\mathbf{n})a(\mathbf{n}).$
(S10)

###  VIII.3 Auxiliary field formalism

In our Monte Carlo simulations, we use the auxiliary-field formalism, where we
replace the nucleon-nucleon interactions with interactions of nucleons with
background auxiliary fields at every lattice point. The main idea behind this
method can be understood as a Gaussian integral formula,

$\displaystyle:\exp\Big{(}-\frac{1}{2}ca_{t}\rho^{2}\Big{)}:=\sqrt{\frac{1}{2\pi}}\int_{-\infty}^{\infty}ds:\exp\Big{(}-\frac{1}{2}s^{2}+\sqrt{-ca_{t}}s\rho\Big{)}:.$
(S11)

We can reproduce the interactions by integrating over the auxiliary fields.
For simulating two-, three- and four-body forces on the lattice
simultaneously, we can write the discrete auxiliary fields in the form,

$\displaystyle:\exp\Bigg{(}-\frac{1}{2}Ca_{t}\rho^{2}-\frac{1}{6}C_{3}a_{t}\rho^{3}-\frac{1}{24}C_{4}a_{t}\rho^{4}\Bigg{)}:=\sum_{k=1}^{N}\omega_{k}:\exp\Big{(}\sqrt{-Ca_{t}}\phi_{k}\rho\Big{)}:,$
(S12)

where $C,C_{3}$, and $C_{4}$ are two-, three- and four-body coefficient
respectively, and the :: symbols indicate normal ordering of operators. We use
the ansatz $N=3$, and $w_{k}>0$ for all $k$. We expand Eq. (S12) up to
$O(\rho^{4})$ and compare both sides order by order to solve for $\omega_{k}$
and $\phi_{k}$. For our case, we only have a two-body interaction
($C_{3}=C_{4}=0$), and the solution is
$\phi_{1}=-\phi_{3}=\sqrt{3},\omega_{1}=\omega_{3}=1/6,\omega_{2}=2/3$. We
then sample the corresponding auxiliary fields with the shuttle algorithm, as
described in Ref. Lu et al. (2019).

###  VIII.4 Eigenvector continuation results for 12C

In Fig. S1 we show the EC ground state energy for 12C with two training points
chosen at $(c_{L},c_{NL})=(0.5,0.5)$ and $(0,1)$, with $L_{t}=600$. There is a
sharp bend in the slope of the ground-state energy, and similar results are
seen for the simulations of 8Be and 16O. This sharp bend in the slope of the
energy coincides with the location of the quantum phase transition between a
Bose gas of alpha particles and a nuclear liquid.

Figure S1: EC emulator results for 12C with two training points. The training
points are at $(c_{L},c_{NL})=(0.5,0.5)$ and $(0,1)$, and are marked as X.

###  VIII.5 Previous quantum Monte Carlo approaches

In previous work such as Ref. Frame (2019), the matrix elements were
calculated by sampling auxiliary field configuration according to some
reference coupling $c_{k}$. The calculations to be performed are then

$\displaystyle H_{ij}(c_{t})$
$\displaystyle=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|H(c_{t})|c_{j}|\cdots|c_{j}|\psi_{I}}}{\braket{\psi_{I}}{c_{k}}{\cdots|c_{k}|c_{k}|\cdots|c_{k}|\psi_{I}}},$
(S13) $\displaystyle N_{ij}(c_{t})$
$\displaystyle=\frac{\braket{\psi_{I}}{c_{i}}{\cdots|c_{i}|c_{j}|\cdots|c_{j}|\psi_{I}}}{\braket{\psi_{I}}{c_{k}}{\cdots|c_{k}|c_{k}|\cdots|c_{k}|\psi_{I}}}.$
(S14)

We introduce a stochastic noise since the expressions in the numerator and
denominator are in general different, but the extra normalization factor
present in both $H$ and $N$ matrices drops out when we solve the generalized
eigenvalue problem
$H(c_{t})\ket{\psi(c_{t})}=E(c_{t})N(c_{t})\ket{\psi(c_{t})}$. However, note
that in this approach the couplings in the numerator and denominator are in
general different, and the energy factors do not cancel out. This means that
we get extra factors of $\exp\\{L_{t}[E(c_{k})-E(c_{i})]/2\\}$ and
$\exp\\{L_{t}[E(c_{k})-E(c_{j})]/2\\}$ . For large $L_{t}$, these factors can
be extremely large or small, making the Monte Carlo sampling difficult and
noisy. Adding to the problem is the fact that the factors are different for
different couplings values $c_{i}$ and $c_{j}$.