Phase Diagram Reconstruction of the Bose-Hubbard Model with a Restricted Boltzmann Machine Wavefunction

The BH Hamiltonian describes the interactions between bosons that can occupy sites in a $d$-dimensional lattice. The model is characterized by a hopping energy $t$, an on-site interaction $U$ and a chemical potential $\mu$, so that the grand-canonical Hamiltonian reads ($\hbar=1$) (Fisher et al., 1989)

where $\hat{a}_{i}$ is the annihilation operator at site $i$, and $\hat{n}_{i}=\hat{a}^{\dagger}_{i}\hat{a}_{i}$ is the number operator. The notation $\langle ij\rangle$ indicates that the sum runs over pairs of neighbor sites in the lattice. The BH model is able to reproduce experimental results in Josephson-junction networks (van Oudenaarden and Mooij, 1996; Baltin and Wagenblast, 1997; Bruder et al., 2005) and in lattices of ultra-cold atoms (Greiner et al., 2003; Jaksch et al., 1998; Bloch et al., 2008; Stöferle et al., 2004; Spielman et al., 2007; Gemelke et al., 2009). The latter offers precise control of the lattice parameters (Grynberg et al., 2000; Kollath et al., 2004; Lugan et al., 2007). Theoretically, a lot of attention has been devoted both to understand the quantum phases of the system (the ones that arise from eq. 1 and from the disordered or extended BH model with longer range interactions) (Freericks and Monien, 1996; Ejima et al., 2012; Cazalilla et al., 2011; Teichmann et al., 2009; Batrouni et al., 1995; Krauth et al., 1991; Krauth and Trivedi, 1991; Niyaz et al., 1991). Another central issue is to calculate the quantum phase transition boundaries (Kühner and Monien, 1998; Ejima et al., 2011; Kühner et al., 2000; Batrouni and Scalettar, 1992; Kashurnikov et al., 1996; Batrouni et al., 1990; Elstner and Monien, 1999; Koller and Dupuis, 2006), whose precision has improved over the years with better calculation techniques and computing power, revealing features such as the re-entrance phenomenon (Pino et al., 2013; Elstner and Monien, 1999), where for particular values of the chemical potential, the system switches between the MI phase to the SF phase, and back to the MI phase before definitely entering the SF phase after an increase of $t/U$.

For simplicity, we will restrict our analyses to the $d=1$ case, where only two quantum phases are possible. When the on-site interaction energy is much larger than the hopping energy, the latter becomes negligible, and the Hamiltonian is written as the sum of independent Hamiltonians for each site $\frac{U}{2}\hat{n}_{i}(\hat{n}_{i}-1)-\mu\hat{n}_{i}$. These one-particle Hamiltonians can be immediately diagonalized by the number basis. The corresponding eigenenergies are $\frac{U}{2}n_{i}(n_{i}-1)-\mu n_{i}$, which reach minimum values for fixed $\mu$ and $U$ at $n_{i}=\max\{0,\lceil\mu/U\rceil\}$ (note that all sites are equivalent). Moreover, in this regime, the expected variance of the local number operator is 0, i.e. $\langle\hat{n}_{i}^{2}\rangle-\langle\hat{n}_{i}\rangle^{2}=0$. This regime characterizes the MI phase. On the other hand, when the hopping energy is much larger than the on-site interaction energy, the latter becomes negligible. Thus, the Hamiltonian can also be written as the sum of independent Hamiltonians, but in momentum space, where $\tilde{a}_{k}=N^{-1/2}\sum_{j=1}^{N}\hat{a}_{j}e^{-ix_{j}p_{k}/\hbar}$ is the boson annihilation operator in the momentum representation. Here, $x_{j}=c\times j$, where $c$ is the lattice constant, and $p_{k}=2\pi k\hbar/(N\times c)$. Each independent Hamiltonian in momentum space ($\sum_{k}(-2t\cos(2\pi k/N)-\mu)\tilde{a}^{\dagger}_{k}\tilde{a}_{k}$) has eigenenergies $-2t\cos(2\pi k/N)-\mu$, which reach their minimum when all bosons condense with 0 momenta. Note that the energies are independent of the states’ occupation, meaning that the ground state is degenerate for any number of particles. This regime is known as the SF phase, characterized by a delocalized wavefunction (formally described by algebraic decaying spatial correlations (Ejima et al., 2012; Pino et al., 2013), which is why the SF phase in 1D is not a true Bose-Einstein condensate). Thus, the ground state has a non-zero expected variance of the local number operator. In fact, the probability distribution for the local occupation is Poissonian, meaning that $\expectationvalue*{n_{i}^{2}}-\expectationvalue*{n_{i}}^{2}=\expectationvalue*{n_{i}}$ (Greiner et al., 2003). For a fixed chemical potential and an infinite number of sites, there exists a continuous phase transition from the MI phase to the SF phase as $U$ decreases with respect to $t$ (with the exception of a range of $\mu$ values in the first Mott lobe, where re-entrance exists).

