Ground state phase diagram of the one-dimensional Bose-Hubbard model from restricted Boltzmann machines

TL;DR

This paper demonstrates that restricted Boltzmann machines can effectively model the ground state phase diagram of the one-dimensional Bose-Hubbard model, capturing phase boundaries between Mott insulator and superfluid phases.

Contribution

It introduces a novel application of restricted Boltzmann machines to bosonic lattice models using NetKet, expanding their applicability in quantum many-body physics.

Findings

Successfully mapped the phase diagram using RBMs

Results agree qualitatively with exact diagonalization

Identified phase boundaries between Mott insulator and superfluid phases

Abstract

Motivated by recent advances in the representation of ground state wavefunctions of quantum many-body systems using restricted Boltzmann machines as variational ansatz, we utilize an open-source platform for constructing such ansatz called NetKet to explore the extent of applicability of restricted Boltzmann machines to bosonic lattice models. Within NetKet, we design and train these machines for the one-dimensional Bose-Hubbard model through a Monte Carlo sampling of the Fock space. We vary parameters such as the strength of the onsite repulsion, the chemical potential, the system size and the maximum site occupancy and use converged equations of state to identify phase boundaries between the Mott insulating and superfluid phases. We compare the average density and the energy to results from exact diagonalization and map out the ground state phase diagram, which agrees qualitatively with previous finding obtained through conventional means.

Figures

Equations

$$H=-t\sum _{⟨ij⟩}{b}_i^{†}{b}_j+\frac{U}{2}\sum _i{n}_i(n_i-1)-\mu \sum _i{n}_i,$$

$$|\psi ⟩=\sum _n{\psi }_n|n⟩,$$

$${\psi }_n=\sum _{\{h\}}{e}^{{\sum }_i{b}_i{\sigma }_i+{\sum }_j{c}_j{h}_j+{\sum }_{ij}{W}_{ij}{\sigma }_i{h}_j}.$$

$$\sum _{\{h\}}{e}^{{\sum }_i{b}_i{\sigma }_i+{\sum }_j{c}_j{h}_j+{\sum }_{ij}{W}_{ij}{\sigma }_i{h}_j}$$

$${\psi }_n=e^{{\sum }_i{b}_i{\sigma }_i+{\sum }_j\log \left(1+e^{c_j+{\sum }_i{W}_{ij}{\sigma }_i}\right)}.$$

$$E=\frac{\text{expectationvalue}\psi |H|\psi }{\text{innerproduct}\psi \psi }.$$

Full text

1 Introduction

In recent years, artificial neural networks have been employed to study the Bose-Hubbard model [1, 2], which describes the many-body systems of interacting bosonic particles on lattices [3]. In Ref. [1], the author approaches the problem using feedforward artificial neural networks where the inputs are Fock states (see below) and the output layer, consisting of two nodes, provides the $Re$ and $Im$ parts of the log of the wavefunction. The networks are trained following conventional supervised learning algorithms for feedforward networks with the cost function taken to be the ground state energy. Ref. [2] expands on this idea and examines the role multiple hidden layers or convolutional layers play in the accuracy of the output and the efficiency of the algorithm.

Restricted Boltzmann machines (RBMs) have also emerged as useful tools among artificial neural networks for providing a superior ansatz for representing quantum wave functions [4, 5, 6, 7, 8], or mimicking thermodynamics of classical and quantum many-body systems [9, 10, 11]. For example, in Ref. [4], it was shown that RBMs lead to lower variational ground state energy for classical and quantum magnetic models than the state-of-the-art ansatz. This in turn motivated several other studies in which the RBMs were used to represent topological [5, 6] or chiral [8] states, draw equivalences between RBMs and tensor network states [7], and to accelerate Monte Carlo simulations [12], to name a few.

