# Specialising Neural-network Quantum States for the Bose Hubbard Model

Michael Y Pei1,2 and Stephen R Clark2 1 Bristol Composites Institute, Faculty
of Engineering, University of Bristol, Bristol BS8 1TR, UK
2 H.H. Wills Physics Laboratory, University of Bristol, Bristol BS8 1TL, UK
<EMAIL_ADDRESS>
###### Abstract

Projected variational wavefunctions such as the Gutzwiller, many-body
correlator and Jastrow ansatzes have provided crucial insight into the nature
of superfluid-Mott insulator transition in the Bose Hubbard model (BHM) in two
or more spatial dimensions. However, these ansatzes have no obvious tractable
and systematic way of being improved. A promising alternative is to use
Neural-network quantum states (NQS) based on Restricted Boltzmann Machines
(RBMs). With binary visible and hidden units NQS have proven to be a highly
effective at describing quantum states of interacting
spin-$\textstyle\frac{1}{2}$ lattice systems. The application of NQS to
bosonic systems has so far been based on one-hot encoding from machine
learning where the multi-valued site occupation is distributed across several
binary-valued visible units of an RBM. Compared to
spin-$\textstyle\frac{1}{2}$ systems one-hot encoding greatly increases the
number of variational parameters whilst also making their physical
interpretation opaque. Here we revisit the construction of NQS for bosonic
systems by reformulating a one-hot encoded RBM into a correlation operator
applied to a reference state, analogous to the structure of the projected
variational ansatzes. In this form we then propose a number of specialisations
of the RBM motivated by the physics of the BHM and the ability to capture
exactly the projected variational ansatzes. We analyse in detail the
variational performance of these new RBM variants for a $10\times 10$ BHM,
using both a standard Bose condensate state and a pre-optimised Jastrow +
many-body correlator state as the reference state of the calculation. Several
of our new ansatzes give robust results as nearly good as one-hot encoding
across the regimes of the BHM, but at a substantially reduced cost. Such
specialised NQS are thus primed tackle bosonic lattice problems beyond the
accuracy of classic variational wavefunctions.

##  1 Introduction

For the past two decades Bose-Einstein condensates (BEC) loaded into an
optical lattice have allowed for the creation and study of strongly correlated
systems of atoms [1, 2, 3, 4], in particular the clean realization [2, 5, 6,
7] of the celebrated superfluid (SF) to Mott insulator (MI) transition in the
Bose-Hubbard model (BHM) [8, 9]. Owing to the experimental flexibility of
cold-atom setups, the BHM and its various extensions serve as exquisite test-
beds for investigating interesting effects such as many-body localisation [10,
11, 12, 13, 14, 15], frustration [16, 17] and fractionalisation effects under
artificial magnetic fields [18, 19, 20, 21]. As such there is extensive demand
for accurate numerical descriptions of the ground states of interacting
quantum lattice systems.

For lattice systems in two spatial dimensions and higher there are very few
viable alternatives to using Variational Monte Carlo (VMC) [22, 23] which can
be both a highly efficient and effective method. For the SF-MI transition in
particular, the application of classic Gutzwiller [24], many-body correlator
[25, 26] and Jastrow [27, 28] projected variational ansatzes has provided
crucial insight into the underlying mechanism. Their success stems from the
ability to judiciously capture key physics of the BHM with only a few
variational parameters whose physical meaning is clear. However, it is not
obvious how to tractably and systematically extend these ansatzes so their
accuracy can be improved and their bias steadily removed.

Recently neural-network quantum states (NQS) have emerged as another class of
highly flexible variational ansatzes with many variants such as restricted
Boltzmann machines (RBM) [29], Deep Boltzmann Machines [30, 31, 32], as well
as convolutional [33, 34, 35, 36], feed-forward [37, 38, 39, 40] and recurrent
neural networks [41]. Like the classic bosonic ansatzes, NQS are highly
flexible and can be applied to any number of spatial dimensions. Since the
original use of the RBM ansatz by Carleo and Troyer [29], these NQS have since
been employed in a wide variety of studies, ranging from quantum tomography
[42, 43, 44, 45, 46] to phase classification of wavefunctions [47, 48, 49, 50,
51, 52] to the treatment of open systems [53, 54, 55, 53]. Of particular
importance to VMC applications is the vast representational power of NQS which
has been used to construct exact representations of a wide variety of states
[56, 30, 31, 57, 58, 59, 60, 61].

The application of NQS to systems with a local on-site dimension $d>2$, such
as spin-1 or bosonic systems, has been attacked with RBMs [62, 63, 64] as well
as convolutional and feedforward neural networks [65, 37, 38]. The use of RBMs
[62, 63] in particular to study the BHM has been limited to relatively small
system sizes. A key reason for this is the adoption of so-called “one-hot” (or
unary) encoding [66] commonly deployed in machine learning to handle
multinomial or categorical input variables. Rather than representing a
physical degree of freedom directly with one visible unit this approach
encodes the possible local physical states into a set of binary visible units.
While this encoding leverages the power of binary RBMs it effectively
multiplies the number of variational parameters by the on-site dimension $d$,
significantly increasing the complexity of the optimisation.

In this work we revisit the use of NQS to describe bosonic systems by first
expressing a one-hot RBM in a projected form $\left|\,\Psi_{\rm
var}({\bm{\lambda}})\right\rangle=\hat{C}({\bm{\lambda}})\left|\,\Phi_{\rm
ref}\right\rangle$ shared by the classical bosonic ansatzes, where
$\hat{C}({\bm{\lambda}})$ is a correlation operator controlled by variational
parameters $\bm{\lambda}$ and is applied to a reference state
$\left|\,\Phi_{\rm ref}\right\rangle$. We then propose five new RBM variants
with considerably fewer variational parameters based on truncating terms,
changing the operator basis, and expanding the values taken by the hidden
unit. Each of these changes is motivated by the physics of the BHM and guided
by how they enhance the ability of the variant RBM to exactly capture the
classic bosonic ansatzes. Through a careful analysis of their variational
performance across the SF to MI regimes of the BHM we demonstrate how this
specialisation can substantially reduce the complexity of NQS for bosons,
whilst retaining the ability of NQS to be systematically refined by increasing
the hidden unit number.

The structure of this paper is as follows. In Sec. 2 we introduce the BHM, its
key properties and useful operators for its description. In Sec. 3 we briefly
outline the VMC approach, before reviewing the classic Gutzwiller, many-body
correlator and Jastrow projected variational wavefunctions. In Sec. 4 NQS are
introduced in their original formulation for
spin-$\textstyle\frac{1}{2}$/qubit/hard-core bosons systems, and a number of
useful properties and exact constructions are discussed. This is followed in
Sec. 5 by a generalisation of NQS to soft-core bosons leading to five newly
proposed specialised ansatzes. The performance of these RBM variants across
the SF to MI regimes are then analysed in detail in Sec. 6, before we conclude
in Sec. 7.

##  2 Bose Hubbard Model

For a system of $L$ sites the BHM Hamiltonian [8] comprises a nearest-
neighbour hopping with amplitude $t>0$ and a contact repulsion with strength
$U>0$ as

$\hat{H}=-t\sum_{\langle
i,j\rangle}\hat{b}^{\dagger}_{i}\hat{b}_{j}+\frac{U}{2}\sum_{i=1}^{L}\hat{n}_{i}(\hat{n}_{i}-1),$
(1)

where $\langle i,j\rangle$ denotes nearest-neighbours on the lattice,
$\hat{b}_{i}$ ($\hat{b}^{\dagger}_{i}$) are canonical bosonic annihilation
(creation) operators for site $i$ and
$\hat{n}_{i}=\hat{b}_{i}^{\dagger}\hat{b}_{i}$ is the corresponding bosonic
number operator. The real-space occupation number states of this system

$\left|\,\bm{n}\right\rangle=\left|\,n_{1},n_{2},\dots,n_{L}\right\rangle\propto(\hat{b}^{\dagger}_{1})^{n_{1}}(\hat{b}^{\dagger}_{2})^{n_{2}}\cdots(\hat{b}^{\dagger}_{L})^{n_{L}}\left|\,\textrm{vac}\right\rangle,$
(2)

are defined by a configuration ${\bm{n}}=(n_{1},n_{2},\dots,n_{L})$ of site
occupations $n_{i}\in\mathbb{N}_{B}=\\{0,1,2,\dots,B\\}$ and forms a basis
that spans the $N$-boson Fock space of the system as
$\mathcal{F}_{N,B}=\\{\left|\,\bm{n}\right\rangle~{}|~{}n_{i}\in\mathbb{N}_{B}~{}{\rm
s.t.}~{}\sum_{i=1}^{L}n_{i}=N\\}$. Technically, an exact description of the
BHM requires all the configurations in $\mathcal{F}_{N,N}$ where the maximum
on-site occupation $B=N$. However, for interactions $U/t>1$ configuration
states $\left|\,\bm{n}\right\rangle$ with large multiple occupations possess a
high energy. In this regime there is a negligible error in numerical
calculations working in the subspace spanned by configurations
$\mathcal{F}_{N,B}$ where the maximum occupation of any site is $B=4$ bosons
[9].

In the number basis the off-diagonal matrix elements of
$\langle{\bm{n}}|\hat{H}|{\bm{n}}^{\prime}\rangle$ are real and non-positive,
meaning it is a stoquastic Hamiltonian [67] whose ground state has real non-
negative amplitudes in the same basis. This restriction is strengthened
further by Feynman’s “no node” theorem which proves that the BHM ground state
for any finite $t$ has positive-definite amplitudes [68].

The competing tendencies of hopping to favour delocalising bosons across the
lattice and interactions to penalise multiple occupation of any site, manifest
for any integer filling $\bar{n}=N/L$ as a well-known SF-MI quantum phase
transition [8, 9]. For $U/t=0$ the ground state of Eq. (1) is a Bose-Einstein
condensate (BEC)

$\left|\,\Phi_{0}\right\rangle=\frac{1}{\sqrt{N!}}\left(\frac{1}{\sqrt{L}}\sum_{i=1}^{L}\hat{b}^{\dagger}_{i}\right)^{N}\left|\,\textrm{vac}\right\rangle,$
(3)

in which all $N=L$ bosons populate the $\textbf{q}=0$ quasi-momentum state.
The Poisson on-site number fluctuations with variance
$\Delta^{2}(n_{i})=\langle\Phi_{0}|\hat{n}_{i}^{2}|\Phi_{0}\rangle-\bar{n}^{2}=\bar{n}^{2}$
and long-ranged off-diagonal coherence,
$\langle\Phi_{0}|\hat{b}^{\dagger}_{i}\hat{b}_{j}|\Phi_{0}\rangle=\bar{n}$,
between any pair of sites $i$ and $j$ of this state are idealised features of
the gapless superfluid phase of the model. In the opposite limit $t/U=0$, the
Hamiltonian Eq. (1) decouples into isolated interacting sites and the ground
state is a real-space Fock state with $\bar{n}$ bosons per site

$\left|\,\Phi_{\infty}\right\rangle=\frac{1}{\sqrt{\bar{n}!}}\prod_{i=1}^{L}(\hat{b}^{\dagger}_{i})^{\bar{n}}\left|\,\textrm{vac}\right\rangle=\left|\,\bar{n},\bar{n},\dots,\bar{n}\right\rangle.$
(4)

In this atomic limit all fluctuations in the on-site occupation are frozen out
with
$\Delta^{2}(n_{i})=\langle\Phi_{\infty}|\hat{n}_{i}^{2}|\Phi_{\infty}\rangle-\bar{n}^{2}=0$
and there are no off-diagonal coherences,
$\langle\Phi_{\infty}|\hat{b}^{\dagger}_{i}\hat{b}_{j}|\Phi_{\infty}\rangle=\bar{n}\,\delta_{ij}$.
These are idealised features of the gapped Mott insulating phase of the model.
Between these two limits, at a critical value of the interaction $U_{c}/t$, a
SF-MI transition occurs [69, 70, 71]. In this work we will focus exclusively
on unit filling $\bar{n}=1$.

It is useful to introduce projectors
$\hat{\mathbbm{P}}_{j}^{[x]}\left|\,\bm{n}\right\rangle=\delta_{xn_{j}}\left|\,{\bm{n}}\right\rangle$
on to each local occupation state $\left|\,x\right\rangle$ of site $j$, with
$x\in\mathbb{N}_{B}$. These projectors can be expressed as order-$B$
polynomials of the local number operator $\hat{n}_{j}$

$\displaystyle\hat{\mathbbm{P}}_{j}^{[x]}=\prod_{n\neq
x}^{B}\frac{1}{(x-n)}\left(\hat{n}_{j}-n\right).$ (5)

The completeness of the local occupation basis means that
$\sum_{x=0}^{B}\hat{\mathbbm{P}}^{[x]}_{j}=\mathbbm{1}$, while the number
operator itself is $\hat{n}_{j}=\sum_{x=1}^{B}x\hat{\mathbbm{P}}^{[x]}_{j}$.

With $\bar{n}=1$ and $U/t\gg 1$ the ground state will be almost entirely
contained in the subspace spanned by configurations $\mathcal{F}_{N,2}$. As
such special designations are given to the projectors on to the local states
$\left|\,0\right\rangle$ as holons $\hat{H}_{j}=\hat{\mathbbm{P}}_{j}^{[0]}$,
$\left|\,1\right\rangle$ as singlons $\hat{S}_{j}=\hat{\mathbbm{P}}_{j}^{[1]}$
and $\left|\,2\right\rangle$ as doublons
$\hat{D}_{j}=\hat{\mathbbm{P}}_{j}^{[2]}$. It is also useful to define a
multiplon projector as

$\displaystyle\hat{M}_{j}=\sum_{x=3}^{B}\hat{\mathbbm{P}}_{j}^{[x]},$ (6)

which accounts for all higher than double occupation states.

In this work we will also focus exclusively on a 2D square
$\sqrt{L}\times\sqrt{L}$ lattices with periodic boundary conditions. A simple
mean-field analysis of this 2D system predicts that the Mott insulator melts
below $U_{c}/t=23.2$. More accurate treatments show that the Mott insulator is
stable at lower interaction strengths, with quantum Monte Carlo calculations
for $6\times 6$ lattices estimating [72] instead that $U_{c}/t=16.4$ which is
in good agreement with a strong-coupling expansion prediction of
$U_{c}/t=16.7$ [73]. Consequently, our attention in Sec. 6 will be on the
region between these estimated critical points.

##  3 Variational insight into the SF-MI transition

