Neural Network Solutions of Bosonic Quantum Systems in One Dimension
Neural networks have been proposed as efficient numerical wavefunction
ansatze which can be used to variationally search a wide range of functional
forms for ground state solutions. These neural network methods are also
advantageous in that more variational parameters and system degrees of freedom
can be easily added. We benchmark the methodology by using neural networks to
study several different integrable bosonic quantum systems in one dimension and
compare our results to the exact solutions. While testing the scalability of
the procedure to systems with many particles, we also introduce using symmetric
function inputs to the neural network to enforce exchange symmetries of
indistinguishable particles.
Paulo F. Bedaque, Hersh Kumar, Andy Sheng
# Neural Network Solutions of Bosonic Quantum Systems in One Dimension

###### Abstract

Neural networks have been proposed as efficient numerical wavefunction ansatze which can be used to variationally search a wide range of functional forms for ground state solutions. These neural network methods are also advantageous in that more variational parameters and system degrees of freedom can be easily added. We benchmark the methodology by using neural networks to study several different integrable bosonic quantum systems in one dimension and compare our results to the exact solutions. While testing the scalability of the procedure to systems with many particles, we also introduce using symmetric function inputs to the neural network to enforce exchange symmetries of indistinguishable particles.

## I Introduction

Solving many-body quantum mechanical systems involves obtaining solutions to a high-dimensional partial differential equation or diagonalizing exponentially large Hamiltonians, both of which pose significant challenges. Thus, analytically solved models are few and far between. [1; 2; 3; 4; 5; 6; 7].

Various methods have been developed to obtain approximate solutions to complex quantum systems. For one such method called Variational Monte Carlo, a variational ansatz for the ground state wavefunction is assumed and an estimate for the ground state energy is obtained by minimizing the energy of the system over the parameters of the chosen approximate wavefunction form. The use of trial wavefunctions has been used to study a broad range of many-body systems, from Fermi and Bose gases to quantum chemistry [8; 9; 10; 11; 12; 13; 14; 15; 16; 17]. In particular, the Slater-Jastrow [18] wavefunction ansatz has been employed widely [19; 20; 21] because, by construction, it captures correlations between interacting particles and incorporates the correct intuitive functional form of the physical system. Tensor networks and density matrix renormalization group methods [22] were shown to efficiently encode low-lying eigenstates of local, gapped Hamiltonian spin systems using only a polynomial number of parameters. This eliminates the need to search through an exponentially-large Hilbert space.

More recently, established machine learning techniques such as neural networks have found their way into quantum physics [23; 24]. In the context of studying ground states of quantum systems, neural networks conveniently provide a numerical wavefunction ansatz for which one can vary to minimize the system energy cost function. The minimization can be done efficiently using established machine learning methods of backpropagation and gradient descent. Another appeal of using neural networks is that they were shown to be universal function approximators [25; 26] and thus can be used to describe any general wavefunction. As a consequence, the number of variational parameters is not intrinsically tied to the function form (i.e. Jastrow) of the ansatz and can be increased easily.

Carleo and Troyer [27] showed that a type of neural network called restricted Boltzmann machines were sufficient to describe ground state wavefunctions of spin systems with a reasonable number of parameters. Later, restricted Boltzmann machines was found to map directly to matrix product state approximations of spin states [28; 29]. For particles interacting via continuous-space potentials, several systems of indistinguishable particles have also been studied using neural networks including the Calogero-Sutherland model [30], nuclear models [31; 32; 33], and chemical systems [34; 35].

Quantum systems involving indistinguishable particles must have ground state wavefunctions which satisfy position exchange (anti)symmetries. This poses a challenge in which the constructed neural network trial wavefunction should satisfy these constraints. In one spatial dimension, Bose symmetry can be satisfied by feeding in ordered coordinates of the bosons into the neural network [30]; however, this method is not easily generalizable to systems of higher dimensions. Pfau et. al. [35] developed the FermiNet neural network structure to impose fermionic antisymmetry to study electron-electron and electron-ion interactions in atoms. The FermiNet structure constructs the most general fermionic wavefunction space for which one can search for the ground state through by taking Slater determinants of permutation-equivariant functions.
The permutation-equivariance of such functions are enforced through permutation-equivariant neural network layers.

In this paper, we instead propose the use of carefully selected symmetric inputs to enforce the Bose symmetry of the neural network ansatz. By using symmetric functions as inputs to the neural network, we are able to construct the most general Bose symmetric function as we increase the number of nodes and layers; this gives us full flexibility in changing the architecture of the network. We benchmark the use of symmetric functions by computing ground state energies and wavefunctions of several one-dimensional, exactly-solvable quantum systems and comparing both quantitative and qualitative results. Our methodology is not restricted to studying exactly-solvable systems and we also explore systems in which no analytical solution is known. We further demonstrate the scalability and versatility of our neural network ansatz by calculating ground states of systems with dozens of particles.