In this section, a review of the method introduced by Carleo and Troyer (2017) is given, where a VMC setup is used to find the ground state, formally written as a trial wavefunction given by an RBM. In the Fock space basis, the ground state of the BH model can be written as Saito (2017)

where $n_{i}$ is the occupation number of the $i$-th site, and $|\Psi(n_{1},n_{2},\ldots)|^{2}$ are the probability amplitudes corresponding to the Fock states $\ket{n_{1},n_{2},\ldots}$. In order to map the wave function to a computer, both the number of sites and the number of possible particles in each site have to be truncated. We will refer to the number of sites as $N$ and to the maximum number of particles in each site as $M-1$. The coefficients $\Psi(n_{1},n_{2},\ldots,n_{N})$ are approximated by an RBM. RBMs are generative neural networks, formally described by a bipartite undirected graph such as the one shown in fig. 1, where there is a layer of visible neurons denoted by ${\bf\it v}$ that are used to input real data, and a layer of hidden neurons denoted by ${\bf\it h}$ that are used as latent variables of the model (Smolensky, 1986). In particular, the wavefunction coefficients take the form $\Psi(n_{1},n_{2},\ldots,n_{N})\approx\psi_{{\bf\it\theta}}(n_{1},\ldots,n_{N})=\sum_{{\bf\it h}}e^{-E_{\text{RBM}}({\bf\it v}({\bf\it n}),{\bf\it h})}$, where $E_{\text{RBM}}({\bf\it v}({\bf\it n}),{\bf\it h})$ is the energy of the RBM (Smolensky, 1986), and ${\bf\it\theta}$ are the variational parameters of the RBM. As a short-hand notation, an occupation configuration is denoted as ${\bf\it n}$, and it is inputted to the visible layer of the RBM. However, the configuration first needs to be one-hot encoded as follows: each occupation $n_{i}$ is encoded into an $M$-component vector whose $j$-th component is $\delta_{j,n_{i}}$, $j=0,1,\ldots,M-1$; then, the vectors for every occupation are concatenated into ${\bf\it v}({\bf\it n})$. Moreover, if the $N_{H}$ hidden neurons are restricted to binary values $1$ or $-1$, then, the approximated wavefunction coefficients can be written as Carleo and Troyer (2017)

where $a_{j},b_{\ell}$ and $W$ are the complex-valued visible bias, hidden bias and connection matrix of the RBM, respectively.

Thus, the approximation of the wavefunction coefficients is done through the adjustment of the parameters ${\bf\it\theta}:\{a_{j},b_{\ell},W_{\ell,j}\}$ that minimize the energy $\bra{\psi_{{\bf\it\theta}}}\hat{H}\ket{\psi_{{\bf\it\theta}}}$. At each step of the minimization, a set $\mathcal{M}$ of configurations ${\bf\it n}$ is sampled from $|\psi_{{\bf\it\theta}}({\bf\it n})|^{2}$ using the Metropolis-Hastings algorithm, so that the energy can be efficiently estimated as (Saito, 2017)

More explicitly, the following steps are carried out to generate the sample $\mathcal{M}$. At the first iteration, a state ${\bf\it n}_{0}$ is randomly proposed. Then, at the $i$-th iteration:

1. 1.

   Under some updating rule, propose a new state ${\bf\it n}^{\prime}_{i}$ from ${\bf\it n}_{i}$.

2. 2.

   With probability $\min\{1,|\psi_{{\bf\it\theta}}({\bf\it n}^{\prime}_{i})/\psi_{{\bf\it\theta}}({\bf\it n}_{i})|^{2}\}$ accept the state ${\bf\it n}^{\prime}_{i}$, i.e. ${\bf\it n}_{i+1}\leftarrow{\bf\it n}^{\prime}_{i}$. If it is not accepted, then ${\bf\it n}_{i+1}\leftarrow{\bf\it n}_{i}$.