Variational Monte Carlo has provided crucial insight into the SF-MI
transition. Here we briefly introduce the method and the classic bosonic
ansatzes.

###  3.1 Variational Monte Carlo

A general approach to VMC, which we follow here, starts by constructing a
family of (unnormalised) quantum states [74] as

$\left|\,\Psi_{\rm
var}({\bm{\lambda}})\right\rangle=\hat{C}({\bm{\lambda}})\left|\,\Phi_{\rm
ref}\right\rangle,$ (7)

where the correlation operator $\hat{C}({\bm{\lambda}})$, or correlator, is
controlled by a vector of variational parameters $\bm{\lambda}$ with $n_{\rm
params}$ elements, and applied to a many-body reference state
$\left|\,\Phi_{\rm ref}\right\rangle$. Variational Monte Carlo works in a
fixed basis where the Hamiltonian is sparse. For the BHM the real-space number
basis $\\{\left|\,\bm{n}\right\rangle\\}$ is a good choice. The correlator
$\hat{C}({\bm{\lambda}})$ is taken to be a diagonal operator in this basis
implying it can be expressed as

$\hat{C}({\bm{\lambda}})=\exp\left[-E({\bm{\lambda}},\\{\hat{n}_{i}\\})\right],$
(8)

defined by a pseudo-classical ‘energy’ function $E$ of the number operators
$\hat{n}_{i}$ defined by the variational parameters $\bm{\lambda}$. So long as
both the matrix elements $E({\bm{\lambda}},{\bm{n}})$ and the amplitudes of
the reference state $\left\langle{\bm{n}}|\Phi_{\rm ref}\right\rangle$ can be
efficiently calculated then we are guaranteed that

$\left\langle{\bm{n}}|\Psi_{\rm
var}({\bm{\lambda}})\right\rangle=e^{-E({\bm{\lambda}},{\bm{n}})}\left\langle\bm{n}|\Phi_{\rm
ref}\right\rangle,$ (9)

is also accessible. For the BHM a common choice of reference state is the
$U=0$ BEC ground state $\left|\,\Phi_{0}\right\rangle$ since it possesses a
constant non-zero amplitudes for every configuration
$\left|\,{\bm{n}}\right\rangle$ in $\mathcal{F}_{N,N}$.

Using $\left\langle{\bm{n}}|\Psi_{\rm var}({\bm{\lambda}})\right\rangle$
standard Monte Carlo methods can be applied to find an estimate of the
expectation value of an observable $\hat{A}$ [75, 23], denoted

$\langle\hat{A}\rangle_{\bm{\lambda}}=\frac{\left\langle\,\Psi_{\rm
var}({\bm{\lambda}})\right|{\hat{A}}\left|\,\Psi_{\rm
var}({\bm{\lambda}})\right\rangle}{\left\langle\Psi_{\rm
var}({\bm{\lambda}})|\Psi_{\rm var}({\bm{\lambda}})\right\rangle}.$ (10)

Formally, basis states $\left|\,\bm{n}\right\rangle$ are sampled according to
the probability distribution

$p({\bm{n}})=\frac{|\left\langle{\bm{n}}|\Psi_{\rm
var}({\bm{\lambda}})\right\rangle)|^{2}}{\sum_{\bm{n}}|\left\langle{\bm{n}}|\Psi_{\rm
var}({\bm{\lambda}})\right\rangle|^{2}},$ (11)

and from each sampled state a local estimator of $\hat{A}$ is constructed

$A({\bm{n}})=\sum_{{\bm{n}}^{\prime}}\frac{\left\langle{\bm{n}}^{\prime}|\Psi_{\rm
var}({\bm{\lambda}})\right\rangle}{\left\langle{\bm{n}}|\Psi_{\rm
var}({\bm{\lambda}})\right\rangle}\left\langle\,\bm{n}\right|\hat{A}\left|\,{\bm{n}}^{\prime}\right\rangle,$
(12)

which together gives
$\langle\hat{A}\rangle=\sum_{\bm{n}}p({\bm{n}})A({\bm{n}})$. The sum over
${\bm{n}}^{\prime}$ in Eq. (12) remains tractable once $\hat{A}$ is sparse in
the basis $\\{\left|\,{\bm{n}}\right\rangle\\}$. The distribution
$p({\bm{n}})$ is efficiently estimated using a Markov-chain algorithm [76],
such as Metropolis-Hastings, where only the ratios of amplitudes between two
states $\left|\,{\bm{n}}\right\rangle$ and
$\left|\,{\bm{n}}^{\prime}\right\rangle$ are needed, and so neither
$p({\bm{n}})$ nor the norm $\left\langle\Psi_{\rm
var}({\bm{\lambda}})|\Psi_{\rm var}({\bm{\lambda}})\right\rangle$ are ever
explicitly required, justifying why we can ignore the normalisation of any
ansatz. Moreover, number conservation can be easily handled by only generating
sequences of configuration states $\left|\,{\bm{n}}\right\rangle$ containing
exactly $N$ particles.

Variational minimisation of $\left|\,\Psi_{\rm
var}({\bm{\lambda}})\right\rangle$ can proceed by evaluating its energy
density $\epsilon({\bm{\lambda}})=\langle\hat{H}\rangle_{\bm{\lambda}}/L$ and
its variance, along with their gradient vectors with respects to parameters
${\bm{\lambda}}$, updating them by a small step along the direction of
steepest descent, and iterating until convergence [23]. More sophisticated
approaches such as modified stochastic optimisation [77] and the ‘linear
method’ [78, 79] are also commonly used. Here we use stochastic
reconfiguration [80] which is well suited to interacting quantum lattice
problems like the BHM. The method involves the construction of a $n_{\rm
params}\times n_{\rm params}$ matrix and solving of a set of linear equations.
Given $n_{\rm samp}$ Monte Carlo samples the error on any matrix element will
scales as $1/\sqrt{n_{\rm samp}}$ [22]. So as a rule of thumb $n_{\rm
samp}>10n_{\rm params}$ is at least required to ensure the sampled matrix is
not rank-deficient [22], giving an overall scaling of the algorithm as
$O(n_{\rm params}^{2})$.

The success of VMC is strongly tied to the judicious use of variational states
whose expressiveness is sufficient to capture the expected physics. We now
briefly review the classic bosonic ansatzes that have been successfully
applied to the SF-MI transition in the BHM.

###  3.2 Gutzwiller correlator

Analogous to its fermionic counterpart [81], the simplest trial state for the
BHM interpolating between the BEC and atomic limit is the so-called bosonic
Gutzwiller ansatz given by

$\left|\,\Psi_{\textrm{G}}(g)\right\rangle=\exp\left(-\frac{g}{2}\sum_{i=1}^{L}(\hat{n}_{i}-\bar{n})^{2}\right)\left|\,\Phi_{0}\right\rangle,$
(13)

controlled by a single variational parameter $g>0$. This on-site correlator
can suppress the amplitude of configurations in the BEC reference state where
any site has an occupation greater than $\bar{n}$, thereby mimicking the
expected effect of interactions. The Gutzwiller approach can be solved
analytically and provides the aforementioned mean-field estimate of the
$\bar{n}=1$ superfluid-Mott insulator transition at $U_{c}/t\approx 5.8z$
[24], where $z$ is the coordination number of the lattice. While exact in
infinite dimensions, a significant drawback of Eq. (13) is that the only
insulating state it can reach is $\left|\,\Phi_{\infty}\right\rangle$ when
$g=\infty$. Consequently, it provides an unrealistic description of the entire
Mott insulating phase $U>U_{c}$ by the atomic limit state with a fixed energy
density $\epsilon=0$ and no local density fluctuations.

###  3.3 Many-body correlator

The perturbative correction to the atomic limit motivates how the Gutzwiller
correlator can be embellished. To first order in $t/U$ we have

$\left|\,\Phi_{U}\right\rangle\approx\left|\,\Phi_{\infty}\right\rangle+\frac{t}{U}\sum_{\langle
i,j\rangle}\hat{b}^{\dagger}_{i}\hat{b}_{j}\left|\,\Phi_{\infty}\right\rangle+O((t/U)^{2}),$
(14)

in which a finite kinetic term induces nearest-neighbour density fluctuations.
For $\bar{n}=1$ this induces a transition
$\left|\,11\right\rangle\mapsto\\{\left|\,20\right\rangle,\left|\,02\right\rangle\\}$
and the fluctuations will be adjacent holon-doublon pairs. A size-extensive
way to target these fluctations is to add to the Gutzwiller ansatz a so-called
“many-body” correlator term

$\left|\,\Psi_{\textrm{MB}}(g,\xi)\right\rangle=\exp\left(-\frac{g}{2}\sum_{i=1}^{L}(\hat{n}_{i}-\bar{n})^{2}-\xi\sum_{i=1}^{L}\hat{Q}_{i}\right)\left|\,\Phi_{0}\right\rangle,$
(15)

where an additional variational parameter $\xi>0$ controls the nearest-
neighbour operator [25, 26]

$\hat{Q}_{i}=\hat{H}_{i}\prod_{\langle
i,j\rangle}(1-\hat{D}_{j})+\hat{D}_{i}\prod_{\langle
i,j\rangle}(1-\hat{H}_{j}),$ (16)

constructed from products of doublon and holon projectors. Given that
$\hat{Q}_{i}$ is non-zero only if a holon (doublon) on site $i$ has no
doublons (holons) neighbouring it, the many-body correlator applies a
suppression of $\exp(-\xi)$ to such configurations. In the limit of
$\xi\rightarrow\infty$ this results in the complete binding of holon-doublon
pairs.

The many-body ansatz $\left|\,\Psi_{\rm MB}(g,\xi)\right\rangle$ provides a
more realistic description of the Mott insulating phase by permitting finite
local density fluctuations and an energy density scaling
$\epsilon\propto-t^{2}/U$. It is found to predict a lower critical interaction
of $U_{c}/t\approx 20.5$, with a discontinuity in the optimal value of $\xi$
suggesting a first-order transition [25]. The ansatz also highlights the
importance of bound holon-doublon pairs in the mechanism driving the SF-MI
transition. A simplistic picture might expect that the transition corresponds
to the complete unbinding of these pairs. However, a more detailed study [26]
with the many-body ansatz has revealed that the transition corresponds to when
distance between doublons $\ell_{\textrm{DD}}$ becomes commensurate with
holon-doublon distance $\ell_{\textrm{HD}}$. As depicted in Fig. 1, the
breakdown of the insulating state occurs not when holon-doublon pairs unbind,
but when their density is sufficiently large that hopping between them induces
mobility.

Figure 1: (a) The ground state energy density $\epsilon$ for the 2D square
lattice BHM plotted for the Bogoliubov weakly interacting superfluid
approximation (dotted curve) [82] and its leading order
$\epsilon_{\textrm{weak}}=-4t+U/2$ (dashed line), as well as the strong-
coupling limit $\epsilon_{\textrm{strong}}=-8t^{2}/U$ [83] (solid line). The
grey region indicates the location of the Mott transition bounded by the
quantum Monte Carlo estimate of $U_{c}/t=16.4$ and the mean-field prediction
of $U_{c}/t=23.2$. The $\times$ symbols denote the values of
$U/t=\\{1,16,23,32\\}$ which are considered in this work in Sec. 6. (b) A
schematic depiction of the conjectured mechanism behind the Mott transition. A
2D lattice is shown that is unit filled except for holes, represented by open
circles, and doublons, represented by two shaded circles. The insulating phase
is characterised by a dilute gas of bound holon-doublon pairs where
$\ell_{\textrm{dd}}\gg\ell_{\textrm{hd}}$. The transition to a superfluid
arises when the increasing density of bound pairs causes
$\ell_{\textrm{dd}}\sim\ell_{\textrm{hd}}$ [26].

###  3.4 Jastrow correlator

Since the many-body correlator is restricted to nearest-neighbours it fails to
modify the long-range correlations of BEC reference state crucial for
describing a genuine insulating state. A well known alternative generalising
the Gutzwiller contribution is the two-body Jastrow correlator [84, 27, 28,
80, 85] giving an ansatz

$\left|\,\Psi_{\textrm{J}}(\\{\nu_{ij}\\})\right\rangle=\exp\left(-\frac{1}{2}\sum_{i,j=1}^{L}\nu_{ij}\hat{n}_{i}\hat{n}_{j}\right)\left|\,\Phi_{0}\right\rangle,$
(17)

defined by translationally invariant pseudo-potential parameters
$\nu_{ij}=\nu(|{\bf R}_{i}-{\bf R}_{j}|)$, where ${\bf R}_{i}$ is the real-
space position vector to site $i$. Here $\nu(0)$ contains the on-site
Gutzwiller factors, while longer-ranged contributions in $\nu_{ij}$ can induce
holon-doublon binding over large distances impeding conduction, but still
allowing local density fluctuations. Of course Eq. (17) can be readily
embellished with the many-body correlator also. Moving to a momentum
representation of the Jastrow terms then gives the Jastrow + many-body ansatz

$\left|\,\Psi_{\textrm{J-MB}}(\\{\nu_{\textbf{q}}\\},\xi)\right\rangle=\exp\left(-\frac{1}{2}\sum_{\textbf{q}\in{\rm
BZ}}\nu_{\textbf{q}}\hat{\rho}_{\textbf{q}}\hat{\rho}_{-{\textbf{q}}}-\xi\sum_{i=1}^{L}\hat{Q}_{i}\right)\left|\,\Phi_{0}\right\rangle,$
(18)

described by the Fourier transformed pseudo-potential

$\nu_{\textbf{q}}=\frac{1}{\sqrt{N}}\sum_{i=1}^{L}\nu(|{\textbf{R}}_{i}|)\exp(\textrm{i}{\textbf{q}}\cdot{\textbf{R}}_{i}),$
(19)

and number fluctuation operator

$\hat{\rho}_{\textbf{q}}=\frac{1}{\sqrt{N}}\sum_{i=1}^{L}\hat{n}_{i}\exp(\textrm{i}\textbf{q}\cdot{\textbf{R}}_{i}),$
(20)

for quasi-momentum q in the 1st Brillioun zone of the lattice. The short-
ranged many-body correlator improves the variational energy for $U\gg t$ but
does not affect the singular behaviour of $\nu_{\bf q}$ at long wavelengths
$|{\bf q}|\rightarrow 0$ which is of direct physical relevance [85]. In any
spatial dimension a Jastrow state with $\nu_{\bf q}\sim 1/|{\bf q}|$ has a
finite condensate fraction, algebraically decaying density-density
correlations and gapless low-energy collective excitations expected of a
weakly interacting superfluid. Variational bounds [86] indicate that stronger
singular behaviour in $\nu_{\bf q}$ should be able to completely deplete the
condensate of $\left|\,\Phi_{0}\right\rangle$ and open a gap in excitations as
$\textbf{q}\rightarrow 0$, consistent with the formation of a bosonic
insulating state.