In Sec. II, we outline the various models for which we study and compare numerical results to. In Sec. III, we explicitly construct the Bose symmetric neural network functions and outline our procedure for computing ground states. Sec. IV discusses our results; we find that our neural network ansatz is able to accurately compute ground state energies and features of the many-body systems considered over their various phases. Then we conclude with Sec. V by exploring prospects of using symmetric function constructions to write general bosonic and fermionic trial wavefunctions in a higher number of spatial dimensions.

## II Models

To show the validity of our symmetric input neural network solutions, we apply our neural-network-based methodologies to solve for the ground state wavefunctions and energies for several different one-dimensional quantum systems. We consider both a system of cold bosons in a harmonic trap and a system of trapped bosons with short and long-range interactions. Both quantum systems are inspired by the Lieb-Liniger model Lieb and Liniger (1963); Lieb and Liniger (1963) -- a homogeneous gas of indistinguishable bosons interacting via a contact delta-function potential which was shown to be exactly solvable in one dimension through the use of Bethe ansatz Lieb and Liniger (1964); Lieb and Liniger (1965); Lieb and Liniger (1966). These Lieb-Liniger-inspired models have rich phase structure and have been studied through the lenses of both finding exact solutions and performing numerical analyses Lieb and Liniger (1964); Lieb and Liniger (1964); Lieb and Liniger (1965); Lieb and Liniger (1965); Lieb (1966); Lieb (1966). Thus, they provide good testbeds for which to compare both quantitative and qualitative ground state results to.

### Cold Bosons in a Harmonic Trap

The first Hamiltonian which we consider describes spinless bosons in a harmonic trap. The interactions of the condensate are approximated by a contact potential. Using units where \(\hbar=1\), the Hamiltonian is

\[\hat{H}=\sum_{i=1}^{N}\left(-\frac{1}{2m}\frac{\partial^{2}}{\partial x_{i}^{ 2}}+\frac{1}{2}m\omega^{2}x_{i}^{2}\right)+\sum_{i<j}^{N}g\delta(x_{i}-x_{j}) \tag{1}\]

In general, we take \(m=\omega=1\). The interaction is repulsive when \(g>0\) and attractive when \(g<0\). When \(g>0\), the interaction strength \(g\) is proportional to the two-boson s-wave scattering length Lieb and Liniger (1964). Although the addition of the harmonic trap causes the system to no longer be analytically solvable Lieb and Liniger (1964) for the general case of many particles, the two-boson case has been solved exactly Lieb and Liniger (1964); Lieb and Liniger (1965).

### Trapped Bosons with Short and Long-Range Interactions

Recently, another quantum system has been shown to be exactly solvable for a general case of \(N\) bosons Lieb and Liniger (1964). The model consists again of particles in harmonic confinement interacting with a delta-function potential of strength \(g\). Additionally, a long-range linear interaction potential of strength \(\sigma\) is added, which, in one dimension, corresponds to either Coulomb repulsion (\(\sigma<0\)) or gravitational attraction (\(\sigma>0\)).

\[\begin{split}\hat{H}&=\sum_{i=1}^{N}\left(-\frac{1 }{2m}\frac{\partial^{2}}{\partial x_{i}^{2}}+\frac{1}{2}m\omega^{2}x_{i}^{2} \right)\\ &+\sum_{i<j}^{N}\left(g\delta(x_{i}-x_{j})+\sigma|x_{i}-x_{j}| \right)\end{split} \tag{2}\]

Again, we take \(m=\omega=1\). The model was found to be integrable when \(\sigma=-m\omega g/2\)1. In that case, \(g>0\) corresponds to both interaction terms being repulsive and \(g<0\) corresponds to both terms being
attractive. The coupling strength is again related to the one-dimensional s-wave scattering length \(a_{s}=-2/(mg)\). The exact expression for the ground state energy of the Hamiltonian (2) was found to be

\[E_{0}=\frac{N\omega}{2}-mg^{2}\frac{N(N^{2}-1)}{24} \tag{3}\]

with a corresponding (non-normalized) ground state wave function

\[\psi_{0}(x_{1},...x_{N})=\prod_{i<j}e^{-|x_{i}-x_{j}|/a_{s}}\prod_{i}e^{-x_{i}^ {2}/(2a_{\rm ho}^{2})} \tag{4}\]

where \(a_{\rm ho}=\sqrt{1/(m\omega)}\) denotes the harmonic oscillator length.

When \(g>0\), the ground state wavefunction is maximized at \(|x_{i}-x_{j}|=Na_{\rm ho}^{2}/a_{s}\) and the bosons tend to spread out with equal spacing about the center of the harmonic trap. When a weak repulsion is turned on, the system forms an incompressible liquid of flat density [44; 45]. In the limit of strong repulsion \(g\gg 1\) the probability of overlap between any two particles drops to zero and the system forms a Wigner crystal [46]. In the other regime, when \(g<0\), the bosons collapse into the center of the trap. In the limit of \(g\rightarrow-\infty\), the ground state wavefunction of (2) becomes

\[\psi_{0}(x_{1},...,x_{N})=\prod_{i<j}e^{-|x_{i}-x_{j}|/a_{s}} \tag{5}\]