To build the sample $\mathcal{M}$ a total of 1000 iterations are performed. The updating rule is provided by NetKet’s Local Metropolis sampler, which picks a site at random and assigns to it a uniformly random occupation, except from the current occupation. Then, through either stochastic gradient descent or stochastic reconfiguration (Sorella et al., 2007), the energy in eq. 4 can be minimized, producing a new set of parameters ${\bf\it\theta}$, as explained by Carleo and Troyer (2017). Both the sampling and minimization of the energy with respect to the RBM parameters are repeated iteratively until the RBM parameters converge, resembling an Expectation Maximization algorithm (Dempster et al., 1977). A schematic representation of the variational Monte Carlo technique is shown in fig. 2.

The blot at the bottom of fig. 2 depicts a region of the Hilbert space that contains the most relevant Fock states $\ket{{\bf\it n}}$ for the ground state of the BH Hamiltonian, with fixed $\mu,t$ and $U$. For instance, the only relevant Fock state within the $m$-th Mott lobe is $\ket{n_{1}=m,n_{2}=m,\ldots}$, whereas in the SF phase there are many relevant Fock states. In that regard, the region of relevant Fock states can contain a variable number of Fock states depending on the parameters of the Hamiltonian that one is intending to solve. Therefore, an issue immediately arises: for an unknown target probability distribution $|\Psi({\bf\it n})|^{2}$, the sampling can be too small or too large with respect to the size of the region of relevant Fock states. If it is small, important information about the ground state might not be taken into account, whereas if it is large, noisy probability from non-relevant states can be taken into account.

Since the VMC minimizes the energy, we check that the energy computed with eq. 4 converges by measuring its variance for the last 500 sampling-optimization steps, as well as by measuring the absolute error of the energy when compared to the exact ground state energy found through Lanczsos diagonalization. For the aforementioned three scenarios, the variance for the last 500 sampling-optimization steps is shown in fig. 3(a)-(c), and the corresponding absolute errors with respect to the exact ground state energy are shown in fig. 3(d)-(f). It is seen that in the Hamiltonian parameter space, the majority of the energies have low-variance, showing convergence towards a value that is in excellent agreement with the exact ground-state energies.

Since the Hilbert spaces are small enough to compute all the probability amplitude coefficients for every state in the Fock basis, we can also compute any observable $\hat{O}$ as

In particular, the absolute error of the energy computed through eq. 5 is shown in fig. 4. It is clear that in the case of 5 sites and 20 neurons, a very large region both in the SF and MI phases shows a strong disagreement between the energy calculated through eqs. 4 and 5, showing that even though the energy converged, the state did not. Another recurrent pattern in the energies calculated through eqs. 4 and 5 is an arc of high absolute errors formed in the left-most side of the Mott lobes, which we will address later.

The lack of convergence of the RBM state for the case of 5 sites and 20 hidden neurons is further confirmed when we measure the overlap between the exact ground state $\ket{\psi_{\text{exact}}}$ and the RBM ground state $\ket{\psi_{{\bf\it\theta}}}$, shown in fig. 5(c). It is now clear why the expected energy with respect to the complete RBM state shown in fig. 4(c) presents large errors when compared to the expected energy with respect to the exact ground state: it appears that the RBM has not learned the ground state in the bottom-right region of the plot, which is a region that covers part of the Mott lobe, as well as part of the SF phase region. Moreover, in the case of 5 sites and 8 hidden neurons, and the case of 8 sites and 11 hidden neurons, the RBM finds difficulty in learning the ground state in the limit between the MI and the SF phase as shown in fig. 5(a) and (b). The difficulty in learning those states, and in general, in treating the ground state near the MI-SF boundary comes from the Kosterlitz-Thouless-like quantum phase transition in 1D systems (Giamarchi, 2003), where an exponentially small Mott gap exists (Ejima et al., 2012). We see once again that there are arcs of low overlap points formed in the left-most side of the Mott lobes. Within the SF phase, there are also lines of low overlap, which appear because of finite size effects. Note that there are as many of these fictitious boundaries as there are sites in the periodic chain under study. Nevertheless, when comparing the 5 sites cases, it is seen that for values of $\mu/U>1$ the states at the MI-SF phases boundary are better learned when 20 hidden neurons are used in the RBM (as in (McBrian et al., 2019), cf. fig. 5(c)) than when only 8 hidden neurons are used (see fig. 5(a)).

An interesting motif is that if one pays very close attention to the energy absolute errors in fig. 3(d)-(f), the lowest errors occur where the overlap in fig. 5(a)-(c) is lowest. Therefore, even though the convergence in the energy is excellent, the state is not learned correctly.

Since the advantage of VMC over exact diagonalization comes for intractably large Hilbert spaces, it is not always possible to compute the probability amplitudes for all of the Fock states basis. In such a case, the RBM state can be built as