In particular numerical optimisation in 2D found [87] an abrupt change in the
behaviour of $\nu_{\bf q}$ at $U_{c}/t=20.6$ from the superfluid scaling
$\nu_{\bf q}\sim 1/|{\bf q}|$, to $\nu_{\bf q}\sim-\log(|{\bf q}|)/|{\bf
q}|^{2}$, signalling a transition to a Mott insulator. While the latter
$\nu_{\bf q}$ is more singular than the $\nu_{\bf q}\sim 1/|{\bf q}|^{2}$
marginal case for the bounds in 2D, it has not been conclusively shown that it
is singular enough to render density-density correlations exponentially
decaying, as expected for a true insulator.

While very insightful there remains a need to systematically generalise the
ansatz $\left|\,\Psi_{\textrm{J-MB}}\right\rangle$. One approach is to simply
increase the complexity of the pseudo-energy function

$\hat{C}(\lambda)=\exp\left(\sum_{i=1}^{L}E^{(1)}_{i}\hat{n}_{i}+\sum_{i,j=1}^{L}E^{(2)}_{ij}\hat{n}_{i}\hat{n}_{j}+\sum_{i,j,k=1}^{L}E^{(3)}_{ijk}\hat{n}_{i}\hat{n}_{j}\hat{n}_{k}+\cdots\right),$
(21)

by including increasingly higher-order and long-range terms. However, the
number of variational parameters quickly proliferates requiring a truncation
in the order of terms. Moreover, the inclusion of these new terms may only
have marginal benefits since they are not able to capture specific high-order
terms like those contained in the many-body correlator. Thus, a potentially
powerful alternative approach for generalising the classic bosonic ansatzes is
to use NQS.

##  4 Neural-network Quantum States

A general NQS can leverage the differing strengths of the numerous network
architectures employed in machine learning. Here we will exclusively focus on
the simplest the RBM form.

###  4.1 Restricted Boltzmann Machines

A Restricted Boltzmann Machine [88, 89] is a classical probabilistic model
comprising a layer of $L$ visible units specified by
${\bm{v}}=\\{v_{1},v_{2},\dots,v_{L}\\}$ with binary variables
$v_{i}\in\mathbb{N}_{1}=\\{0,1\\}$ representing the input data, and a layer of
$M$ hidden units specified by ${\bm{h}}=\\{h_{1},h_{2},\dots,h_{M}\\}$ binary
variables $h_{\mu}\in\mathbb{N}_{1}$ which are marginalised internal
variables. The Boltzmann in RBM refers to the parameterisation of the model in
terms of an effective energy function

$\mathcal{E}_{{\bm{\lambda}}}({\bm{v}},{\bm{h}})=-\sum_{i=1}^{L}a_{i}v_{i}-\sum_{\mu=1}^{M}b_{\mu}h_{\mu}-\sum_{\mu=1}^{M}\sum_{i=1}^{L}w_{\mu
i}h_{\mu}v_{i},$ (22)

with $n_{\textrm{params}}=M+L+ML$ parameters
${\bm{\lambda}}=\\{{\bm{a}},{\bm{b}},{\bm{w}}\\}$ including $L$ visible biases
in vector ${\bm{a}}$, $M$ hidden biases in vector ${\bm{b}}$ and $M\times L$
interaction weights in matrix ${\bm{w}}$. It is also useful to denote
${\bm{w}}_{\mu}$ as the vector of $L$ weights for a specific hidden unit
$\mu$. The restricted in RBM refers to bipartite structure of the interactions
in Eq. (22) which are fully connected between the visible and hidden variables
only, as depicted in Fig. 2. The (unnormalised) probability of a given visible
configuration ${\bm{v}}$ is then given by

$P_{{\bm{\lambda}}}({\bm{v}})=\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left[-\mathcal{E}_{\bm{\lambda}}({\bm{v}},{\bm{h}})\right].$
(23)

Although the joint probability of ${\bm{v}}$ and ${\bm{h}}$ is Boltzmann, once
the hidden units are marginalised $P_{{\bm{\lambda}}}({\bm{v}})$ is capable of
describing much more complex visible unit distributions with increasing $M$.
The complexity of an RBM is often quantified by its hidden unit density
$\alpha=M/L$.

Figure 2: The bipartite graph of interaction weights $\bm{w}$ between hidden
and visible units in an RBM with edges shown as solid arcs. The biases
$\bm{a}$ and $\bm{b}$ on each unit are depicted as self-loop edges.

The generalisation of RBMs to a variational ansatz for quantum systems is most
directly accomplished by considering an $L$ site hard-core boson system. The
visible configuration ${\bm{v}}$ then specifies the real-space occupation
number basis $\left|\,\bm{v}\right\rangle$ via hard-core boson number
operators $\hat{v}_{i}$ whose binary eigenvalues in this basis are
$\hat{v}_{i}\left|\,\bm{v}\right\rangle=v_{i}\left|\,\bm{v}\right\rangle$. The
RBM correlator is then constructed [29] by elevating the classical visible
variables $v_{i}$ to operators $\hat{v}_{i}$ as

$\displaystyle\qquad\left|\,\Psi_{\rm RBM}({{\bm{\lambda}}})\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left(\sum_{i=1}^{L}a_{i}\hat{v}_{i}+\sum_{\mu=1}^{M}b_{\mu}h_{\mu}+\sum_{\mu=1}^{M}\sum_{i=1}^{L}w_{\mu
i}h_{\mu}\hat{v}_{i}\right)\left|\,\Phi_{+}\right\rangle,$ $\displaystyle=$
$\displaystyle\exp\left(\sum_{i=1}^{L}a_{i}\hat{v}_{i}\right)\prod_{\mu=1}^{M}\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}}_{\mu})\left|\,\Phi_{+}\right\rangle,$

where we have identified each hidden unit’s factor as
$\hat{\Upsilon}_{i}(b_{\mu},{\bm{w}}_{\mu})$ and applied the full correlator
to the uniform reference state

$\left|\,\Phi_{+}\right\rangle=\sum_{{\bm{v}}\in\mathbb{N}_{1}^{L}}\left|\,\bm{v}\right\rangle.$
(25)

If a restriction to a fixed number $N$ of hard-core bosons is required then
$\left|\,\Phi_{+}\right\rangle$ can be projected into this sector on-the-fly
within Monte Carlo sampling. Finally, in the context of modelling quantum
states the parameters $\bm{\lambda}$ are permitted to be complex-valued. A
restriction to real-valued amplitudes $\left\langle\bm{v}|\Psi_{\rm
RBM}({{\bm{\lambda}}})\right\rangle$ on all configurations requires the
visible biases $\bm{a}$ be real, but the hidden bias $\bm{b}$ and weights
$\bm{w}$ can be either real or imaginary. Positive-definite amplitudes are
guaranteed by restricting all parameters $\bm{\lambda}$ to be real.

###  4.2 Representation properties

The terms and operators that appear in the hard-core RBM correlator have
important representational consequences. In particular they involve
$\hat{v}_{i}$ (explicitly) and $\mathbbm{1}$ (trivially), which together form
a complete basis of local diagonal operators at site $i$. This gives rise to
two related properties. First, the presence of the visible bias term
$\sum_{i=1}^{L}a_{i}\hat{v}_{i}$ ensures that the RBM correlator can describe
an arbitrary product correlator

$\hat{C}_{\rm
prod}({\bm{c}})=\prod_{i=1}^{L}\left[\mathbbm{1}+(c_{i}-1)\hat{v}_{i}\right],$
(26)

without using any hidden units by simply setting the visible biases to

$a_{i}=\log(c_{i}-1).$ (27)

This product representability property means that, sensibly, hidden units are
only required in an RBM to capture non-local correlations. Second, the hidden
bias $\sum_{\mu=1}^{M}b_{\mu}h_{\mu}$ and interaction terms
$\sum_{\mu=1}^{M}\sum_{i=1}^{L}w_{\mu i}h_{\mu}\hat{v}_{i}$ ensure that the
states an RBM ansatz can describe are in fact independent on the choice of
local diagonal operator to couple to. Specifically, if in place of
$\hat{v}_{i}$ we instead constructed the RBM correlator in Eq. (4.1) by
coupling to a different local diagonal operator

$\hat{s}_{i}=\chi\mathbbm{1}+(\kappa-\chi)\hat{v}_{i},$ (28)

with eigenvalues $\chi\neq\kappa$, then the new RBM parameters follow from the
original ones directly as

$a_{i}\mapsto\frac{a_{i}}{\kappa-\chi},\quad b_{\mu}\mapsto
b_{\mu}-\frac{\chi}{\kappa-\chi}\sum_{j=1}^{L}w_{\mu i},\quad w_{\mu
i}\mapsto\frac{w_{\mu i}}{\kappa-\chi}.$ (29)

A common alternative choice, for example, is for $\hat{s}_{j}$ to be a spin-
like operator with eigenvalues $\chi=1$ and $\kappa=-1$ 111Note exactly the
same reasoning applies to the classical hidden variables $h_{\mu}$. They can
also be taken as having any pair of values $h_{i}\in\\{\chi,\kappa\\}$ with a
simple mapping of the RBM parameters.. This labelling freedom property
together with product representability will shortly guide our attempts to
generalising the RBM ansatz to soft-core boson systems.

A powerful feature of the RBM ansatz its ability to be systematically improved
by increasing the number of hidden units $M$. Indeed, if we let $M\sim 2^{L}$
then $\left|\,\Psi_{\rm RBM}({{\bm{\lambda}}})\right\rangle$ can describe
exactly any arbitrary state of hard-core bosons [56, 59], demonstrating that
the RBM ansatz is exhaustive. However, for practical numerical calculations we
instead expect a scaling $M\sim\textrm{poly}(L)$ so that
$\left\langle{\bm{v}}|\Psi_{\rm RBM}({\bm{\lambda}})\right\rangle$ can be
sampled efficiently in VMC.

Since their introduction as a wavefunction ansatz in Ref. [29], considerable
work has gone into formulating efficient exact NQS representations of
spin-$\textstyle\frac{1}{2}$ quantum states as an RBM. This has so far
included cluster [90], graph [30, 59] and hypergraph [57], Jastrow [58, 91,
59], stabiliser [61, 92, 93, 57, 94] and XS-stabiliser [57] states. Many of
these states can be described [61] with just $M=L$ hidden units ($\alpha=1$)
making them very compact.

###  4.3 Capturing exactly Jastrow and many-body correlators

To connect the RBM ansatz to the classic bosonic ansatzes discussed earlier in
Sec. 3 it is instructive to rewrite it in the form $\left|\,\Psi_{\rm
RBM}({{\bm{\lambda}}})\right\rangle=\exp[-E({\bm{\lambda}},\\{\hat{v_{i}}\\})]\left|\,\Phi_{+}\right\rangle$
with a pseudo-energy function

$\displaystyle E({\bm{\lambda}},\\{\hat{v_{i}}\\})$ $\displaystyle=$
$\displaystyle-\sum_{i=1}^{L}a_{i}\hat{v}_{i}-\sum_{\mu=1}^{M}\log\left[1+\exp(\hat{\theta}_{\mu})\right],$
(30) $\displaystyle=$
$\displaystyle-\sum_{i=1}^{L}a_{i}\hat{v}_{i}-\sum_{\mu=1}^{M}\left(\frac{1}{2}\hat{\theta}_{\mu}+\frac{1}{8}\hat{\theta}_{\mu}^{2}-\frac{1}{192}\hat{\theta}_{\mu}^{4}+\dots\right),$
(31)

where each hidden unit contributes a series expansion of a physical operator

$\displaystyle\hat{\theta}_{\mu}=b_{\mu}+\sum_{i=1}^{L}w_{\mu i}\hat{v}_{i}.$
(32)

From this we see that $E({\bm{\lambda}},\\{\hat{v_{i}}\\})$ naturally contains
high-order interactions, analogous to those in Eq. (21), albeit with a
complicated parameterisation in terms of ${\bm{b}}$ and ${\bm{w}}$.

Some simple but useful constructions are possible if we permit diverging
parameters [61]. To perfectly correlate a particular hidden variable $h_{\mu}$
with a given visible unit’s occupation $v_{i}$ its bias and weight are set
accordingly as

$\displaystyle b_{\mu}=-\mathcal{S},\quad w_{\mu
j}=\left\\{\begin{array}[]{cc}2\mathcal{S}&j=i\\\ \nu_{ij}&j\neq
i\end{array}\right.,$

where we will formally take the limit $\mathcal{S}\rightarrow\infty$. The
action alone of the $\mu$th hidden unit on an arbitrary basis state
$\left|\,\bm{v}\right\rangle$ can then be expressed as

$\displaystyle\qquad\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}}_{\mu})\left|\,\bm{v}\right\rangle$
$\displaystyle=$ $\displaystyle\exp\left(\mathcal{S}v_{i}+\sum_{j\neq
i}^{L}\nu_{ij}v_{i}v_{j}\right)\times$ (34)
$\displaystyle\qquad\left(1+\exp\left[-\mathcal{S}\left\\{1+\sum_{j\neq
i}^{L}\frac{\nu_{ij}}{\mathcal{S}}v_{j}-2\sum_{j\neq
i}^{L}\frac{\nu_{ij}}{\mathcal{S}}v_{i}v_{j}\right\\}\right)\right].$

After shifting the bias of the $i$th visible unit as $a_{i}\mapsto
a_{i}-\mathcal{S}$ we find

$\displaystyle\qquad\qquad\lim_{\mathcal{S}\rightarrow\infty}\exp\left[(a_{i}-\mathcal{S})\hat{v}_{i}\right]\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}_{\mu}})$
$\displaystyle=$
$\displaystyle\exp\left(a_{i}\hat{v}_{i}\right)\exp\left(\sum_{j\neq
i}^{L}\nu_{ij}\hat{v}_{i}\hat{v}_{j}\right).$ (35)

Thus, a single hidden unit can capture exactly arbitrary Jastrow density-
density effective energy terms between one fixed visible unit $i$ and the rest
of the system, as depicted in Fig. 3(a).

Figure 3: A depiction of the weights and biases required to construct for
hard-core bosons a (a) Jastrow correlator as in Eq. (35) and (b) many-body
correlator as in Eq. (37). In both cases $\mathcal{S}\gg 1$.