resembling a McGuire bound-state solution [39] describing a bright soliton [47]. It is surprising that despite the fact that different signs of the interaction strength \(g\) yield such different behaviors, the expression (3) for the ground state energy of the system implies that the ground state energy does not depend on the sign of \(g\).

## III Methods

Given a Hamiltonian describing a quantum system of \(N\) indistinguishable bosons, we would like to solve the time independent Schrodinger equation for the ground state energy and wavefunction. To do this, we employ the variational principle; any arbitrary wavefunction \(\psi\) must have an expectation value of the energy at least as large as the ground state energy.

\[\langle E\rangle_{0}\leq\langle E(\theta)\rangle=\frac{\int dX\:\psi^{\dagger }(X,\theta)\hat{H}\psi(X,\theta)}{\int dX\:\psi^{\dagger}(X,\theta)\psi(X, \theta)} \tag{6}\]

Here \(X\) stands for the positions of the \(N\) bosons. Typically, one begins with an ansatz for the wavefunction \(\psi(X,\theta)\), where the \(\theta\) are its variational parameters and then minimize \(\langle E(\theta)\rangle\) with respect to \(\theta\) to obtain an estimate for the ground state of the quantum system.

### Neural Network Wavefunction Antsatze

The quality of the estimate for the ground state energy and wavefunction depends both on the number of variational parameters one uses as well as the particular functional form of \(\psi\). Although one could successfully make well-posed guesses as to the correct functional form of the ground state wavefunction based on expectations on the underlying physics of the system, we would like to more systematically search the space of potential ground state wavefunctions. In this context, the application of neural networks is powerful.

#### iii.1.1 Neural Networks

We define a neural network with \(L\) layers as a function which takes \(n_{1}=N\) inputs and produces \(n_{L}=M\) outputs -- defined in the following way:

\[\begin{split}& I_{i+1}=O_{i}=f(\mathbf{W}_{i}\vec{I}_{i}+\vec{b}_{ i}),\ 1\leq i\leq L-1\\ & O_{L}=\mathbf{W}_{L}\vec{I}_{L}+\vec{b}_{L}\end{split} \tag{7}\]

The inputs \(I_{i}\) at layer \(i\) are multiplied by an \(n_{i+1}\times n_{i}\) matrix of weights \(\mathbf{W}_{i}\)2. A bias vector \(\vec{b}_{i}\) of length \(n_{i+1}\) is added before applying an activation function \(f\) to each resulting vector element. The output \(O_{i}\) of the \(i\)-th layer then becomes the input for the \((i+1)\)-th layer, where the same procedure is repeated with a different set of weights and biases. While in general, a different activation function can be used at each layer, we choose to keep the same function throughout. To ensure that the neural network produces an output whose range is over all reals, an activation function is _not_ used in the final output layer.

Footnote 2: For clarity, \(n_{1}=N\) and \(n_{L}=M\); \(n_{i}\) is the number of nodes at layer \(i\).

Neural networks are useful in the context of writing down a variational ansatz for the wavefunction because multilayer neural networks have been shown to be universal function approximators [25; 26]. The advantage is that, without needing to change the architecture of the neural network
function, we can explore a general space of functions to find the ground state wavefunction by varying the parameters \(\theta=\{\mathbf{W}_{i},\vec{b}_{i}\}\). The efficiency of the functional expressivity for many-body wavefunctions, however, is not clear and is a property we seek to explore. We denote our ansatz neural network function as \(\mathcal{A}(I_{1},...,I_{N})\); it takes in \(N\) inputs and returns a single output (\(M=1\)).

#### ii.1.2 Enforcing Ground State Symmetries

A neural network ansatz for the ground state wavefunction of a system of \(N\) indistinguishable bosons should satisfy the following two constraints. Firstly, the wavefunction must obey Bose exchange symmetry.

\[\psi(...,x_{i},...,x_{j},...)=\psi(...,x_{j},...,x_{i},...) \tag{8}\]

where the \(x_{n}\) are the positions of the bosons. Thus we must construct our neural network such that it is fully symmetric under exchange of particle positions. To do this, we look for a bijection between the set of particle positions \(\{x_{i}\}\) and a set of functions of \(\{x_{i}\}\) which are symmetric under \(x_{i}\leftrightarrow x_{j}\). The simplest such bijection generates the set of elementary symmetric polynomials

\[e_{i}=\sum_{1\leq j_{1}<j_{2}<...<j_{i}\leq N}x_{j_{1}}x_{j_{2}}...x_{j_{i}} \tag{9}\]

For instance, for \(N=3\):

\[\begin{split} e_{1}&=x_{1}+x_{2}+x_{3},\\ e_{2}&=x_{1}x_{2}+x_{1}x_{3}+x_{2}x_{3},\\ e_{3}&=x_{1}x_{2}x_{3}.\end{split} \tag{10}\]

We can see this by noting that the \(e_{i}\) are the coefficients of a degree-\(N\) polynomial over an auxiliary variable \(t\) whose \(N\) roots are the positions \(x_{i}\) -- namely