where $\mathcal{M}^{\prime}$ is a sample of Fock states sampled with the Metropolis-Hastings algorithm, with an acceptance probability of $\min\{1,p_{\text{GC}}({\bf\it n}_{i+1})/p_{\text{GC}}({\bf\it n}_{i})\}$, where $p_{\text{GC}}({\bf\it n})=\mathcal{Z}^{-1}\exp(-\expectationvalue{\hat{H}}{{\bf\it n}})$ ²²2If instead the sample is built from the distribution $\absolutevalue{\psi_{{\bf\it n}}}^{2}$, then the acceptance probability becomes $\min\{1,|\psi_{{\bf\it\theta}}({\bf\it n}_{i+1})/\psi_{{\bf\it\theta}}({\bf\it n}_{i})|^{2}\}$ and the state build through eq. 5 inherits the RBM state problems (data not shown). is the probability associated with the grand canonical ensemble (note that the term $\mu\expectationvalue*{\hat{N}}$ has already been introduced in the Hamiltonian). This strategy was used to generate a sample $\mathcal{M}^{\prime}$ of up to 2048 Fock states yielding a state $\ket*{\tilde{\psi}_{{\bf\it\theta}}}$ for every point in the phase diagram. The overlap between the exact ground state $\ket{\psi_{\text{exact}}}$ and the sampled RBM ground state $\ket*{\tilde{\psi}_{{\bf\it\theta}}}$ is shown in fig. 6 for the three studied scenarios.

Comparing fig. 5 with fig. 6, it is seen that the retrieved sampled state “cleans” the RBM state. In fact, the Fock states whose probability amplitudes were badly learned were removed, and only the Fock states relevant for the actual ground state were left. Despite this cleaning, the larger the Hilbert space, the more states have to be sampled to consider all relevant Fock states, as it is seen that 2048 Fock states are insufficient to capture the ground state in the case of 8 sites shown in fig. 6(b). Nonetheless, note that the low overlap in the arc from fig. 5(b) almost completely disappears after the cleaning, cf. fig. 6(b).

On the other hand, an important feature of fig. 6(b) is that the overlap of the sampled RBM state diminishes as $t$ gets larger. This happens because larger values of $t/U$ imply larger delocalization of the ground wavefunction, thus, involving more Fock states. Moreover, the Hilbert space size for 8 sites consists of $5^{8}=390625$ Fock states, which is why only sampling 2048 Fock states results in poor representations of the ground state, especially in the SF phase. Sampling more Fock states eventually reconstructs the exact ground state with very high overlap, except at the boundary between the MI and SF phases (data not shown).

The phase diagram of the BH model can be reconstructed by measuring quantities in the sampled RBM ground state that exhibit the phase transition. As mentioned before, we chose the variance of the local number operator, in particular, of the first site. In fig. 7(a) and (b), the order parameter measured with the exact ground state is shown for 5 and 8 sites, respectively. It is seen that in the Mott insulator phase, the variance of the number of bosons in the first lattice site is near to 0, but not exactly 0 because of finite size effects. On the other hand, fig. 7(d) and (e) show the order parameter for the sampled RBM ground state with 2048 Fock states for the cases of 8 and 20 hidden neurons for 5 sites, which show excellent agreement with their exact counterpart fig. 7(a). However, at the boundary between the MI and SF phases, there are absolute errors that could be as high as 0.15, which come from the difficulty of learning the ground state near the phase transition boundary. In spite of these differences at the boundaries, it is clear that the learned RBM ground state mimics the re-entrance found in finite 1D chains (Buonsante and Vezzani, 2007; Park et al., 2004). Finally, fig. 7(c) shows the order parameter for the exact RBM ground state for 8 sites and 11 hidden neurons. Figure 7(f) also shows the order parameter for 8 sites and 11 hidden neurons but with a different sampling limit. Instead of fixing a number of Fock states to be sampled (2048 and 4096 were insufficient, data not shown), we build the sample $\mathcal{M}^{\prime}$ by accepting states with the Metropolis-Hastings algorithm until there is a run of 400 consecutive state proposals that do not raise a new accepted state into $\mathcal{M}^{\prime}$. In this case, the state is cleaned (note that the arc of badly learned order parameter almost completely disappears). A clear indicator of the number of Fock states to represent the ground state emerges: for very low values of $t$, within the first Mott lobe, only one sampled state ($\ket{1,1,1,1,1}$) is needed to reproduce the exact ground state; in the SF phase, up to 30000 states are sampled before hitting a 400 streak of no new accepted states into the sample.