Here, we use RBMs to represent ground state wavefunctions of the one-dimensional Bose-Hubbard model and study the extent to which its equilibrium properties can be captured in this architecture for a range of model parameters. It has been shown that multi-valued RBMs, such as the one we have used here, can represent a wide class of many-body entangled states efficiently [13], and to the best of our knowledge, our results are the first to demonstrate that. We perform the training using Monte Carlo sampling with the same goal of minimizing the ground state energy as in Refs. [1, 2]. We employ NetKet [14], an open-source platform for solving quantum many-body problems using artificial intelligence. We show that RBMs can be trained to represent the ground state of the bosonic system with a good degree of accuracy. We find that the energy and density both in the Mott insulating and superfluid phases rapidly converge to final values during the training. By calculating equations of state for different values of the interaction strength in the grand canonical ensemble, we map out the phase diagram of the model and show that it agrees well with one obtained in the past directly from quantum Monte Carlo simulations on larger systems.

2 The Bose-Hubbard Model

The Bose-Hubbard Hamiltonian is expressed as

[TABLE]

where $b_{i}^{\dagger}$ and $b_{j}$ are the bosonic creation and annihilation operators, respectively, $\left<ij\right>$ indicates that $i$ and $j$ are nearest neighbors, $t$ is the corresponding hopping amplitude, which we set to unity to serve as the unit of energy throughout the paper, $U$ is the strength of the onsite repulsion, $n_{i}=b_{i}^{\dagger}b_{i}$ is the site occupation operator, and $\mu$ is the chemical potential. In this study, we consider the periodic one-dimensional (1D) geometry.