\[f(t)=(t-x_{1})(t-x_{2})...(t-x_{N}) \tag{11}\]

For ease of implementation however, we note that the Newton-Girard identities define a further bijection between the set of \(\{e_{i}\}\) to the set of power-sums

\[\xi_{i}=\sum_{k=1}^{N}\,x_{k}^{i} \tag{12}\]

Since the set of \(\{\xi_{i}\}\) are also fully symmetric under exchange of \(x_{i}\), we can construct a Bose symmetric neural network function by feeding in the \(\{\xi_{i}\}\) as inputs (\(I_{i}=\xi_{i}\)). In practice, for large \(N\), \(|\xi_{i}|\) can easily become extremely large; thus, we normalize the initial particle positions by an approximate width \(w\) of the position-space ground state wavefunction 3

Footnote 3: For example, if the system were a harmonic oscillator with \(m=\omega=1\), we choose \(w\approx 2\) such that the majority of the nonzero region of the Gaussian ground state is contained within the selected width.

\[\tilde{\xi}_{i}=\sum_{k=1}^{N}\left(\frac{x_{k}}{w}\right)^{i} \tag{13}\]

such that the \(\tilde{\xi}_{i}\) are at most \(\mathcal{O}(1)\) throughout the search for the ground state wavefunction.

\[\psi(x_{1},...,x_{N})\sim\mathcal{A}(\tilde{\xi}_{1},...,\tilde{\xi}_{N}) \tag{14}\]

The \(\{\tilde{\xi}_{i}\}\) retains information about the original particle coordinates but loses information about their ordering, as desired.

Secondly, the ground state wavefunction of a system of indistinguishable bosons can be taken to be real and positive. Thus, we constrain our neural network weights and biases to be real-valued. We can also, without loss of generality, write the wavefunction ansatz as

\[\psi(x_{1},...,x_{N})=e^{-\mathcal{A}(\tilde{\xi}_{1},...,\tilde{\xi}_{N})} \cdot e^{-\Omega\sum_{i=1}^{N}\,x_{i}^{2}} \tag{15}\]

Since we are looking for bound ground state solutions of the Hamiltonians models we consider, we want to search the space of wavefunctions which vanish at infinity. We multiply our neural network function by the Gaussian factor to ensure that this property is satisfied even when the neural network function itself may be relatively flat over its inputs for a given set of weights and biases. Because the bosons in the quantum systems we study are confined in a harmonic trap with \(m=\omega=1\), we choose \(\Omega=0.5\). The neural network function \(\mathcal{A}\), then, is purely responsible for capturing the complexities of particle interactions. At zero inter-particle interaction strength, \(\mathcal{A}\) should minimize trivially to zero everywhere in the ground state.

### Neural Network Architecture

To search for the ground state energy and wavefunction of a quantum system, we first construct a neural network wavefunction ansatz of the form in (15). For the neural network, we choose to use the continuously-differentiable CELU as our activation
function \(f\).

\[\text{CELU}(x)=\begin{cases}x&x>0\\ e^{x}-1&x\leq 0\end{cases} \tag{16}\]

Furthermore, we initialize the real-valued weights and biases of the neural network to be \(\mathcal{O}(1)\) and nonzero -- drawn randomly from a normal distribution centered around zero with standard deviation \(\sigma=1\). Using \(\langle E\rangle\) as a cost function, we perform gradient descent using the Adam optimizer [48] to approach the ground state. We use Adam step sizes of between \(10^{-3}\) and \(10^{-6}\), decreasing the step size as the minimization approaches closer to convergence. Using relatively large step sizes in certain regions of parameter space close to convergence can lead to being stuck in orbit cycles and exhibit oscillatory behavior in the cost function [49]. The other hyperparameters of the Adam optimizer describing the decay rates of the moments and ensuring numerical stability are kept at their default values: \(\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=10^{-8}\).

Because the integrals involved in obtaining \(\langle E\rangle\) and \(\langle\partial E/\partial\theta\rangle\) are difficult to compute exactly for a complicated neural network function and for many particles \(N\), we estimate the energy and the gradient using Metropolis importance sampling, drawing from the distribution \(\pi(X)\sim\psi(X,\theta)^{2}/\int dX\;\psi(X,\theta)^{2}\). However, the Hamiltonians which we study contain potential terms which involve delta functions in position. This creates a zero-overlap problem in the Monte Carlo evaluation of this integral in which none of the sampled position configurations induce a contribution to the energy or gradient from those delta function terms. We introduce a remedy in Appendix A. Although computing the gradient using Monte Carlo methods is not as precise as an analytic calculation, the stochasticity tends to aid in avoiding being stuck in local minima during the minimization [50; 51; 52].