The construction leading to Eq. (35) is limited to two-body interactions since
it exploits those built into the RBM formalism. Another related construction
can generate many-body interactions by fixing the $\mu$th hidden unit’s bias
and its weights with a set of visible units $\Omega$ as

$\displaystyle b_{\mu}=-|\Omega|\mathcal{S}+\log(e^{-\xi}-1),\quad w_{\mu
i}=\left\\{\begin{array}[]{cl}\mathcal{S}&i\in\Omega\\\ 0&{\rm
otherwise}\end{array}\right.,$

where $\xi>0$ is a real parameter222In this construction our hidden bias
$b_{\mu}$ is complex. However, we can keep it a single real parameter within
variational calculations by introducing a fixed imaginary part $b_{\mu}\mapsto
b_{\mu}+{\rm i}\pi$.. In this case we find

$\displaystyle\qquad\lim_{\mathcal{S}\rightarrow\infty}\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}}_{\mu})$
$\displaystyle=$
$\displaystyle\left[1+(e^{-\xi}-1)\lim_{\mathcal{S}\rightarrow\infty}\exp\left(\mathcal{S}\sum_{i\in\Omega}\hat{v}_{i}-|\Omega|\mathcal{S}\right)\right],$
(37) $\displaystyle=$
$\displaystyle\exp\left(-\xi\prod_{i\in\Omega}\hat{v}_{i}\right).$

Thus, a single hidden unit can also generate specific high-order interaction
terms involving any set $2\leq|\Omega|\leq L$ of visible units, as depicted in
Fig. 3(b). In practise the constructions Eq. (35) and Eq. (37) can both be
realised to good accuracy with finite $\mathcal{S}\sim O(10)$ since their
errors scale as $\exp(-\mathcal{S})$.

###  4.4 Imposing translational invariance

For a periodic $\sqrt{L}\times\sqrt{L}$ 2D lattice there are $L$ distinct
combinations of translations along the $x$\- and $y$-axes labelled as
$s\in\\{1,2,\dots,L\\}$ which implement a site index map $i\mapsto s(i)$ and
define a translation operator as
$\hat{T}_{s}\hat{v}_{i}\hat{T}_{s}^{\dagger}=\hat{v}_{s(i)}$. To enforce that
the RBM ansatz is translationally invariant we require that

$\displaystyle\hat{T}_{s}\left|\,\Psi_{\rm RBM}({\bm{\lambda}})\right\rangle$
$\displaystyle=\left|\,\Psi_{\rm
RBM}({\bm{\lambda}})\right\rangle,\quad\forall~{}s.$ (38)

Since the reference state $\left|\,\Phi_{+}\right\rangle$ obeys these
constraints already an immediate consequence is that the physical biases are
independent of the site index, so $a_{i}\mapsto a$, leaving

$\displaystyle\hat{T}_{s}\left|\,\Psi_{\rm
RBM}({{\bm{\lambda}}})\right\rangle$ $\displaystyle=$
$\displaystyle\exp\left(a\sum_{i=1}^{L}\hat{v}_{i}\right)\prod_{\mu=1}^{M}\hat{T}_{s}\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}}_{\mu})\hat{T}^{\dagger}_{s}\,\left|\,\Phi_{+}\right\rangle,$
$\displaystyle=$
$\displaystyle\exp\left(a\sum_{i=1}^{L}\hat{v}_{i}\right)\prod_{\mu=1}^{M}\hat{\Upsilon}_{\mu}(b_{\mu},{\bf
T}_{s}{\bm{w}}_{\mu})\,\left|\,\Phi_{+}\right\rangle,$

where ${\bf T}_{s}{\bm{w}}_{\mu}$ is the vector of $L$ weights $w_{\mu s(i)}$.
Working in a fixed number sector $N=\sum_{i=1}^{L}v_{i}$ renders the
contribution of the visible bias an irrelevant overall constant.
Translationally invariant states within the RBM ansatz thus arise from subsets
of hidden unit correlators $\hat{\Upsilon}_{\mu}(b_{\mu},{\bm{w}}_{\mu})$
mapping between themselves under translation. A non-trivial class of such
states are constructed from $\alpha$ arbitrary hidden units with weights
$w_{\eta i}$ and biases $b_{\eta}$ for $\eta\in\\{1,2,\dots,\alpha\\}$. For
each $\eta$ there are a further $L-1$ hidden units with an identical bias
$b_{\eta}$ and weights which are translates $w_{\eta s(i)}$ for each $s$. In
this case the RBM ansatz is

$\displaystyle\left|\,\Psi_{\rm RBM}({{\bm{\lambda}}})\right\rangle$
$\displaystyle=$
$\displaystyle\exp\left(a\sum_{i=1}^{L}\hat{v}_{i}\right)\prod_{\eta=1}^{\alpha}\prod_{s^{\prime}=1}^{L}\hat{\Upsilon}_{\eta
s^{\prime}}(b_{\eta},{\bf
T}_{s^{\prime}}{\bm{w}}_{\eta})\,\left|\,\Phi_{+}\right\rangle,$

defined by $n_{\rm params}=1+\alpha+\alpha L$ variational parameters, but
formally has $M=\alpha L$ hidden units in total. Owing to the group property
of translations we now have [29]

$\displaystyle\quad\hat{T}_{s}\left|\,\Psi_{\rm
RBM}({{\bm{\lambda}}})\right\rangle$ $\displaystyle=$
$\displaystyle\exp\left(a\sum_{i=1}^{L}\hat{v}_{i}\right)\prod_{\eta=1}^{\alpha}\prod_{s^{\prime}=1}^{L}\hat{\Upsilon}_{\eta
s^{\prime}}(b_{\eta},{\bf T}_{s}{\bf
T}_{s^{\prime}}{\bm{w}}_{\eta})\,\left|\,\Phi_{+}\right\rangle=\left|\,\Psi_{\rm
RBM}({{\bm{\lambda}}})\right\rangle.$

As a result, our construction of the Jastrow correlator in Eq. (35) and many-
body correlator in Eq. (37) each only require $\alpha=1$ hidden unit once
translational invariance is explicitly enforced in this way. The generality of
this method of imposing translational symmetry means it applies to all
variants of the RBM ansatz we will consider.

##  5 Extension to bosonic systems

The generalisation of the $N$-particle projection of the hard-core boson
reference state $\left|\,\Phi_{+}\right\rangle$ to soft-core bosons is the BEC
reference state $\left|\,\Phi_{0}\right\rangle$ introduced in Sec. 3. Note
that the operators $\hat{v}_{i}$ involved in the RBM construction for hard-
core bosons act both as a local number operator and a projector on to the
local occupied configuration $\left|\,1\right\rangle$ of a site $i$. In the
soft-core boson case either of these properties provide possible routes for
generalisation.

###  5.1 Naive number operator replacement

The first and simplest generalisation is to directly replace each hard-core
number operator $\hat{v}_{j}$ with a soft-core boson number operator
$\hat{n}_{j}$ yielding an RBM-inspired correlator

$\hat{C}_{0}=\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left(\sum^{L}_{i=1}\tilde{a}_{i}\hat{n}_{i}+\sum_{\mu=1}^{M}\tilde{b}_{\mu}h_{\mu}+\sum_{\mu=1}^{M}\sum_{\mu=1}^{L}\tilde{w}_{\mu
i}h_{\mu}\hat{n}_{i}\right),$ (39)

with exactly the same number of parameters. We tested the translationally
invariant version of this ansatz with $\alpha=1$ for both real parameters and
imaginary hidden bias and weights. In the case of real parameters the
stochastic optimisations did not converge, suggesting an expressivity failure.
This observation is in agreement with those in Ref. [38], where similar
difficulties were encountered. In the case of some imaginary parameters
successful minimisation was possible, but the performance for all but $U/t\ll
1$ was generically poor due to the lack of positivity in the variational
state. The origin of this behaviour can be linked to this ansatz not having
product representability. Indeed, there is no local mechanism for this ansatz
to suppress multiple occupation of a site as implemented by the Gutzwiller
correlator. Instead of using $\hat{C}_{0}$ we will construct improved variants
of it by considering the alternative route to soft-core bosons via projection
operators.

###  5.2 Projector expansion (NQS-OH)

A more comprehensive generalisation of RBM to soft-core bosons is to introduce
independent visible biases $a^{[x]}_{i}$ and interaction weights $w^{[x]}_{\mu
i}$ for each projector $\hat{\mathbbm{P}}^{[x]}_{i}$ on to the local
occupation states $\left|\,x\right\rangle$ of the sites $i$. By exploiting
local completeness one of the projector can be chosen to be omitted from this
expansion. In the hard-core boson case this was the projector on to the empty
state $\left|\,0\right\rangle$. Following this here we arrive at a generalised
RBM correlator

$\hat{C}_{\rm
OH}=\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left(\sum_{x=1}^{B}\sum^{L}_{i=1}a^{[x]}_{i}\hat{\mathbbm{P}}^{[x]}_{i}+\sum_{\mu=1}^{M}b_{\mu}h_{\mu}+\sum_{x=1}^{B}\sum_{\mu=1}^{M}\sum_{i=1}^{L}w^{[x]}_{\mu
i}h_{\mu}\hat{\mathbbm{P}}^{[x]}_{i}\right).$ (40)

We now have $n_{\rm params}=BL+M+BML$ variational parameters. This ansatz is
completely equivalent333One technical difference here is that we have
exploited the local completeness relation to reduce the parameter count by
$L+ML$ compared to conventional one-hot encoding. to the one-hot or unary RBM
encoding scheme commonly used in machine learning to handle discrete multi-
valued visible variables [95]. Crucially, since it is built on a locally
complete set of projectors this ansatz immediately inherits the product
representability and labelling freedom properties of the hard-core boson RBM.
Moreover, the ansatz reduces back to the hard-core boson RBM in the limit
$B=1$ where $\hat{\mathbbm{P}}^{[1]}_{i}\mapsto\hat{v}_{i}$.

Figure 4: The one-hot RBM encoding for a bosonic system where each hidden unit
couples independently to $B$ projection operators
$\hat{\mathbbm{P}}^{[x]}_{i}$ on the local state $\left|\,x\right\rangle$ of a
each physical site $i$.

As it stands the ansatz couples operators with binary eigenvalues to binary-
valued hidden units. However, labelling freedom permits us to change our basis
of diagonal local operators from
$\\{\mathbbm{1},\mathbbm{P}_{i}^{[1]},\mathbbm{P}_{i}^{[2]},\dots,\mathbbm{P}_{i}^{[B]}\\}$.
A natural alternative choice of basis is
$\\{\mathbbm{1},\hat{n}_{i},\hat{n}_{i}^{2},\dots,\hat{n}_{i}^{B}\\}$ giving
an equivalent correlator

$\hat{C}_{\rm
OH}=\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left(\sum_{x=1}^{B}\sum^{L}_{i=1}\tilde{a}^{[x]}_{i}\hat{n}^{x}_{i}+\sum_{\mu=1}^{M}b_{\mu}h_{\mu}+\sum_{x=1}^{B}\sum_{\mu=1}^{M}\sum_{i=1}^{L}\tilde{w}^{[x]}_{\mu
i}h_{\mu}\hat{n}^{x}_{i}\right).$ (41)

Note that by excluding $\hat{\mathbbm{P}}^{[0]}_{i}$ from our basis the hidden
bias is unchanged in this transformation. Inverting Eq. (5) the linear
relation between the two sets of visible bias and weights can be extracted,
for example with $B=2$ we have

$\displaystyle\tilde{\bm{a}}^{[1]}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(4{\bm{a}}^{[1]}-{\bm{a}}^{[2]}\right),\quad~{}\tilde{\bm{a}}^{[2]}=\frac{1}{2}\left({\bm{a}}^{[2]}-2{\bm{a}}^{[1]}\right),$
(42) $\displaystyle\tilde{\bm{w}}^{[1]}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(4{\bm{w}}^{[1]}-{\bm{w}}^{[2]}\right),\quad\tilde{\bm{w}}^{[2]}=\frac{1}{2}\left({\bm{w}}^{[2]}-2{\bm{w}}^{[1]}\right).$
(43)

The spin-1 version of Eq. (41) was introduced and exploited already in Ref.
[60].

Later in Sec. 6 we will investigate the performance of one-hot encoding.
Formally, it is the most powerful ansatz we consider in this work. However, by
fully discriminating every local configuration $\left|\,x\right\rangle$ of a
site $n_{\rm params}$ scales with $B$ making its computational cost
prohibitive. This also compounds the general criticism of RBMs that the
physical interpretation of their parameters is often opaque. For this reason
it is insightful to consider versions of both Eq. (40) and Eq. (41) in which
numerous terms are removed and their equivalence is broken. A guiding
principle in proposing these variants will be the ability to efficiently
reproduce Gutzwiller, Jastrow and many-body correlators from Sec. 3 while
offering systematic extensions beyond them.

###  5.3 Holon-doublon variant (NQS-HD)

The physics of the BHM in the $U/t>1$ limit at unit-filling indicates that the
dominant contributions to the ground state arise from holon-doublon
fluctuations about the atomic limit $\left|\,1,1,1,\dots,1\right\rangle$ Mott
state. It is thus sensible to use local completeness to remove the singlon
projector $\hat{S}_{i}$ from the basis and then to truncate Eq. (40) to the
biases and couplings to just the holon $\hat{H}_{i}$ and doublon $\hat{D}_{i}$
projectors. This gives a new RBM correlator

$\displaystyle\hat{C}_{\rm HD}$ $\displaystyle=$
$\displaystyle\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left[\sum^{L}_{i=1}\left(a^{[0]}_{i}\hat{H}_{i}+a^{[2]}_{i}\hat{D}_{i}+a^{[m]}_{i}\hat{M}_{i}\right)+\sum_{\mu=1}^{M}b_{\mu}h_{\mu}\right.$
(44)
$\displaystyle\qquad\qquad\qquad+\left.\sum_{\mu=1}^{M}\sum_{i=1}^{L}\left(w^{[0]}_{\mu
i}h_{\mu}\hat{H}_{i}+w^{[1]}_{\mu j}h_{\mu}D_{i}\right)\right].$

This ansatz has $n_{\rm params}=3L+M+2ML$. By construction the ansatz only
describes correlations between holons and doublons. The inclusion of biases
${\bm{a}}^{[m]}$ associated to the multiplon operator ensures the ansatz can
uniformly suppress configurations $\bm{n}$ with $n_{i}\geq 3$.

The holon-doublon variant is less expensive than one-hot encoding, but can
still capture all the classic variational ansatzes. Within the subspace
$\mathcal{F}_{N,3}$ this ansatz can exactly reproduce the Gutzwiller
correlator with no hidden units by using visible biases