Other quantities can be used as an order parameter, which more explicitly relate to quantum correlations such as entanglement. For instance, a partial trace carried out over all degrees of freedom except the first site yields a reduced density matrix $\rho_{1}(t,\mu)$, from which the linear entropy $S(t,\mu)=1-\Tr{\rho^{2}_{1}(t,\mu)}$ can be measured (the von Neumann entanglement entropy can also be used, e.g., Ejima et al. (2012)). Figure 8 shows the absolute errors between the exact and the sampled RBM ground states for the three studied scenarios, where the errors at the MI-SF boundary become large.

In all the studied scenarios, there are problems for learning the ground state at the MI-SF boundaries, as well as at the mini-plateaus boundaries within the SF phase. In order to understand the differences between the learned ground state and the one obtained through exact diagonalization, we performed a study of the composition of the ground states for 5 sites and 8 hidden neurons. For that reason, we examined the probability amplitudes of the ground state (both RBM learned and exact) at 11 different points in the $t$-$\mu$ space, fixing $t=0.1$, and with $\mu$ at the middle and border of each plateau in the MI and SF phases. We also examined the RBM and exact ground states for very small $t$ at the middle of the first Mott lobe, as indicated by the black dot in fig. 9, where the ground state for the RBM and for exact diagonalization was $\ket{1,1,1,1,1}$ with an associated probability amplitude of 99.99%, as expected.

Note that the BH chain is invariant (up to a phase) under displacements and inversions, i.e. there are displacement and inversion operators that act as follows on Fock states: $\hat{T}\ket{n_{1},\ldots,n_{N-1},n_{N}}=e^{i\phi}\ket{n_{N},n_{1},\ldots,n_{N-1}}$ for displacement, and $\hat{I}\ket{n_{1},n_{2},\ldots,n_{N-1},n_{N}}=e^{i\varphi}\ket{n_{N},n_{N-1},\ldots,n_{2},n_{1}}$ for inversion. Therefore, in order to perform a tomography, we must take into account that all Fock states that belong to the same rung defined by displacement and inversion operations are equivalent. Now, each Fock state can be brought to a canonical Fock state through a consecutive application of displacement and inversion operators onto the original Fock state. This canonical Fock state is selected as the lexicographically smallest one, after every possible application of displacement and inversion operators, as in (Choo et al., 2018). In fig. 9, we show the probability amplitude distribution of the Fock state manifold rungs for the exact and the RBM ground states. More explicitly, the exact ground state is

where $i$ indexes the rungs, and $R_{i}$ is the $i$-th rung. The probability amplitude corresponding to a rung is, therefore, the sum of the probability amplitudes of all of its Fock states, and the bars from fig. 9 show those rungs’ probability amplitudes. Reading the plot from right to left, i.e. starting with the smallest value of $\mu$, we see that within the first Mott lobe, both the RBM and the exact ground states show very similar distributions over three rungs of the one-filling manifold: $\ket{1,1,1,1,1}$, $\ket{0,1,1,1,2}$ and $\ket{0,1,1,2,1}$. Increasing $\mu$ up to the first MI-SF boundary (slightly within the SF phase), we see that the probability amplitude distribution starts to differ between the exact and the RBM ground states, and the RBM struggles to identify if the ground state is now in an excitation manifold above the one-filling manifold, represented by the rung $\ket{1,1,1,1,2}$, or in the MI phase which is represented by the rung $\ket{1,1,1,1,1}$. Increasing again $\mu$ in order to be at the middle of the first plateau within the SF phase shows that both the RBM and the exact ground states have similar probability amplitude distributions, even though the RBM assigns a probability to other rungs (not shown in the plot because their contribution is less than 0.01). It is clear that this first plateau is mostly represented by the rung $\ket{1,1,1,1,2}$ which is in the 6-th excitation manifold. If we continue to increase $\mu$ we see that near the boundaries between the MI-SF phases and between the plateaus within the SF phase, the RBM and the exact ground states exhibit differences in the probability amplitude distributions. On the contrary, in the middle of those plateaus, the probability amplitude distributions from both RBM and exact ground states are very similar (which is also seen in fig. 5(a)). Moreover, each time the tomography moves onto a new plateau (with a higher value of $\mu$), the excitation manifold increases by one, up to the 10-th excitation manifold, which corresponds to the two-filling manifold at the second Mott lobe, characterized by the rung $\ket{2,2,2,2,2}$ as well as the rungs $\ket{1,2,2,2,3}$ and $\ket{1,2,3,1,3}$.