In general, there has not been a systematic solution to choosing an optimal number of layers and nodes of a neural network for any given problem [53]. Here, we arbitrarily choose to construct our neural network function with between four to ten layers. Although a more careful analysis in choosing the number of layers has yet to be done, we decide to use more layers when solving quantum systems of more particles since we expect that more layers are needed to capture the ground state of a more complicated system. Keeping the same number of layers in solving a particular system, we systematically increase the number of nodes at each layer until we find that our prediction for the ground state energy ceases to decrease within error bars given by the finite Monte Carlo sampling. Fig. 1 shows convergence in \(\langle E(\beta)\rangle\) as we increase the total number of variational parameters (weights and biases) in the neural network. We find that an increasing number of parameters for systems of more particles is needed; however, this dependence appears to be slower than exponential in number of particles.

Finally, to perform the neural network computations and gradient descent minimizations, we utilize Google's JAX [54] package in Python. JAX has built-in automatic differentiation, just-in-time compilation speed-up, and vectorization features which make it an appealing and convenient choice to run neural network computations with.

## IV Results

To check that our variational neural network wavefunction ansatz (15) is universal enough to capture the ground state behaviors of many different quantum systems, we use our ansatz to compute the ground states in various phases of the Hamiltonian models described in Sec. II. We compare our results to the analytic solutions of the exactly solvable models. We then also explore an unsolved regime of one of these models.

Figure 1: Shows the relative error of the ground state energies for the quantum system of \(N\) bosons described by the Hamiltonian in (2) as the number of neural network parameters \(\beta\) is increased. Eight layers were used in the neural network throughout. The statistical error bars from the Monte Carlo calculations on some points are too small to be seen. \(\langle E\rangle_{0}\) here is taken to be the value of the ground state energy computed with the largest \(\beta\).
### Cold Bosons in a Harmonic Trap

Preliminarily, we first study a one-dimensional system of two indistinguishable bosons (\(N=2\)) in a harmonic trap, interacting via a contact potential (II.1). An exact expression for the ground state energy of the system was given by Busch et. al. [40]. Using our neural network ansatz, we compute the ground state energies of the system at various interaction strengths \(g\). Fig. 2 shows a comparison of the two results. As we increase the number of parameters \(\beta\) in the neural network ansatz, we find that, over the wide range of \(g\), the neural network ground state energies converge towards agreement with those of the analytic solution to within error bars. The mean values are also consistent with the exact expression to approximately one to two percent precision. Throughout the final calculations for this model, we fixed our neural network to have four layers and roughly \(5\cdot 10^{4}\) variational parameters. Although this likely uses more parameters than needed to describe the simpler two-boson system, the results give confidence that our Bose symmetric neural network ansatz, as desired, is able to capture the ground state physics of the different systems without needing to specify a particular functional form for the ansatz or adjusting the architecture of the neural network.

### Trapped Bosons with Short and Long-Range Interactions

Next, we examine the model in which the bosons interact with a long-range linear potential along with the delta-function potential whose properties are summarized in Sec. II.2. We are able to use the neural network ansatz to study both the

Figure 3: Shows ground state behaviors of two bosons interacting with short and long-range interactions described in II.2 in the exactly-solvable regime \(\sigma=-g/2\). We plot the position-space probability density; brighter colors indicate a larger value of \(\psi(x,y)^{2}\) and the bosons have a higher probability of being there. The blue lines which run through \(y=-x\) are shown to guide the eye. The top shows the system with a repulsive potential (\(g=1.5\)) while the bottom shows the system with an attractive potential (\(g=-1.0\)).

Figure 2: Shows the ground state energies of the quantum system described in II.1. We plot the results computed from our neural network ansatz (red) against the exact results (black) given by Busch et. al. [40]. The statistical error bars from the Monte Carlo integration of the energy in most cases are too small to be seen. Convergence towards the exact result is seen as the number of parameters \(\beta\) in the neural network ansatz is increased.
exactly-solvable regime of \(\sigma=-g/2\) solved by Beau et. al [7] as well as the regime \(\sigma=-g\) in which no analytic solution is known.

#### iii.1.1 Exactly-solvable Regime (\(\sigma=-g/2\))

We again begin with a simpler system of just two particles. Using our neural network ansatz, we compute the ground state energies and wavefunctions for various values of the coupling in the exactly-solvable regime \(g=-\sigma/2\) and compare the results to the analytic solution. Fig. 3 displays the ground state wavefunction of the system while Fig. 4 shows a comparison of the ground state energies. We find that the ground states given by the neural network ansatz give the same behavior as described by the exact ground state wavefunction (4). When \(g>0\) and the interactions are both repulsive, the two bosons separate an equal distance apart from the center of the harmonic trap. Alternatively, when \(g<0\) and the interactions are both attractive, the bosons collapse towards the center of the trap and form a bright soliton. Quantitatively, the ground state energies computed using the variational method also agree with the analytic expression (3) to within two-percent statistical errors. The unintuitive quadratic dependence of the ground state energy on \(g\) is also confirmed.

Having shown the ability of our neural network ansatz to accurately compute ground state wavefunctions and energies of smaller, simpler systems, we further demonstrate the versatility and generalizability of the neural network methodology by computing ground state properties of the quantum system in (2) with fifty indistinguishable bosons 4. We study both a repulsive system (\(g=0.05\)) and an attractive system (\(g=-0.03\)). Using a neural network ansatz with \(\beta\approx 2\cdot 10^{5}\), we find that the computed ground state energies fall within one percent of those predicted by