$\displaystyle
a^{[0]}_{i}=a^{[2]}_{i}=-\mbox{$\textstyle\frac{1}{2}$}g\quad{\rm and}\quad
a^{[m]}_{i}=-2g.$ (45)

Similarly for configurations in $\mathcal{F}_{N,2}$ it can exactly reproduce
the Jastrow density-density correlator following the construction introduced
in Eq. (35), and the many-body correlator by modifying the construction in Eq.
(37). Details are given in A. However, in both cases parameters are spread
across holon and doublon weights for a pair of hidden units. One might expect
that there are even more economical variants.

###  5.4 Number variant - (NQS-A)

We can equally consider truncating one-hot encoding in the basis of number
operators. Indeed, in the context of Eq. (41) we can now appreciate the naive
ansatz $\hat{C}_{0}$ in Eq. (39) as being an extreme truncation to just
$\hat{n}_{i}$. Its inability to provide local occupation suppression is easily
overcome by retaining the $\hat{n}_{i}^{2}$ bias as

$\hat{C}_{\rm
A}=\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left[\sum_{i=1}^{L}\left(\tilde{a}^{[1]}_{i}\hat{n}_{i}+\tilde{a}^{[2]}_{i}\hat{n}_{i}^{2}\right)+\sum_{\mu=1}^{M}\left(b_{\mu}h_{\mu}+\sum_{\mu=1}^{L}\tilde{w}^{[1]}_{\mu
i}h_{\mu}\hat{n}_{i}\right)\right].$ (46)

This ansatz has $n_{\rm params}=2L+M+ML$. The Gutzwiller correlator is exactly
reproduced by this variant in the limit of no hidden units using
$\tilde{a}_{i}^{[1]}=g$ and
$\tilde{a}_{i}^{[2]}=-\mbox{$\textstyle\frac{1}{2}$}g$. It is capable of
generating two-body and multi-body density-density interaction terms but has
no simple exact parameterisation of the Jastrow or the many-body correlator.

###  5.5 Hidden unit expansion - (NQS-B)

The expansion in the number basis couples binary-valued hidden units to
operators which are multi-valued. A novel modification is to expand the hidden
unit to be identically multi-valued, as $h_{\mu}\in\mathbb{N}_{B}$, and
exploit it by including non-linear terms involving $h_{\mu}$. The simplest is
a quadratic hidden bias leading to a generalisation of Eq. (46) as

$\hat{C}_{\rm
B}=\sum_{{\bm{h}}\in\mathbb{N}_{B}^{M}}\exp\left[\sum_{i=1}^{N}\left(\tilde{a}^{[1]}_{i}\hat{n}_{i}+\tilde{a}^{[2]}_{i}\hat{n}^{2}_{i}\right)+\sum_{\mu=1}^{M}\left(b^{[1]}_{\mu}h_{\mu}+b^{[2]}_{\mu}h_{\mu}^{2}+\sum_{i=1}^{L}\tilde{w}^{[1]}_{\mu
i}h_{\mu}\hat{n}_{i}\right)\right],$ (47)

which has only a modest increase in parameters to $n_{\rm params}=2L+2M+ML$.
The inclusion of quadratic biases make this ansatz very similar in form to
Gaussian RBMs for continuous visible and hidden variables [96]. An immediate
consequence of this modification is that the bosonic occupation of any site
can now be perfectly correlated with a multi-valued hidden variable. This
allows one hidden unit to exactly describe Jastrow correlator of that site
with the rest, as outlined in B. Nonetheless this ansatz is still unable to
easily capture the many-body correlator motivating one final variant.

###  5.6 Quadratic number variant - (NQS-C)

The next logical step in generalising the ansatz in Eq. (47) is to include
weights that couple hidden units to $\hat{n}_{i}^{2}$ as

$\displaystyle\hat{C}_{\rm C}$
$\displaystyle=\sum_{{\bm{h}}\in\mathbb{N}_{B}^{M}}\exp\left[\sum_{i=1}^{N}\left(\tilde{a}^{[1]}_{i}\hat{n}_{i}+\tilde{a}^{[2]}_{i}\hat{n}^{2}_{i}\right)+\sum_{\mu=1}^{M}\left(b^{[1]}_{\mu}h_{\mu}+b^{[2]}_{\mu}h_{\mu}^{2}\right)\right.$
$\displaystyle\qquad\qquad\qquad\quad+\left.\sum_{\mu=1}^{M}\sum_{i=1}^{L}\left(\tilde{w}^{[1]}_{\mu
i}h_{\mu}\hat{n}_{i}+\tilde{w}^{[2]}_{\mu
i}h_{\mu}\hat{n}^{2}_{i}\right)\right].$ (48)

This ansatz now has $n_{\rm params}=2L+2M+2ML$. It continues to capture
Gutzwiller with no hidden units and translationally invariant Jastrow with
just one hidden unit. The additional quadratic interaction terms now allow
this ansatz to directly discriminate holons and doublons and thus capture
exactly the many-body correlator in the subspace $\mathcal{F}_{N,2}$. However,
like the holon-doublon variant, it requires a pair of hidden units444The
ansatz does not fully leverage the extra bandwidth of the hidden unit to
reduce this to a single unit. This requires $h^{2}_{\mu}\hat{n}_{i}$ and
$h^{2}_{\mu}\hat{n}^{2}_{i}$ interactions proliferating the number of
variational parameters. following a modified version of the construction in
Eq. (37), as described in C.

Wavefunction | Shorthand | Parameter sets | $n_{\rm params}$
---|---|---|---
Jastrow Many-body | J-MB | $\nu_{ij},\xi$ | $\propto L$
One-hot | OH | $a^{[x]},{\bm{b}},{\bm{w}}^{[x]},~{}x=1,2,\dots,B$ | $B+\alpha+\alpha BL$
Holon-doublon | HD | $a^{[0]},a^{[2]},a^{[m]},{\bm{b}},{\bm{w}}^{[0]},{\bm{w}}^{[2]}$ | $3+\alpha+2\alpha L$
Quadratic number bias | A | $\tilde{a}^{[1]},\tilde{a}^{[2]},{\bm{b}},\tilde{\bm{w}}^{[1]}$ | $2+\alpha+\alpha L$
Hidden quadratic bias | B | $\tilde{a}^{[1]},\tilde{a}^{[2]},{\bm{b}}^{[1]},{\bm{b}}^{[2]},\tilde{\bm{w}}^{[1]}$ | $2+2\alpha+\alpha L$
Quadratic interaction | C | $\tilde{a}^{[1]},\tilde{a}^{[2]},{\bm{b}}^{[1]},{\bm{b}}^{[2]},\tilde{\bm{w}}^{[1]},\tilde{\bm{w}}^{[2]}$ | $2+2\alpha+2\alpha L$
Table 1: A table of the translationally invariant variational wavefunctions
presented in the results section, including shorthands, types and number of
parameters in the ansatz. Recall that $L$ is the total number of lattice
sites, $B$ is the maximum bosonic occupation of a site, and $\alpha$ is the
hidden unit density of the RBM.

##  6 Results

To explore the performance of our proposed RBM variants we focus on four
regimes of the BHM: the weakly interacting SF at $U/t=1$; the vicinity of the
transition point $U/t=16$ predicted by quantum Monte Carlo [72] and strong
coupling expansions [73]; the vicinity of the transition point $U/t=23$
predicted by mean-field theory [24]; and the strongly interacting MI at
$U/t=32$. These points were highlighted earlier in Fig. 1(a). We gauge the
effectiveness of an ansatz by comparing their groundstate energy density
$\epsilon$ against the number of variational parameters $n_{\textrm{params}}$
which controls the computational cost of the optimisation. For all calculation
we consider a $10\times 10$ lattice with periodic boundary conditions.

We analyse two scenarios. The first case uses the standard BEC reference state
$\left|\,\Phi_{0}\right\rangle$ for all NQS ansatzes, mirroring that used by
the classic bosonic ansatz. Given the lack of structure in
$\left|\,\Phi_{0}\right\rangle$ this allows the broad abilities of our
proposed NQS ansatzes to describe the BHM regimes to be revealed and compared
against the Jastrow and J-MB ansatzes. The second case uses the best classic
bosonic ansatz, namely a pre-optimised J-MB state $\left|\,\Psi_{\rm
J-MB}\right\rangle$ [85, 87], as the reference state for all NQS ansatzes.

Figure 5: Variational energy density $\epsilon$ of the for BHM on a $10\times
10$ lattice obtained by the NQS-C, -HD and -OH ansatzes applied to a BEC
reference $\left|\,\Phi_{0}\right\rangle$ for (a) $U/t=1$, (b) $U/t=16$, (c)
$U/t=23$ and (d) $U/t=32$. For NQS-C and -HD the two data points are for
$\alpha=1$ and $\alpha=2$ and are shown connected with solid lines to guide
the eye. Key comparators are the energy densities for Jastrow (dashed line)
and Jastrow + many-body correlator (dashed-dotted line) ansatzes. The region
between these results is also shaded to guide the eye.

###  6.1 Case 1: Bose-condensate (BEC) reference state

We start by considering NQS-C, -HD and -OH shown in Fig. 5. In principle these
are the most flexible but also most expensive ansatzes proposed. Given the
specialisation of NQS-HD to holon-doublon correlations, dominant for $U/t\gg
1$, its comparatively poor performance in Fig. 5(a) for $U/t=1$ is not
surprising. For increasing interactions both NQS-HD and -C leverage the
potential for $\alpha>1$ to improve their description. For the strongly
interacting SF regime $U/t=16$, shown in Fig. 5(b), NQS-C produces a 2%
improvement over the J-MB ansatz with $\alpha=2$. Beyond this in the MI regime
all NQS ansatzes struggle to match J-MB ansatzes or even improve on plain
Jastrow. The most expensive $\alpha=1$ NQS-OH ansatz has the order of
magnitude increase in $n_{\textrm{params}}$ compared to $\left|\,\Psi_{\rm
J-MB}\right\rangle$, but this fails to translate into a consistently more
accurate description in any regime due to the complexity of its numerical
optimisation.

Figure 6: Variational energy density $\epsilon$ of the for BHM on a $10\times
10$ lattice obtained by the NQS-A and -B ansatzes applied to a BEC reference
$\left|\,\Phi_{0}\right\rangle$ for (a) $U/t=1$, (b) $U/t=16$, (c) $U/t=23$
and (d) $U/t=32$. For both ansatzes data for $\alpha=1,2,3$ and $4$ is shown
connected with solid lines to guide the eye. Key comparators are the energy
densities for Jastrow (dashed line) and Jastrow + many-body correlator
(dashed-dotted line) ansatzes. The region between these results is also shaded
to guide the eye.

Moving to the ansatzes NQS-A and -B in Fig. 6 we see a similar picture of
improvements in the SF regime and difficulties in the MI regime. Although
these simpler ansatz show improvement with increasing $\alpha$, and even
$\alpha=4$ can be reached, the optimisation is often found to display shallow
plateauing behaviour. Nonetheless, these results confirm that all proposed NQS
correlators can reshape the structureless BEC reference
$\left|\,\Phi_{0}\right\rangle$ into a reasonable ground state in proximity to
the J-MB. The difficulty in improving on the $\left|\,\Psi_{\rm
J-MB}\right\rangle$ highlights how accurate and compact this ansatz is.

Interestingly, numerical optimisation of NQS-B with $\alpha=1$ can be found to
converge to a solution closely resembling the exact NQS description of
Jastrow. Specifically, each hidden unit has a dominant weight to a distinct
visible unit, while the weights to all other sites mimic the Jastrow pseudo-
potential. In Fig. 7(a) the pseudo-potential is shown for $U/t=1$ along with
the NQS-B $\alpha=1$ weights in Fig. 7(b), demonstrating they are very
similar, but not identical. Unlike the exact construction, the dominant weight
is not permitted to diverge during the optimisation. However, this deviation
is actually seen to improve the description, as displayed in Fig. 6(a) where
$\epsilon$ is marginally improved over J-MB. The weights for additional hidden
units with NQS-B no longer display any dominant coupling. This suggests NQS-B
first forms a Jastrow-like state and then modifies it with successive hidden
units. Motivated by this observation we consider now in detail the second case
where an NQS correlator is applied to a structured reference state.

Figure 7: (a) For $U/t=1$ the pseudo-potential $\nu_{ij}$ of a Jastrow + many-
body state shown for site $i=(5,5)$ over sites $j=(x,y)$ over a $10\times 10$
lattice. Translationally invariance means that changing $i$ simply shifts the
pattern around the lattice. (b) The corresponding weights $w^{[1]}_{\mu j}$
for $\mu=i$ of the NQS-B $\alpha=1$ ansatz. The dominant diagonal weight
$w^{[1]}_{ii}$ strongly correlating the $\mu=i$ hidden unit with visible site
$n_{i}$ occupation has been removed from the plot.

###  6.2 Case 2: Jastrow + Many-body (J-MB) correlator reference state

Given our aim is to improve on the classic variational ansatz it makes sense
to apply the NQS correlators directly to a reference state that is a pre-
optimised Jastrow + many-body state for each $U/t$. To our knowledge this NQS
combination has not been explored before in the literature. Modifying the
properties of a structured reference state is a very different optimisation
task to the first scenario.

Figure 8: For $U/t=1$: (a) The variational energy density $\epsilon$ obtained
using all five variants NQS-A, -B, -C, -HD and -OH ansatzes applied to a
Jastrow + many-body reference state. For all but the OH ansatz, data is shown
for $\alpha=1$ and $2$ connected by a solid line to guide the eye. The key
comparator is the energy density for Jastrow + many-body correlator (dashed-
dotted line) ansatz and the shaded region above it to the Jastrow energy
density (not shown) is included to guide the eye. (b) The probability $P_{n}$
of occupation for a local configuration state $\left|\,n\right\rangle$ along
with the deviation with J-MB ($\times 10$). (c) The spectrum $p_{\gamma}$ of
the single-particle correlation matrix. (d) The corresponding single-particle
$\langle\hat{b}^{\dagger}_{x,y}\hat{b}_{c}\rangle$, (e) the density-density
$\langle\hat{n}_{x,y}\hat{n}_{c}\rangle$ and (f) the doublon-holon
$\langle\hat{D}_{x,y}\hat{H}_{c}\rangle$ correlations for $x,y$ coordinates of
a site in the $10\times 10$ lattice and $c=(6,6)$ fixed. The results reported
are for NQS-OH ansatz which is the best performing ansatz at $U/t=1$.