For large enough interactions, the ground state of the system is the interaction-driven Mott insulating phase while it remains a superfluid for smaller interactions. This was shown in a seminal paper by Fisher et al. [3] where the universality class of the model was explored and a mean-field solution was presented. Later, the model received much attention in the 1990’s, mostly confirming numerically the transition in different dimensions [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. The Mott insulator (MI) to superfluid (SF) transition was also observed in experiments with cold atoms in optical lattices [26].

3 Restricted Boltzmann Machine

3.1 The Architecture

The RBM is a two-layer fully-connected artificial neural network. The architecture is shown in Fig. 1. It has one visible layer, and one hidden layer. Each layer consists of several nodes with an external bias parameter associated with them ($b_{i}$ for the visible and $c_{i}$ for the hidden layer). Each node in a layer is connected to all the nodes in the other layer and a weight parameter ($W_{ij}$) is associated to every connection.

The RBM lends itself well to representing ground state wavefunctions of the Bose-Hubbard model because a probability distribution of possible Fock states of the system can then be calculated directly from the RBM parameters after training. This is similar to representing thermodynamics of a many-body system using an RBM where a probability distribution of microstates can be extracted [10].

Within NetKet, we define a wave function ansatz for the $q=0$ momentum sector using a RBM with spin-1/2 hidden nodes and permutation symmetry. The input layer can be visualized as a two-dimensional $L\times(N_{max}+1)$ array of binary numbers, where each row represents a site, and a 1 in column number $n$ indicates that there are $n$ bosons at that site. This so-called “one-hot” encoding of site occupations is not unique. For example, an alternative would be to use $L$ non-binary input nodes whose values can vary from [math] to $N_{max}$ instead. One advantage of the one-hot representation over the latter is that the number of network parameters is much larger and therefore, the training will presumably be more refined. We choose $L$ and $N_{max}$ values as large as 7 and 5, respectively, and a hidden layer with twice as many nodes as in the visible layer.

3.2 The Training

We write the ground state wave function as a linear combination of the Fock state basis

[TABLE]

where $\ket{n}=\ket{n_{1},n_{2},\dots,n_{L}}$ are the particle number Fock states with site occupancies $n_{i}\in[0,N_{max}]$. The goal is to create a RBM that can stochastically represent $\ket{\psi}$ by producing the correct probability, $|\psi_{n}|^{2}$ for every input Fock state $\ket{n}$ through the input layer. For the Bose-Hubbard model, the coefficients $\psi_{n}$ are known to be real and positive and are written for the RBM as

[TABLE]

where $\sigma_{i}$ and $h_{i}$ denote the input node and hidden node binary values, respectively.

The hidden degrees of freedom can be integrated out exactly for the simple architecture shown in Fig. 1,

[TABLE]

which results in [4, 9]

[TABLE]

The goal in the training of the RBM is to optimize its weights and biases in order to arrive at a probability distribution for $\psi_{n}$ so that $\ket{\psi}$ represents the ground state. This is accomplished through the minimization of the expectation value of the Hamiltonian,

[TABLE]

NetKet adopts a Metropolis algorithm for importance sampling of $\ket{n}$’s during which $E$ is minimized also using stochastic techniques. Since we are working in the grand canonical ensemble, we choose local moves in which a site is picked at random and its occupation is proposed to change to a value in the range $[0,N_{max}]$. The move is accepted with the probability min$[1,|\psi_{n^{\prime}}|^{2}/|\psi_{n}|^{2}]$, where $n^{\prime}$ and $n$ are the new and old states, respectively. This guarantees that after the training is completed, the probability distribution is $P_{n}=|\psi_{n}|^{2}/\sum_{m}|\psi_{m}|^{2}=|\psi_{n}|^{2}/\innerproduct{\psi}{\psi}$. The minimization of energy is done using the Stochastic Reconfiguration algorithm [27] with AdaMax optimizer in NetKet and the default parameters. Our sampler performs 1000 sweeps before sampling the probability distribution and updating the RBM parameters. More details about the method and the training parameters can be found in Ref. [14].

4 Results

In Figs. 2 and 3, we show the average energy, $E$, and the average density, $\rho$, as the training of our RBM progresses for two sets of values for the interaction strength and the chemical potential. The system with $\mu/U=0.56$ (blue lines) is a Mott insulator with a density of one particle per site and the system with $\mu/U=1.00$ (orange lines) is a superfluid for the system size we have considered ($L=6$). Here, $N_{max}=5$ is much larger than the average density. We observe that the energies and densities quickly converge to their final values regardless of whether the system is in the MI or SF phase. However, we find that the average density has larger fluctuations over the training steps inside the SF phase while it is pinned to unity when the system is in the first Mott lobe.

For comparison, we also show in Figs. 2 and 3 results from exact diagonalization as horizontal dashed and dotted lines. We find that the error in both the energy and the density is larger in the superfluid phase. The fluctuations in the density seem to extend to very long training times and point to the limitations of the method and the fact that more complex networks or better optimization may be needed to capture the physics in the superfluid phase.

Taking long-time values of the converged density as we vary the chemical potential at a fixed interaction strength, we obtain the equation of state. Figure 4 shows this property for $U/t=15$ and several different combinations of $L$ and $N_{max}$. We find the expected behavior [3, 15] where the density is pinned at integer numbers of bosons per site as we cross multiple SF and MI phase boundaries. Two density plateaus corresponding to the first and second Mott lobes are clearly visible in this plot. Increasing the system size from $4$ to $7$ seems to lead to smoother curves, signaling better convergence during the training, but no appreciable change in the onsets of the Mott regions. Moreover, the results do not change significantly by increasing $N_{max}$ from $4$ to $5$. These observations do not mean that our results are not suffering from finite size effects. Much larger clusters have been shown to give significantly different results in conventional treatments of the model [15, 28].

By training different RBMs as we sweep the chemical potential at increasing values of the interaction strength, we are able to map out the ground state phase diagram of the model in the interaction-chemical potential space. The findings are summarized in Fig. 5. We extend the results up to the second Mott lobe. They qualitatively agree with those obtained using Monte Carlo algorithms and much larger system sizes (see e.g. Ref. [28]). As with other techniques, the accurate locating of the phase boundary near the tip of the Mott lobes are the most challenging due to the Kosterlitz-Thouless nature of the transition [29].

5 Summary

Using the open-source package NetKet, recently developed for the study of quantum many-body systems through machine learning algorithms, we design and train RBMs to represent the ground state wavefunction of the 1D Bose-Hubbard model on small chains. We demonstrate the degree of applicability of the RBM ansatz in various regions of the parameter space as MI or SF phases set in. We show that the optimization of free parameters of the RBMs during training in order to minimize the ground state energy leads to physical properties, including expectation values of the Hamiltonian and a phase diagram, that agree with known exact results from conventional treatments of the model.

6 Acknowledgements

KM acknowledges support from Undergraduate Research Grants at San Jose State University. EK acknowledges support from the National Science Foundation (NSF) under Grant No. DMR-1609560.

7 References