Figure 5: Shows the ground state density functions for the system of fifty bosons with \(\sigma=-g/2\). The top panel displays the normalized local density while the bottom displays the normalized pair correlation. The points are heights of the position/distance binning and are linearly interpolated to guide the eye. The dark/light red show the density functions from the analytic/neural network wavefunction respectively in the attractive case while the dark/light blue show the density function sfrom the analytic/neural network wavefunction in the repulsive case.

Figure 4: Shows a comparison of the ground state energies of the system (2) with two bosons at potential strengths \(\sigma=-g/2\). The black line shows the analytic expression for the ground state energy \(E_{0}\) from Beau et. al. (3) and the red points are ground state energies computed from the neural network ansatz.
the exact \(E_{0}\) (3); the calculations yielded for \(g=0.05\), \(E_{0}^{\rm NN}=12.083(72)\) as compared with \(E_{0}=11.984375\) and for \(g=-0.03\), \(E_{0}^{\rm NN}=20.367(76)\) as compared with \(E_{0}=20.314375\).

To comprehensively check of the ground state features of the system, we also calculate the local density profile

\[n(x)=\int dx_{2}...dx_{N}\psi_{0}(x=x_{1},x_{2},...,x_{N})^{2} \tag{17}\]

and the density-density correlation function

\[g_{2}(y)=\int dx_{1}...dx_{N}\:\delta(y-|x_{1}-x_{2}|)\psi_{0}(x_{1},...,x_{N}) ^{2} \tag{18}\]

for both the analytic \(\psi_{0}\) (4) and the ground state wavefunction given by the neural network ansatz \(\psi_{0}^{\rm NN}\)5. The density functions \(n(x)\) and \(g_{2}(y)\) are computed by generating histograms of positions \(x_{1}\) and distances \(x_{1}-x_{2}\) respectively -- both taken from Monte Carlo samples drawn from the distribution \(\pi(X)\sim\psi_{0}(X)^{2}/\int dX\psi_{0}(X)^{2}\). Fig. 5 shows the density functions from both \(\psi_{0}\) and \(\psi_{0}^{\rm NN}\) plotted on top of one another. We find close agreement between the two in both phases of the model (\(g>0\) and \(g<0\)). This gives evidence that the neural network ansatz is not only able to provide an accurate estimate for the ground state energy but can also faithfully encode information about the ground state wavefunction.

Footnote 5: The choices \(x=x_{1}\) and \(\delta(y-|x_{1}-x_{2}|)\) in \(n(x)\) and \(g_{2}(y)\) respectively are arbitrary because \(\psi_{0}\) is Bose symmetric.

#### iii.2.2 Non-integrable Regime (\(\sigma=-g\))

The variational approach using neural network ansatze is not restricted to computing integrable systems. We also apply the methodology to explore the system of bosons interacting with short and long-range interactions (2) with symmetric coupling strengths \(\sigma=-g\) for which no exact ground state has been found. We then compare our computed ground state properties to those of the exactly solvable system \(\sigma=-g/2\) (see [7]). The results show that although the ground state energies behave quantitatively very differently, the general phase structures of the two regimes are similar in nature.

We begin by looking at a system of four bosons with a strong repulsive coupling \(g=-\sigma=2.5\) and another with an attractive coupling \(g=-\sigma=-1.0\). Fig. 6 shows the local density and the pair-correlation function for both four-boson systems. Just as with the exactly-solvable model, a crystal structure is clearly seen in the case of strong repulsion. The local density shows four distinct equally spaced peaks. The pair-correlation function shows that the probability of two particles being in the same position is close to zero. Instead, the correlation function is maximal at distinct separations corresponding to \(y=\pm|x_{1}-x_{2}|=nd,n\in\mathds{Z}\) where \(d\) is the inter-particle spacing between two nearest neighbors in the crystal. In the case of attraction on the other hand, both the density profile and the pair density peak strongly at \(y=0\); again, all of the particles clump towards the center of the harmonic trap and form a bright soliton [47].

We further use our neural network ansatz to search for the ground states of systems of ten and

Figure 6: Shows the ground state density functions for the system of four bosons. The top panel displays the normalized local density while the bottom displays the normalized pair correlation. Each panel contains results for the attractive (blue) and the repulsive (red) systems.
fifty particles at various coupling strengths \(g=-\sigma\). Fig. 7 shows the computed ground state energies as a function of \(g\) while Fig. 8 shows the single-particle and pair density profiles. We find that here, the energies exhibit the same asymptotic behavior of \(E_{0}\rightarrow-\infty\) as \(|g|\rightarrow\infty\) as in the exactly-solvable model. However, the ground state energy is no longer independent of the sign of \(g\) and changes asymmetrically as the interactions in the system flip from repulsion to attraction. When \(g>0\), the energy drops more rapidly than in the exactly-solvable case while when \(g<0\), the energy of the resulting soliton initially rises before dropping again at higher attractive coupling.