In Fig. 8(a) we display the variational energy for all five variants for
$U/t=1$. We find that all NQS ansatz consistently improve on the already high
quality reference state $\left|\,\Psi_{\rm J-MB}\right\rangle$, but the
improvement is small at 0.02%, similar to the best achieved with BEC reference
in Fig. 5(a) and Fig. 6(a). This reflects how $\left|\,\Psi_{\rm
J-MB}\right\rangle$ is only marginally altered from
$\left|\,\Phi_{0}\right\rangle$ in this regime. Nonetheless the J-MB reference
state has now enabled NQS-HD to perform well in this regime.

To examine the ground states beyond energy further insightful physical
observables are plotted in Fig. 8 for the NQS-OH ansatz, which is the best
performing for $U/t=1$. In Fig. 8(b) the probability of local occupation
states are shown. For $U/t=1$ a close approximation to a Poisson distribution
is seen, as expected for a SF. In Fig. 8(c) we plot eigenspectrum $p_{\gamma}$
of the single-particle $\langle\hat{b}^{\dagger}_{i}\hat{b}_{j}\rangle$
correlations displayed in Fig. 8(d). For $U/t=1$ we see the expected
characteristics of a SF, with long-range off-diagonal order in
$\langle\hat{b}^{\dagger}_{i}\hat{b}_{j}\rangle$ for the SF and a 96%
condensate fraction. We show in Fig. 8(e) the density-density
$\langle\hat{n}_{i}\hat{n}_{j}\rangle$ and Fig. 8(f) the doublon-holon
$\langle\hat{D}_{i}\hat{H}_{j}\rangle$ correlations. These indicate
significant density fluctuations and unbound doublon-holon pairs. The
differences in these NQS correlations compared to J-MB alone are rather small
and not visible on the scale of these plots.

Figure 9: For $U/t=32$: (a) The variational energy density $\epsilon$ obtained
using all five variants NQS-A, -B, -C, -HD and -OH ansatzes applied to a
Jastrow + many-body reference state. For all but the OH ansatz, data is shown
for $\alpha=1$ and $2$ connected by a solid line to guide the eye. The key
comparator is the energy density for Jastrow + many-body correlator (dashed-
dotted line) ansatz and the shaded region above it to the Jastrow energy
density (not shown) is included to guide the eye. (b) The probability $P_{n}$
of occupation for a local configuration state $\left|\,n\right\rangle$ along
with the deviation with J-MB ($\times 10$). (c) The spectrum $p_{\gamma}$ of
the single-particle correlation matrix. (d) The corresponding single-particle
$\langle\hat{b}^{\dagger}_{x,y}\hat{b}_{c}\rangle$, (e) the density-density
$\langle\hat{n}_{x,y}\hat{n}_{c}\rangle$ and (f) the doublon-holon
$\langle\hat{D}_{x,y}\hat{H}_{c}\rangle$ correlations for $x,y$ coordinates of
a site in the $10\times 10$ lattice and $c=(6,6)$ fixed. The results reported
are for NQS-HD $\alpha=2$ ansatz which is the best performing ansatz at
$U/t=32$.

The same quantities are plotted in Fig. 9 for $U/t=32$. All NQS ansatzes now
consistently improve on the J-MB reference state, as shown in Fig. 9(a). Deep
in the MI regime a more significant energy improvement of 3% is obtained. The
best performing ansatz is now the NQS-HD at $\alpha=2$, although this is only
marginally cheaper and better than NQS-OH with $\alpha=1$. The local number
distribution reported in Fig. 9(b) is a close approximation to a unit-filled
atomic distribution. The single-particle correlation spectrum in Fig. 9(c)
shows a flat distribution of $p_{\gamma}$ giving a small condensate fraction
$n_{c}=2.3$. As expected for a quench deep into the MI regime compared to Fig.
8(d) we see that $\langle\hat{b}^{\dagger}_{i}\hat{b}_{j}\rangle$ becomes
short-ranged in Fig. 9(d), local density fluctuations
$\langle\hat{n}^{2}_{i}\rangle$ in $\langle\hat{n}_{i}\hat{n}_{j}\rangle$ are
suppressed in Fig. 9(e), and the binding of holon-doublon to nearest-neighbour
pairs in $\langle\hat{D}_{i}\hat{H}_{j}\rangle$ in Fig. 9(f).

The effect of NQS optimisation on top of J-MB has caused the redistribution of
probability from the state $\left|\,1\right\rangle$ to
$\left|\,0\right\rangle$ and $\left|\,2\right\rangle$ in Fig. 9(b) causing an
increase in the local density fluctuations, and a slight shift of population
$p_{\gamma}$ into the dominant modes in Fig. 9(c). The NQS is thus able to
describe a MI with stronger short- and long-ranged density fluctuations than
J-MB alone permits.

Figure 10: For $U/t=16$: (a) The variational energy density $\epsilon$
obtained using all five variants NQS-A, -B, -C, -HD and -OH ansatzes applied
to a Jastrow + many-body reference state. For all but the OH ansatz, data is
shown for $\alpha=1$ and $2$ connected by a solid line to guide the eye. The
key comparator is the energy density for Jastrow + many-body correlator
(dashed-dotted line) ansatz and the shaded region above it to the Jastrow
energy density (not shown) is included to guide the eye. (b) The probability
$P_{n}$ of occupation for a local configuration state $\left|\,n\right\rangle$
along with the deviation with J-MB ($\times 10$). (c) The spectrum
$p_{\gamma}$ of the single-particle correlation matrix. (d) The corresponding
single-particle $\langle\hat{b}^{\dagger}_{x,y}\hat{b}_{c}\rangle$, (e) the
density-density $\langle\hat{n}_{x,y}\hat{n}_{c}\rangle$ and (f) the doublon-
holon $\langle\hat{D}_{x,y}\hat{H}_{c}\rangle$ correlations for $x,y$
coordinates of a site in the $10\times 10$ lattice and $c=(6,6)$ fixed. The
results reported are for NQS-OH ansatz which is the best performing ansatz at
$U/t=16$.

Moving to the boundaries of the critical region we display in Fig. 10 the
results for $U/t=16$. The variational improvement over J-MB shown in Fig.
10(a) is the greatest of all calculations at 4%. This was achieved by the NQS-
OH ansatz, although the performance of all variants was similarly good.
Despite the improvement at $U/t=16$, Fig. 10(b)-(c) show only a mild changes
in the local occupation distribution and the single-particle correlation
spectrum from J-MB. The condensate fraction of $28\%$ and the long-ranged
correlations visible in Fig. 10(d)-(f) are all consistent with $U/t=16$ being
in the SF regime, just as J-MB predicts.

Figure 11: For $U/t=23$: (a) The variational energy density $\epsilon$
obtained using all five variants NQS-A, -B, -C, -HD and -OH ansatzes applied
to a Jastrow + many-body reference state. For all but the OH ansatz, data is
shown for $\alpha=1$ and $2$ connected by a solid line to guide the eye. The
key comparator is the energy density for Jastrow + many-body correlator
(dashed-dotted line) ansatz and the shaded region above it to the Jastrow
energy density (not shown) is included to guide the eye. (b) The probability
$P_{n}$ of occupation for a local configuration state $\left|\,n\right\rangle$
along with the deviation with J-MB ($\times 10$). (c) The spectrum
$p_{\gamma}$ of the single-particle correlation matrix. (d) The corresponding
single-particle $\langle\hat{b}^{\dagger}_{x,y}\hat{b}_{c}\rangle$, (e) the
density-density $\langle\hat{n}_{x,y}\hat{n}_{c}\rangle$ and (f) the doublon-
holon $\langle\hat{D}_{x,y}\hat{H}_{c}\rangle$ correlations for $x,y$
coordinates of a site in the $10\times 10$ lattice and $c=(6,6)$ fixed. The
results reported are for NQS-HD $\alpha=2$ ansatz which is the best performing
ansatz at $U/t=23$.

Finally, the plots for $U/t=23$ are shown in Fig. 11. Here the variational
improvement over J-MB reported in Fig. 11(a) reaches 3.2%, with the NQS-HD
ansatz the best performer. This NQS ansatz makes only very small changes to
both the local occupation distribution in Fig. 11(b) and the single-particle
correlation spectrum in Fig. 11(c) of J-MB. The small condensate fraction of
3.7% and the short-ranged correlations visible in Fig. 10(d)-(f) are all
consistent with $U/t=23$ being in the MI regime, again as J-MB predicts.

Across both the SF and MI regimes the improvement in the variational energy
$\epsilon$ originates from NQS increasing the neighbouring single-particle
correlations and increasing the local density fluctuations above those of the
reference J-MB state. These modifications to the state reduce $\epsilon$ by
describing an elevated negative kinetic energy that outweighs the
commensurately increased positive interaction energy. These changes are
quantitatively significant when interactions are strong, showing that the NQS
and J-MB combination can capture a more complex competition between kinetic
and interaction energies.

##  7 Conclusion

Jastrow and many-body correlators are hugely successful at capturing some of
the keys physics of the BHM. Yet despite enormous improvements over a mean-
field description there are still discrepancies between their predictions for
the critical point compared to other approaches. The SF-MI critical point
$U_{c}/t\approx 20.5$ predicted from the simple many-body ansatz
$\left|\,\Psi_{\rm MB}\right\rangle$ in Eq. (15) is essentially unchanged by
an analysis using the more powerful and encompassing Jastrow + many-body
ansatz $\left|\,\Psi_{\rm J-MB}\right\rangle$ in Eq. (18). Here we have
explored using NQS as a means of systematically improving these variational
wavefunctions. We have found substantial benefits can be gained by introducing
some physical specificity to the generic one-hot encoding RBM ansatz for
bosons. This not only reduces the number of variational parameters to be
optimised but also makes their function more transparent. We introduced a
number of truncated variants guided by their ability to exactly capture
Gutzwiller, Jastrow and many-body correlators, whilst also providing a
systematic means of expanding beyond those cases. Despite analytically exact
reproductions of the classic ansatzes being possible with the NQS variants, we
find that numerical optimisation of NQS do not easily converge to them. The
one exception being NQS-B ansatz which can “learn” the Jastrow state.

We found that NQS correlators, while in principle powerful enough to describe
BHM ground states from a simple BEC reference state, performed much better
when inducing smaller refinements on top of a pre-optimised J-MB reference
state. Specifically, both the consistency of the numerical optimisation and
the final variational energies obtained were improved in this scenario. As
might be expected, the best performance was extracted from the most
expressive, and also most expensive, ansatzes such as NQS-HD and NQS-OH.
Outside the weakly-interacting regime NQS-HD with J-MB reference is the stand-
out choice. Yet, the balance of performance verses optimisation complexity
across $U/t$ also favours the use of the much simpler NQS-B with a J-MB
reference state. Moreover, the simplicity of this ansatz allows an increase in
complexity to be explored with much larger $\alpha$ than exploited here.

Conclusive evidence for how NQS combined with a J-MB reference state changes
the prediction of the BHM critical point would require a more exhaustive scan
across $U/t$ along with finite size scaling, beyond the scope of this current
work. However, data from the judicious points explored here indicate that NQS
for the parameters used, while subtly modifying the J-MB reference state’s
local properties, do not appear to change the regime the state would be
identified as compared to J-MB alone. This suggests the critical point will
remain $U_{c}/t\approx 20.5$ as determined by the reference state. It would be
interesting if future work exploiting a much larger $\alpha$ can refute this
observation.

## Acknowledgements

S.R.C. gratefully acknowledge financial support from UK’s Engineering and
Physical Sciences Research Council (EPSRC) under the grants EP/P025110/2 and
EP/T028424/1. M.P. also acknowledges the University of Bristol’s Advanced
Computing Research Centre (ACRC) for the use of their High Performance
Computing facility (BluePebble) to perform the VMC calculations presented.

##  Appendix A Jastrow and many-body correlator in NQS-HD

The holon-doublon variant RBM correlator is

$\displaystyle\hat{C}_{\rm HD}$ $\displaystyle=$
$\displaystyle\sum_{{\bm{h}}\in\mathbb{N}_{1}^{M}}\exp\left[\sum^{L}_{i=1}\left(a^{[0]}_{i}\hat{H}_{i}+a^{[2]}_{i}\hat{D}_{i}+a^{[m]}_{i}\hat{M}_{i}\right)+\sum_{\mu=1}^{M}b_{\mu}h_{\mu}\right.$
(49)
$\displaystyle\qquad\qquad\qquad+\left.\sum_{\mu=1}^{M}\sum_{i=1}^{L}\left(w^{[0]}_{\mu
i}h_{\mu}\hat{H}_{i}+w^{[1]}_{\mu j}h_{\mu}D_{i}\right)\right].$

For configurations in $\mathcal{F}_{N,2}$ this variant can exactly reproduce
the Jastrow density-density correlator following the construction introduced
in Eq. (35). Specifically, in this subspace we have

$\nu_{ij}\hat{n}_{i}\hat{n}_{j}=\nu_{ij}\left(1-\hat{H}_{i}-\hat{H}_{j}+\hat{D}_{i}+\hat{D}_{j}-\hat{H}_{i}\hat{D}_{j}+\hat{H}_{i}\hat{H}_{j}-\hat{D}_{i}\hat{H}_{j}+\hat{D}_{i}\hat{D}_{j}\right).$
(50)

For a fixed $i$ we then set the biases as