The qualitative features, however, again match the behaviors of the system at \(\sigma=-g/2\). From the local density profile, we see that as we turn on a repulsive potential (\(g>0\)), the bosons begin to spread symmetrically about the center of the harmonic trap. As the repulsion becomes stronger, the particle density shows a mesa-like feature of finite width \(w\). This resembles, for example, the behavior of nuclei at saturation density. The addition of particles to the system does not change the density but instead increases the droplet size. The two-body correlation for many particles at weak repulsion increases monotonically from the maximal separation \(y=\pm w\) inwards towards \(y=0\) since more pairs of particles have smaller rather than larger separations. In the vicinity of zero separation, the probability of finding two particles in the same location is decreased and we expect a slight dip in the pair-correlation function. This is clearly seen in the system with ten particles where the repulsion was stronger; the dip is not seen in the fifty particle system because the weaker re

Figure 8: Shows the ground state density functions for systems of ten and fifty bosons. In each panel, the (light) blue shows the (repulsive) attractive system for \(N=10\) and the (light) red shows the (repulsive) attractive system for \(N=50\).

Figure 7: Shows a comparison of the ground state energies of the Hamiltonian (2) with N = 10 and 50 bosons. The red points are energies computed from the neural network ansatz at the non-integrable regime \(\sigma=-g\). The statistical error are too small to be seen. The black line shows the analytic expression for \(E_{0}\) when \(\sigma=-g/2\) for comparison.
pulsion allows the bosons to retain a tendency of remaining close to one another. On the other hand, we again find that if we instead turn on an attractive interaction between the particles, the local density as well as the pair-correlation peak sharply at the center of the harmonic trap and zero separation respectively. This structure is again characteristic to that of a bright soliton.

## V Conclusions

Neural networks and machine learning techniques provide powerful tools for the Variational Monte Carlo approach to solving many-body quantum systems. We used a neural network to parametrize a variational wavefunction ansatz for the ground state of systems of indistinguishable bosons in one dimension. Then using the ansatz, we were able to successfully numerically compute the ground state energies and wavefunctions of several quantum systems of various qualitative phases, ranging from two to fifty particles. The advantages of using a neural network ansatz include being able to search the most general space of possible ground state wavefunctions without needing to limit to specific functional forms. Neural networks also have convenient scalability; it is straightforward both to increase the number of variational parameters in the ansatz as well as increase the number of particles in the system. The application of such ansatze also does not need to be restricted to solvable systems and can be used to explore properties of non-integrable systems.

Many important many-body systems, including condensates, atoms, and nuclei involve indistinguishable particles. Wavefunctions of such systems must satisfy Fermi or Bose exchange symmetry. We advocated and demonstrated that using symmetric functions of the original particle coordinates as inputs to the neural network aptly enforces Bose exchange symmetry in one dimension. The symmetric inputs are chosen to be bijective to the original set of coordinates and so we do not limit the generality of the wavefunction space spanned by the neural network. In addition, we propose a method to deal with delta-function interactions without the need to regularize them.

The use of symmetric inputs to enforce exchange symmetries of the system can be generalized to use to study systems in a higher number of spatial dimensions in a straightforward way. To enforce Bose symmetry in three dimensions, one needs to find a bijection between the set of vector coordinates \(\{\vec{r}_{i}\},i=1,...,N\) and symmetric functions \(\{s_{j}(r_{1}^{x},r_{1}^{y},r_{1}^{z},...,r_{N}^{x},r_{N}^{y},r_{N}^{z})\},j= 1,...,M\) where \(M\) may be greater than \(3N\). The \(\{s_{j}\}\) must be symmetric in exchanges \(\vec{r}_{n}\leftrightarrow\vec{r}_{m}\) but not in exchanges of individual coordinate directions (i.e \(r_{n}^{x}\leftrightarrow r_{m}^{x}\)). Given the set of \(\{s_{j}\}\), one should be able to reconstruct the set of \(\{r_{i}^{x},r_{i}^{y},r_{i}^{z}\}\) including the fact that particular sets of three coordinates specify the position of individual particles. To write a Fermi anti-symmetric function in three dimensions, one theoretically simply needs to take a single Slater determinant of functions \(\phi_{i}(\vec{r}_{j};\{\vec{r}_{j}/\})\), where \(\phi_{i}\) is _Bose_ symmetric under exchange of any \(\vec{r}_{n}\rightarrow\vec{r}_{m},m,n\neq j\)[35]. Thus if one finds a bijection to a set of symmetric functions necessary for enforcing Bose symmetry in three dimensions, one can use the same bijection to enforce Fermi symmetry. Work to find such a mapping is being pursued.

###### Acknowledgements.

 We thank Andrei Alexandru and Scott Lawrence for helpful discussions and comments. This work was supported in part by the U.S. Department of Energy, Office of Nuclear Physics under Award Numbers DE-SC0021143, DE-FG02-93ER40762.

## Appendix A Importance Sampling for Delta-function Observables

Consider the Lieb-Liniger Hamiltonian [1; 2] for a system of \(N\) indistinguishable bosons