$a^{[0]}_{j}=\left\\{\begin{array}[]{ll}-\sum_{k\neq
i}\nu_{ik}-\mathcal{S}&j=i\\\ -\nu_{ij}&j\neq i\end{array}\right.,\quad
a^{[2]}_{j}=\left\\{\begin{array}[]{ll}\sum_{k\neq
i}\nu_{ik}-\mathcal{S}&j=i\\\ \nu_{ij}&j\neq i\end{array}\right.,$ (51)

which accounts for the onsite terms in Eq. (50). Next, we introduce two hidden
units, one labelled $\mu_{H}$, correlated to the holon occupation on site $i$,
and another labelled $\mu_{D}$, correlated to the doublon occupation on site
$i$. The biases and weights then follow as

$\displaystyle b_{\mu_{H}}=-\mathcal{S},\quad w^{[0]}_{\mu_{H}j}$
$\displaystyle=$ $\displaystyle\left\\{\begin{array}[]{ll}2\mathcal{S}&j=i\\\
\nu_{ij}&j\neq i\end{array}\right.,\quad
w^{[2]}_{\mu_{H}j}=\left\\{\begin{array}[]{ll}0&j=i\\\ -\nu_{ij}&j\neq
i\end{array}\right.,$ $\displaystyle b_{\mu_{D}}=-\mathcal{S},\quad
w^{[0]}_{\mu_{D}j}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&j=i\\\ -\nu_{ij}&j\neq
i\end{array}\right.,\quad
w^{[2]}_{\mu_{D}j}=\left\\{\begin{array}[]{ll}2\mathcal{S}&j=i\\\
\nu_{ij}&j\neq i\end{array}\right..$

After taking the limit $\mathcal{S}\rightarrow\infty$ hidden unit $\mu_{H}$
generates the $\hat{H}_{i}\hat{H}_{j}$, $\hat{H}_{i}\hat{D}_{j}$ terms, while
$\mu_{D}$ generates the $\hat{D}_{i}\hat{H}_{j}$, $\hat{D}_{i}\hat{D}_{j}$
terms. Altogether the correlator $\exp(\sum_{j\neq
i}\nu_{ij}\hat{n}_{i}\hat{n}_{j})$ is reproduced inside the subspace
$\mathcal{F}_{N,2}$.

The many-body correlator can be described exactly with this ansatz by
modifying the construction introduced in Eq. (37). Specifically, for a fixed
site $i$ a correlator $\exp(-\xi\hat{Q}_{i})$ is represented by a pair of
hidden units with the following bias and weights

$\displaystyle
b_{\mu_{H}}=-\mathcal{S}+\log(e^{-\xi}-1),~{}w^{[0]}_{\mu_{H}j}=\left\\{\begin{array}[]{cl}\mathcal{S}&j=i\\\
0&{\rm
otherwise}\end{array}\right.,~{}w^{[2]}_{\mu_{H}j}=\left\\{\begin{array}[]{cl}-\mathcal{S}&j\in\langle
i,j\rangle\\\ 0&{\rm otherwise}\\\ \end{array}\right.,$ $\displaystyle
b_{\mu_{D}}=-\mathcal{S}+\log(e^{-\xi}-1),~{}w^{[0]}_{\mu_{D}j}=\left\\{\begin{array}[]{cl}-\mathcal{S}&j\in\langle
i,j\rangle\\\ 0&{\rm
otherwise}\end{array}\right.,~{}w^{[2]}_{\mu_{D}j}=\left\\{\begin{array}[]{cl}\mathcal{S}&j=i\\\
0&{\rm otherwise}\\\ \end{array}\right..$

The hidden unit $\mu_{H}$ generates the correlator
$\exp[-\xi\hat{H}_{i}\prod_{j\in\langle i,j\rangle}(1-\hat{D}_{j})]$ while
$\mu_{D}$ generates $\exp[-\xi\hat{D}_{i}\prod_{j\in\langle
i,j\rangle}(1-\hat{H}_{j})]$.

Figure 12: A depiction of the weights and biases required for the two hidden
units realising (a) a Jastrow correlator and (b) a many-body correlator. In
(a) $\chi=\nu_{ij}+\nu_{ik}$. In both cases $\mathcal{S}\gg 1$.

##  Appendix B Jastrow correlator in NQS-B

Defining the contribution of the $\mu$th hidden unit to the NQS-B ansatz as

$\displaystyle\hat{\Upsilon}_{\mu}(b^{[1]}_{\mu},b^{[2]}_{\mu},\tilde{{\bm{w}}}^{[1]}_{\mu})=\sum_{h_{\mu}\in\mathbb{N}_{B}}\exp\left[b^{[1]}_{\mu}h_{\mu}+b^{[2]}_{\mu}h_{\mu}^{2}+\sum_{j=1}^{L}\tilde{w}^{[1]}_{\mu
j}h_{\mu}\hat{n}_{j}\right],$

then for a fixed site $i$ the construction in Eq. (35) generalises to biases
and weights

$\tilde{b}^{[1]}_{\mu}=0,\quad\tilde{b}^{[2]}_{\mu}=-\mathcal{S},\quad\tilde{w}^{[1]}_{\mu
j}=\left\\{\begin{array}[]{cc}2\mathcal{S}&j=i\\\ \nu_{ij}&j\neq
i\end{array}\right.,$ (56)

which reduce this contribution to

$\displaystyle\qquad\hat{\Upsilon}_{\mu}(b^{[1]}_{\mu},b^{[2]}_{\mu},\tilde{{\bm{w}}}^{[1]}_{\mu})\left|\,\bm{n}\right\rangle=\exp[\mathcal{S}n_{i}^{2}]\sum_{h_{\mu}\in\mathbb{N}_{B}}\exp\left[-\mathcal{S}(n_{i}-h_{\mu})^{2}\right]\exp\left[\sum_{j\neq
i}^{L}\nu_{ij}h_{\mu}n_{j}\right]\left|\,\bm{n}\right\rangle.$

Since
$\lim_{\mathcal{S}\rightarrow\infty}\exp\left[-\mathcal{S}(n_{i}-h_{\mu})^{2}\right]=\delta_{h_{\mu}n_{i}}$
then, after shifting the quadratic bias of site $i$ as
$\tilde{a}^{[2]}_{i}\mapsto\tilde{a}^{[2]}_{i}-\mathcal{S}$, we find that

$\displaystyle\lim_{\mathcal{S}\rightarrow\infty}\exp[(\tilde{a}^{[2]}_{i}-\mathcal{S})\hat{n}_{i}^{2}]\hat{\Upsilon}_{\mu}(b_{\mu},\tilde{{\bm{w}}}^{[1]}_{\mu})=\exp[\tilde{a}^{[2]}_{i}\hat{n}_{i}^{2}]\exp\left[\sum_{j\neq
i}^{L}\nu_{ij}\hat{n}_{i}\hat{n}_{j}\right],$

and thus exactly reproduce the Jastrow correlations of site $i$ to all other
sites with one hidden unit. The simplicity of this extension to the
construction in Eq. (35) strongly justifies expanding the bandwidth of the
hidden unit and introducing its own quadratic bias.

##  Appendix C Many-body correlator in NQS-C

The variant RBM with a quadratic interaction is

$\displaystyle\hat{C}_{\rm C}$
$\displaystyle=\sum_{{\bm{h}}\in\mathbb{N}_{B}^{M}}\exp\left[\sum_{i=1}^{N}\left(\tilde{a}^{[1]}_{i}\hat{n}_{i}+\tilde{a}^{[2]}_{i}\hat{n}^{2}_{i}\right)+\sum_{\mu=1}^{M}\left(b^{[1]}_{\mu}h_{\mu}+b^{[2]}_{\mu}h_{\mu}^{2}\right)\right.$
$\displaystyle\qquad\qquad\qquad\quad+\left.\sum_{\mu=1}^{M}\sum_{i=1}^{L}\left(\tilde{w}^{[1]}_{\mu
i}h_{\mu}\hat{n}_{i}+\tilde{w}^{[2]}_{\mu
i}h_{\mu}\hat{n}^{2}_{i}\right)\right].$ (57)

Using a pair of hidden units a modified version of the construction in Eq.
(37) this variant can describe the many-body correlator in the subspace
$\mathcal{F}_{N,2}$. This is done with the following biases and weights

$\displaystyle\quad\begin{array}[]{l}b^{[1]}_{\mu_{H}}=\log(\tau)\\\
b^{[2]}_{\mu_{H}}=0\end{array},\qquad\quad~{}\tilde{w}^{[1]}_{\mu_{H}j}=\left\\{\begin{array}[]{cl}0&j=i\\\
\mathcal{S}&j\in\langle i,j\rangle\\\ 0&{\rm
otherwise}\end{array}\right.,~{}\tilde{w}^{[2]}_{\mu_{H}j}=\left\\{\begin{array}[]{cl}-\mathcal{S}&j=i\\\
-\mathcal{S}&j\in\langle i,j\rangle\\\ 0&{\rm otherwise}\end{array}\right.,$
$\displaystyle\quad\begin{array}[]{l}b^{[1]}_{\mu_{D}}=-12\mathcal{S}+\log(\tau)\\\
b^{[2]}_{\mu_{D}}=0\end{array},~{}\tilde{w}^{[1]}_{\mu_{D}j}=\left\\{\begin{array}[]{cl}4\mathcal{S}&j=i\\\
3\mathcal{S}&j\in\langle i,j\rangle\\\ 0&{\rm
otherwise}\end{array}\right.,~{}\tilde{w}^{[2]}_{\mu_{D}j}=\left\\{\begin{array}[]{cl}-\mathcal{S}&j=i\\\
-\mathcal{S}&j\in\langle i,j\rangle\\\ 0&{\rm otherwise}\\\
\end{array}\right..$

where $\tau$ is a root of the equation $\sum_{p=1}^{B}x^{p}-e^{-\xi}+1=0$. To
unravel this construction consider the contribution of hidden unit $\mu_{H}$

$\displaystyle\hat{\Upsilon}_{\mu_{H}}(b^{[1]}_{\mu_{H}},b^{[2]}_{\mu_{H}},\tilde{\bm{w}}^{[1]}_{\mu_{H}},\tilde{\bm{w}}^{[2]}_{\mu_{H}})\left|\,\bm{n}\right\rangle=$
$\displaystyle\sum_{h_{\mu_{H}}\in\mathbb{N}_{B}}\tau^{h_{\mu_{H}}}\exp\left[-\mathcal{S}\left(n_{i}^{2}+\sum_{j\in\langle
i,j\rangle}n_{j}(n_{j}-1)\right)h_{\mu_{H}}\right]\left|\,\bm{n}\right\rangle.$

For the terms with $h_{\mu_{H}}>0$ the $\exp[\cdots]$ factor is unity only if
there is a holon on site $i$ and neigbouring sites are only combinations of
holons or singlons. Otherwise, in the limit $\mathcal{S}\rightarrow\infty$,
this factor is zero. Altogether the hidden unit reproduces the correlator
$\exp[-\xi\hat{H}_{i}\prod_{j\in\langle i,j\rangle}(1-\hat{D}_{j})]$ within
$\mathcal{F}_{N,2}$. Similar logic applies for hidden unit $\mu_{D}$ whose
contribution is

$\displaystyle\hat{\Upsilon}_{\mu_{D}}(b^{[1]}_{\mu_{D}},b^{[2]}_{\mu_{D}},\tilde{\bm{w}}^{[1]}_{\mu_{D}},\tilde{\bm{w}}^{[2]}_{\mu_{D}})\left|\,\bm{n}\right\rangle=$
$\displaystyle\sum_{h_{\mu_{D}}\in\mathbb{N}_{B}}\tau^{h_{\mu_{D}}}\exp\left[-\mathcal{S}\left((n_{i}-2)^{2}+\sum_{j\in\langle
i,j\rangle}(n_{j}-1)(n_{j}-2)\right)h_{\mu_{D}}\right]\left|\,\bm{n}\right\rangle,$

which reproduces the correlator $\exp[-\xi\hat{D}_{i}\prod_{j\in\langle
i,j\rangle}(1-\hat{H}_{j})]$ within $\mathcal{F}_{N,2}$.

## References

## References

  * [1] Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885 URL https://link.aps.org/doi/10.1103/RevModPhys.80.885
  * [2] Greiner M, Mandel O, Esslinger T, Hänsch T W and Bloch I 2002 Nature 415 39 URL https://doi.org/10.1038/415039a
  * [3] Jaksch D, Bruder C, Cirac J I, Gardiner C W and Zoller P 1998 Phys. Rev. Lett. 81 3108 URL http://dx.doi.org/10.1103/PhysRevLett.81.3108
  * [4] Jördens R, Strohmaier N, Günter K, Moritz H and Esslinger T 2008 Nature 455 204–207 URL https://doi.org/10.1038/nature07244
  * [5] Gemelke N, Zhang X, Hung C L and Chin C 2009 Nature 460 995 URL https://doi.org/10.1038/nature08244
  * [6] Spielman I B, Phillips W D and Porto J V 2007 Phys. Rev. Lett. 98 080404 URL https://link.aps.org/doi/10.1103/PhysRevLett.98.080404
  * [7] Bakr W S, Peng A, Tai M E amd Ma R, Simon J, Gillen J I, Fölling S, Pollet L and Greiner M 2010 Science 329 547 URL https://www.science.org/doi/abs/10.1126/science.1192368
  * [8] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 Phys. Rev. B 40 546 URL https://link.aps.org/doi/10.1103/PhysRevB.40.546
  * [9] Sachdev S 2011 Quantum Phase Transitions, 2nd ed. (Cambridge University Press) URL https://doi.org/10.1017/CBO9780511973765
  * [10] Schreiber M, Hodgman S S, Bordia P, Luschen H P, Fischer M H, Vosk R, Altman E, Schneider U and Bloch I 2015 Science 349 842 URL http://dx.doi.org/10.1126/science.aaa7432
  * [11] Lüschen H P, Bordia P, Scherg S, Alet F, Altman E, Schneider U and Bloch I 2017 Phys. Rev. Lett. 119 260401 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.260401
  * [12] Choi J y, Hild S, Zeiher J, Schauss P, Rubio-Abadal A, Yefsah T, Khemani V, Huse D A, Bloch I and Gross C 2016 Science 352 1547 URL http://dx.doi.org/10.1126/science.aaf8834
  * [13] Wahl T B, Pal A and Simon S H 2019 Nature Physics 15 164 URL https://doi.org/10.1038/s41567-018-0339-x
  * [14] Lukin A, Rispoli M, Schittko R, Tai M E, Kaufman A M, Choi S, Khemani V, Léonard J and Greiner M 2019 Science 364 256 URL https://www.science.org/doi/abs/10.1126/science.aau0818
  * [15] Rispoli M, Lukin A, Schittko R, Kim S, Tai M E, Léonard J and Greiner M 2019 Nature 573 385 URL https://doi.org/10.1038/s41586-019-1527-2
  * [16] Struck J, Ölschläger C, Targat R L, Soltan-Panahi P, Eckardt A, Lewenstein M, Windpassinger P and Sengstock K 2011 Science 333 996–999 (Preprint https://www.science.org/doi/pdf/10.1126/science.1207239) URL https://www.science.org/doi/abs/10.1126/science.1207239
  * [17] Saugmann P, Vargas J, Kiefer Y, Hachman M, Eichberger R, Hemmerich A and Larson J 2022 Phys. Rev. A 106(4) L041302 URL https://link.aps.org/doi/10.1103/PhysRevA.106.L041302
  * [18] Tai M E, Lukin A, Rispoli M, Schittko R, Menke T, Borgnia D, Preiss P M, Grusdt F, Kaufman A M and Greiner M 2017 Nature 546 519 URL http://dx.doi.org/10.1038/nature22811
  * [19] Łącki M, Pichler H, Sterdyniak A, Lyras A, Lembessis V E, Al-Dossary O, Budich J C and Zoller P 2016 Phys. Rev. A 93 013604 URL https://link.aps.org/doi/10.1103/PhysRevA.93.013604
  * [20] Palmer R N and Jaksch D 2006 Phys. Rev. Lett. 96 180407 URL https://link.aps.org/doi/10.1103/PhysRevLett.96.180407
  * [21] Hafezi M, Sørensen A S, Demler E and Lukin M D 2007 Phys. Rev. A 76 023613 URL https://link.aps.org/doi/10.1103/PhysRevA.76.023613
  * [22] Becca F and Sorella S 2017 Quantum Monte Carlo Approaches for Correlated Systems (Cambridge University Press) URL https://doi.org/10.1017/9781316417041
  * [23] Gubernatis J, Kawashima N and Werner P 2016 Quantum Monte Carlo Methods: Algorithms for Lattice Models (Cambridge University Press) URL https://doi.org/10.1017/CBO9780511902581
  * [24] Zwerger W 2003 Journal of Optics B: Quantum and Semiclassical Optics 5 S9 URL https://dx.doi.org/10.1088/1464-4266/5/2/352
  * [25] Yokoyama H and Ogata M 2008 Journal of Physics and Chemistry of Solids 69 3356 URL http://dx.doi.org/10.1016/j.jpcs.2008.06.087
  * [26] Yokoyama H, Miyagawa T and Ogata M 2011 Journal of the Physical Society of Japan 80 084607 URL http://dx.doi.org/10.1143/JPSJ.80.084607
  * [27] McMillan W L 1965 Phys. Rev. 138 A442 URL https://link.aps.org/doi/10.1103/PhysRev.138.A442
  * [28] Ceperley D, Chester G V and Kalos M H 1977 Phys. Rev. B 16 3081 URL https://link.aps.org/doi/10.1103/PhysRevB.16.3081
  * [29] Carleo G and Troyer M 2017 Science 355 602606 URL https://science.sciencemag.org/content/355/6325/602
  * [30] Gao X and Duan L M 2017 Nature Communications 8 662 URL https://www.nature.com/articles/s41467-017-00705-2
  * [31] Carleo G, Nomura Y and Imada M 2018 Nature Communications 5322 URL https://www.nature.com/articles/s41467-018-07520-3
  * [32] He H, Zheng Y, Bernevig B A and Sierra G 2019 ArXiv preprint, [cond-mat.str-el] 1910.13454 URL https://arxiv.org/abs/1910.13454
  * [33] Choo K, Neupert T and Carleo G 2019 Phys. Rev. B 100 125124 URL https://link.aps.org/doi/10.1103/PhysRevB.100.125124
  * [34] Irikura N and Saito H 2020 Phys. Rev. Research 2 013284 URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013284
  * [35] Schmitt M and Heyl M 2020 Phys. Rev. Lett. 125 100503 URL https://link.aps.org/doi/10.1103/PhysRevLett.125.100503
  * [36] Liang X, Dong S J and He L 2021 Phys. Rev. B 103 035138 URL https://link.aps.org/doi/10.1103/PhysRevB.103.035138
  * [37] Saito H and Kato M 2018 Journal of the Physical Society of Japan 87 014001 (Preprint https://doi.org/10.7566/JPSJ.87.014001) URL https://doi.org/10.7566/JPSJ.87.014001
  * [38] Choo K, Carleo G, Regnault N and Neupert T 2018 Phys. Rev. Lett. 121 167204 URL https://link.aps.org/doi/10.1103/PhysRevLett.121.167204
  * [39] Luo D and Clark B K 2019 Phys. Rev. Lett. 122 226401 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.226401
  * [40] Adams C, Carleo G, Lovato A and Rocco N 2021 Phys. Rev. Lett. 127 022502 ISSN 1079-7114 URL http://dx.doi.org/10.1103/PhysRevLett.127.022502
  * [41] Levine Y, Sharir O, Cohen N and Shashua A 2019 Phys. Rev. Lett. 122 065301 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.065301
  * [42] Tubman N M 2016 ArXiv preprint, [cond-mat.str-el] 1609.08142
  * [43] Torlai G, Mazzola G, Carrasquilla J, Troyer M, Melko R and Carleo G 2018 Nature Physics 14 447 URL https://doi.org/10.1038/s41567-018-0048-5
  * [44] Torlai G, Timar B, van Nieuwenburg E P L, Levine H, Omran A, Keesling A, Bernien H, Greiner M, Vuletić V, Lukin M D, Melko R G and Endres M 2019 Phys. Rev. Lett. 123 230504 URL https://link.aps.org/doi/10.1103/PhysRevLett.123.230504
  * [45] Carrasquilla J, Torlai G, Melko R G and Aolita L 2019 Nature Machine Intelligence 1 155 URL http://dx.doi.org/10.1038/s42256-019-0028-1
  * [46] Neugebauer M, Fischer L, Jäger A, Czischek S, Jochim S, Weidemüller M and Gärttner M 2020 Phys. Rev. A 102 042604 URL https://link.aps.org/doi/10.1103/PhysRevA.102.042604
  * [47] Gao J, Qiao L F, Jiao Z Q, Ma Y C, Hu C Q, Ren R J, Yang A L, Tang H, Yung M H and Jin X M 2018 Phys. Rev. Lett. 120 240501 URL https://link.aps.org/doi/10.1103/PhysRevLett.120.240501
  * [48] Dong X Y, Pollmann F and Zhang X F 2019 Phys. Rev. B 99 121104 URL https://link.aps.org/doi/10.1103/PhysRevB.99.121104
  * [49] Carrasquilla J and Melko R G 2017 Nature Physics 13 431 URL https://doi.org/10.1038/nphys4035
  * [50] Schindler F, Regnault N and Neupert T 2017 Phys. Rev. B 95 245134 URL https://link.aps.org/doi/10.1103/PhysRevB.95.245134
  * [51] Harney C, Pirandola S, Ferraro A and Paternostro M 2020 New Journal of Physics 22 045001 URL https://dx.doi.org/10.1088/1367-2630/ab783d
  * [52] Hsu Y T, Li X, Deng D L and Das Sarma S 2018 Phys. Rev. Lett. 121 245701 URL https://link.aps.org/doi/10.1103/PhysRevLett.121.245701
  * [53] Yoshioka N and Hamazaki R 2019 Phys. Rev. B 99 214306 URL https://link.aps.org/doi/10.1103/PhysRevB.99.214306
  * [54] Nagy A and Savona V 2019 Phys. Rev. Lett. 122 URL http://dx.doi.org/10.1103/PhysRevLett.122.250501
  * [55] Hartmann M J and Carleo G 2019 Phys. Rev. Lett. 122 250502 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.250502
  * [56] Le Roux N and Bengio Y 2008 Neural Computation 20 1631 (Preprint https://doi.org/10.1162/neco.2008.04-07-510) URL https://doi.org/10.1162/neco.2008.04-07-510
  * [57] Lu S, Gao X and Duan L M 2019 Phys. Rev. B 99 155136 URL https://link.aps.org/doi/10.1103/PhysRevB.99.155136
  * [58] Glasser I, Pancotti N, August M, Rodriguez I D and Cirac J I 2018 Phys. Rev. X 8 011006 URL https://link.aps.org/doi/10.1103/PhysRevX.8.011006
  * [59] Clark S R 2018 Journal of Physics A: Mathematical and Theoretical 51 135301 URL http://dx.doi.org/10.1088/1751-8121/aaaaf2
  * [60] Pei M Y and Clark S R 2021 Entropy 23 879 URL https://www.mdpi.com/1099-4300/23/7/879
  * [61] Pei M Y and Clark S R 2021 J. Phys. A: Math. Theor. 54 405304 URL https://dx.doi.org/10.1088/1751-8121/ac1f3d
  * [62] McBrian K, Carleo G and Khatami E 2019 Journal of Physics: Conference Series 1290 012005 URL https://doi.org/10.1088/1742-6596/1290/1/012005
  * [63] Vargas-Calderón V, Vinck-Posada H and González F A 2020 Journal of the Physical Society of Japan 89 094002 URL http://dx.doi.org/10.7566/JPSJ.89.094002
  * [64] Nomura Y 2020 Journal of the Physical Society of Japan 89 054706 (Preprint https://doi.org/10.7566/JPSJ.89.054706) URL https://doi.org/10.7566/JPSJ.89.054706
  * [65] Saito H 2017 Journal of the Physical Society of Japan 86 093001 (Preprint https://doi.org/10.7566/JPSJ.86.093001) URL https://doi.org/10.7566/JPSJ.86.093001
  * [66] Guo C and Berkhahn F 2016 ArXiv preprint, [cs.LG] 1604.06737 URL https://arxiv.org/abs/1604.06737
  * [67] Bravyi, Sergey and DiVincenzo, David P and Oliveira, Roberto and Terhal, Barbara M 2008 Quantum Information and Computation 8 0361 URL https://doi.org/10.26421/QIC8.5-1
  * [68] Feynman R 1998 Statistical Mechanics: A Set Of Lectures Advanced Books Classics (Avalon Publishing) ISBN 9780813346106 URL https://books.google.co.uk/books?id=Ou4ltPYiXPgC
  * [69] Lewenstein M, Sanpera A and Ahufinger V 2012 Ultracold Atoms in Optical Lattices: Simulating Quantum Many-Body Systems (Oxford University Press) URL https://doi.org/10.1093/acprof:oso/9780199573127.001.0001
  * [70] Coleman P 2015 Introduction to Many Body Physics (Cambridge University Press) URL https://doi.org/10.1017/CBO9781139020916
  * [71] Stoof H T C, Gubbels K B and Dickerscheid D B M 2009 Ultracold Quantum Fields Theoretical and Mathematical Physics (Springer Netherlands) URL https://doi.org/10.1007/978-1-4020-8763-9
  * [72] Krauth W and Trivedi N 1991 Europhysics Letters (EPL) 14 627 URL https://doi.org/10.1209/0295-5075/14/7/003
  * [73] Elstner N and Monien H 1999 Phys. Rev. B 59 12184 URL https://link.aps.org/doi/10.1103/PhysRevB.59.12184
  * [74] Gros C 1989 Ann. Phys. 189 53 URL http://www.sciencedirect.com/science/article/pii/0003491689900778
  * [75] Foulkes W M C, Mitas L, Needs R J and Rajagopal G 2001 Rev. Mod. Phys. 73 33 URL https://link.aps.org/doi/10.1103/RevModPhys.73.33
  * [76] Krauth W 2006 Statistical Mechanics: Algorithms and Computations Oxford Masters Series in Physics: Statistical, Computational and Theoretical Physics (Oxford University Press) URL https://doi.org/10.1093/oso/9780198515357.001.0001
  * [77] Lou J and Sandvik A W 2007 Phys. Rev. B 76 104432 URL https://link.aps.org/doi/10.1103/PhysRevB.76.104432
  * [78] Nightingale M P and Melik-Alaverdian V 2001 Phys. Rev. Lett. 87 043401 URL https://link.aps.org/doi/10.1103/PhysRevLett.87.043401
  * [79] Toulouse J and Umrigar C J 2007 J. Chem. Phys. 126 084102 (Preprint http://dx.doi.org/10.1063/1.2437215) URL http://dx.doi.org/10.1063/1.2437215
  * [80] Sorella S 2001 Phys. Rev. B 64 024512 URL https://link.aps.org/doi/10.1103/PhysRevB.64.024512
  * [81] Gutzwiller M C 1963 Phys. Rev. Lett. 10 159 URL https://link.aps.org/doi/10.1103/PhysRevLett.10.159
  * [82] van Oosten D, van der Straten P and Stoof H T C 2001 Phys. Rev. A 63 URL http://dx.doi.org/10.1103/PhysRevA.63.053601
  * [83] Freericks J K and Monien H 1994 Europhysics Letters (EPL) 26 545 URL http://dx.doi.org/10.1209/0295-5075/26/7/012
  * [84] Jastrow R 1955 Phys. Rev. 98 1479 URL https://link.aps.org/doi/10.1103/PhysRev.98.1479
  * [85] Capello M, Becca F, Fabrizio M and Sorella S 2007 Phys. Rev. Lett. 99 056402 URL http://dx.doi.org/10.1103/PhysRevLett.99.056402
  * [86] Pitaevskii L and Stringari S 2016 Bose-Einstein Condensation and Superfluidity Oxford Science Publications (Oxford University Press) URL https://doi.org/10.1093/acprof:oso/9780198758884.001.0001
  * [87] Capello M, Becca F, Fabrizio M and Sorella S 2008 Phys. Rev. B 77 144517 URL http://dx.doi.org/10.1103/PhysRevB.77.144517
  * [88] Smolensky P 1986 Information Processing in Dynamical Systems: Foundations of Harmony Theory (Cambridge, MA, USA: MIT Press) p 194 in Parallel Distributed Processing: Explorations in the Microstructure of Cognition, Vol. 1: Foundations.
  * [89] Hinton G E 2002 Neural Comput. 14 1771 URL https://doi.org/10.1162/089976602760128018
  * [90] Deng D L, Li X and Das Sarma S 2017 Phys. Rev. B 96 195145 URL https://link.aps.org/doi/10.1103/PhysRevB.96.195145
  * [91] Kaubruegger R, Pastori L and Budich J C 2018 Phys. Rev. B 97 195136 URL https://link.aps.org/doi/10.1103/PhysRevB.97.195136
  * [92] Zhang Y H, Jia Z A, Wu Y C and Guo G C 2018 ArXiv preprint, [quant-ph] 1809.08631 URL https://arxiv.org/abs/1809.08631
  * [93] Zheng Y, He H, Regnault N and Bernevig B A 2019 Phys. Rev. B 99 155129 URL https://link.aps.org/doi/10.1103/PhysRevB.99.155129
  * [94] Jia Z A, Zhang Y H, Wu Y C, Kong L, Guo G C and Guo G P 2019 Phys. Rev. A 99 012307 URL https://link.aps.org/doi/10.1103/PhysRevA.99.012307
  * [95] Hinton G E and Salakhutdinov R R 2006 Science 31 504 URL https://www.science.org/doi/abs/10.1126/science.1127647
  * [96] Mehta P, Bukov M, Wang C H, Day A G, Richardson C, Fisher C K and Schwab D J 2019 Physics Reports 810 1–124 URL https://www.sciencedirect.com/science/article/pii/S0370157319300766