\[\hat{H}=\sum_{i=1}^{N}-\frac{1}{2}\frac{\partial^{2}}{\partial x_{i}^{2}}+ \sum_{i<j}^{N}\,g\delta(x_{i}-x_{j}) \tag{10}\]

and a positive and real wavefunction

\[\psi(X)=e^{-A(x_{1},...,x_{N})} \tag{11}\]

We would like to compute the energy associated with the wavefunction.

\[\langle E\rangle=\frac{\int dx_{1}..dx_{N}\,\psi\hat{H}\psi}{\int dx_{1}...dx _{N}\,\psi^{2}}=K+V \tag{12}\]

where

\[K=\frac{\int dx_{1}...dx_{N}\,e^{-2A(X)}\sum_{i}\,\frac{1}{2}\left(\frac{ \partial^{2}A}{\partial x_{i}^{2}}-\left(\frac{\partial A}{\partial x_{i}} \right)^{2}\right)}{\int dx_{1}...dx_{N}\,e^{-2A(X)}}\]

\[V=\frac{\int dx_{1}...dx_{N}\,e^{-2A(X)}\sum_{i<j}g\delta(x_{i}-x_{j})}{\int dx _{1}...dx_{N}\,e^{-2A(X)}} \tag{13}\]
The kinetic term \(K\) can be computed easily using importance sampling

\[K=\left\langle\,\sum_{i=1}^{N}\,\frac{1}{2}\left(\frac{\partial^{2}A}{\partial x _{i}^{2}}-\left(\frac{\partial A}{\partial x_{i}}\right)^{2}\right)\,\right\rangle _{\psi^{2}} \tag{10}\]

where \(\langle...\rangle_{\psi^{2}}\) denotes an average over samples drawn from the distribution

\[\Pi_{\psi^{2}}(x_{1},...x_{N})=\frac{e^{-2A(X)}}{\int dx_{1}...dx_{N}\,e^{-2A( X)}} \tag{11}\]

Using the same importance sampling technique to compute the potential term \(V\) leads to a zero-overlap problem as described in Sec. III.2. Instead we restructure the integral involved in \(V\) to become amenable to importance sampling using the following manipulations.

We first explicitly integrate out the \(x_{2}\) coordinate using the delta functions. The choice of \(x_{2}\) is arbitrary due to the Bose exchange symmetry of the wavefunction.

\[\begin{split} V&=g\frac{N(N-1)}{2}I\\ I&=\frac{\int dx_{1}dx_{3}...dx_{N}\,e^{-2A(\tilde{ X})}}{\int dx_{1}dx_{2}...dx_{N}\,e^{-2A(X)}}\end{split} \tag{12}\]

where \(\tilde{X}\) denotes the set of positions where \(x_{2}\) is replaced by \(x_{1}\) (i.e \(\{x_{1},x_{1},x_{3},...,x_{N}\}\)). We use the Bose symmetry of the wavefunction to eliminate the sum over pairs of particles. Next, we rewrite the integral \(I\) as

\[I=\frac{\int dx_{1}...dx_{N}\,\frac{\exp(-2A(\tilde{X}))}{\exp(-2A(X))} \mathcal{D}(x_{2})\exp(-2A(X))}{\int dx_{1}...dx_{N}\,\exp(-2A(X))} \tag{13}\]

where \(\mathcal{D}(x_{2})\) is any function such that

\[\int_{-\infty}^{\infty}dx_{2}\,\mathcal{D}(x_{2})=1 \tag{14}\]

In practice, we choose a Gaussian form

\[\mathcal{D}(x_{2})=\frac{1}{\alpha\sqrt{\pi}}e^{-x_{2}^{2}/\alpha^{2}} \tag{15}\]

The parameter \(\alpha\), which is different for each different \(\psi\), is chosen arbitrarily such that \(\mathcal{D}(y_{*})\approx 10^{-5}\) where \(y_{*}\) is the maximum value of \(|x_{2}|\) over many position samples taken from \(\Pi_{\psi^{2}}(X)\). In this way, \(\mathcal{D}(x_{2})\) is neither too wide such that we emphasize unimportant regions of \(\psi^{2}\) not sampled by the Monte Carlo method nor too narrow such that we don't capture all of the important regions of \(\psi^{2}\).

Now, we can also compute \(V\) by drawing samples from \(\Pi_{\psi^{2}}(X)\), just as we did for \(K\).

\[V=g\frac{N(N-1)}{2}\Big{\langle}\,\frac{e^{-2A(x_{1},x_{1},x_{3},...,x_{N})}}{ e^{-2A(x_{1},x_{2},...,x_{N})}}\mathcal{D}(x_{2})\Big{\rangle}_{\psi^{2}} \tag{16}\]

We note that this method of avoiding the zero-overlap problem occasionally produces large error estimates and fails to converge when \(e^{-2A(X)}\) is small for many different sampled \(X\). However, the method is still reliable because this is precisely the region of \(e^{-2A(X)}\) that is sampled the least by the Metropolis algorithm.
