Multi-body wave function of ground and low-lying excited states using
  unornamented deep neural networks
We propose a method to calculate wave functions and energies not only of the
ground state but also of low-lying excited states using a deep neural network
and the unsupervised machine learning technique. For systems composed of
identical particles, a simple method to perform symmetrization for bosonic
systems and antisymmetrization for fermionic systems is also proposed.
Tomoya Naito, Hisashi Naito, Koji Hashimoto
A simple method for multi-body wave function of ground and low-lying excited states using deep neural network

###### Abstract

We propose a method to calculate wave functions and energies not only of the ground state but also of low-lying excited states using a deep neural network and the unsupervised machine learning technique. For systems composed of identical particles, a simple method to perform symmetrization for bosonic systems and antisymmetrization for fermionic systems is also proposed.

+
Footnote †: preprint: RIKEN-iTHEMS-Report-23

KUNS-2954

## I Introduction

Atoms, molecules, and solids are composed of many electrons and ions, and atomic nuclei are composed of many nucleons. In principle, once the Schrodinger equation of these systems is solved, most properties can be described. However, in practice, they are quantum many-fermion systems, which are difficult to solve directly. Hence, it has been one of the important issues to solve the Schrodinger equation for the quantum many-fermion efficiently and accurately; in fact, many numerical methods, for instance, the Faddeev calculation [1], several methods for few-body systems [2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12], the quantum Monte Carlo method (QMC) including the variational Monte Carlo and diffusion Monte Carlo (DMC) methods [13; 14; 15; 16; 17], the configuration interaction method [18; 19; 20], the coupled cluster method [21; 22; 23; 24], and the density functional theory (DFT) [25; 26; 27], have been proposed in recent decades. These methods can be classified into two categories: methods based on the diagonalization of the Hamiltonian and on the variational principle.

Among above, DFT and QMC are classified into the latter. The variational principle [28] guarantees that the ground-state energy \(E_{\rm gs}\) of a Hamiltonian \(H\) satisfies

\[E_{\rm gs}=\inf\frac{\langle\Psi|H|\Psi\rangle}{\langle\Psi|\Psi\rangle}, \tag{1}\]

where all the possible functions are considered in the infimum. The minimizer corresponds to the ground-state wave function. Indeed, it is scarcely possible to consider all the possible functions when minimizing the energy expectation value; hence, in practice, the calculation accuracy of a method based on the variational principle depends on the ansatz of a trial wave function. In other words, the size of the space of trial wave functions determines the calculation accuracy.

For instance, a trial wave function of DFT is a Slater determinant, which is the simplest antisymmetric trial wave function. Owing to the simpleness of the ansatz, the numerical cost is drastically reduced, while it is known that inter-particle correlation is missing [29].

In the QMC calculation, a Jastrow-type trial wave function [30] is often used. A Jastrow-type wave function \(|\Psi\rangle\) consists of a single- (or sometimes multi-) Slater determinant \(|\Phi_{0}\rangle\) and a correlation factor \(F\), \(|\Psi\rangle=F\,|\Phi_{0}\rangle\). With assuming \(F\) as a symmetric function, \(|\Psi\rangle\) satisfies antisymmetry. Owing to the introduction of the factor \(F\), inter-particle correlations are described better than a single Slater determinant. Nevertheless, most QMC calculations optimize mainly \(F\), while \(|\Phi_{0}\rangle\) is optimized only around an initial ansatz [31]. In addition, an ansatz is introduced even for \(F\); hence, accuracy also depends on the ansatz. Recently, based on the QMC calculation, a deep neural network (DNN) has been used for the ansatz of a trial wave function [32; 33; 34]. Since deep neural networks span much wider space, calculation accuracy is much improved.

Another problem of variational-principle-based methods is calculation of excitation spectra. Excitation spectra are important quantities of molecules and atomic nuclei, while the variational principle [Eq. (1)] obtains the ground state only. Hence, another technique is needed to calculate excited states based on the variational principle. Indeed, a method to calculate low-lying excited states using the DMC calculation was proposed [35] by considering the orthogonality of wave functions, while its application has been still limited [36], whereas calculation of excited states on top of the DFT ground state has been widely performed by using the random phase approximation or some other techniques [37; 38; 39; 40; 41], while there is still room to be improved. Recently, low-lying excited states
were also obtained in Ref. [42] using the combination of the DNN, QMC, and the orthogonal condition.

In this paper, we propose a new method to calculate energies and wave functions of the ground state and low-lying excited states based on the variational principle. The ground-state wave function is assumed to be a DNN, which, in principle, is able to represent any function [43; 44]. Using an essence of the machine learning technique--the minimization of the loss function--, they are directly optimized by using the machine learning technique. References [45; 46] also assumed the ground-state wave function as a deep neural network that were optimized by using the machine learning technique: Ref. [45] performed calculation of few-body bosonic systems and Ref. [46] performed calculation of the simplest realistic system--a deuteron. Although these papers are pioneering works of unsupervised machine learnings for quantum many-body problems, the fermion antisymmetrization was not considered and excited states were not studied, while in quantum many-body problems both ground and excited states of many-fermionic systems are interesting in general. Recently, Ref. [47] proposed a method to obtain the ground state of the many-body Schrodinger equation for fermionic systems by using a tensor neural network, while its implementation is involved and the antisymmetrization is, indeed, not perfectly guaranteed.

In this paper, on top of the method in Refs. [45; 46] a simple method of the antisymmetrization for many-fermion systems or the symmetrization for many-boson systems is introduced. Then, low-lying excited states are sequentially calculated by using the orthogonality conditions and the variational principles. In this method, there is no need to discover a DNN architecture to generate (anti)symmetric wave functions. The symmetrization is put at the level of loss functions. Furthermore, the symmetrization and the antisymmetrization are implemented in almost the same way and the wave function perfectly satisfies (anti)symmetry. Thanks to the simplicity of the implementation, the numerical cost is quite small. We show that our method works successfully for popular examples in bosonic and fermionic quantum mechanical systems, providing a fundamental basis of the DNN method for quantum mechanics.

This paper is organized as follows: Section II is devoted to calculation of the ground-state. The novel antisymmetrization is introduced. Section III is devoted to the calculation of low-lying excited states. All the calculations are performed in a MacBook Pro with the Apple M1 chip (MacBook Pro (13-inch, M1, 2020): MacBookPro17,1) and \(16\,\mathrm{GB}\) memory. Section IV gives a summary of this paper.

## II Ground-state calculation

In this section, the ground-state wave function and energy are calculated by using a DNN and the machine learning technique. Throughout the paper, a machine learning software named Tensorflow[48] is used.

### Network structure and machine learning technique

In general, a wave function of a \(d\)-dimensional \(N\)-particle system is a function of the spatial coordinates of all the particles \(\mathbf{r}_{j}=(r_{j1},r_{j2},\ldots,r_{jd})\) (\(j=1,\ 2,\ \ldots,\ N\)). Here, to sake simplicity, we neglect the spin and isospin dependence of wave functions and \(\mathbf{R}\) denotes \(\mathbf{R}=(\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{N})=(r_{11},r_{12},\ldots,r_{1d}, r_{21},r_{22},\ldots,r_{2d},\ldots,r_{N1},r_{N2},\ldots,r_{Nd})\).

In this work, the wave function is represented by a deep neural network with \(Nd\)-input units that corresponds to the spatial coordinate \(\mathbf{R}\) and one-output unit that corresponds to the value of the wave function \(\psi\left(\mathbf{R}\right)\). Between the input and output layers, there are hidden layers. Each unit is connected to all the units just one before or after layers. The schematic figure of the deep neural network is shown in Fig. 1. In this paper, the "softplus" function

\[\mathrm{softplus}\left(x\right)=\log\left(1+e^{x}\right) \tag{2}\]

is used for an activation function and the adam optimizer [49] is used for the optimization process.

As the normal procedure of the numerical calculation of the wave function, the spatial coordinate is discretized. Each point is treated as a batch of the machine learning. In other words, if the spatial coordinate of each direction is discretized with \(M\) meshes, the batch size is the

Figure 1: Schematic figure of the deep neural network representing a one-dimensional three-body system.
same as the number of meshes, \(M^{dN}\). The mini-batch technique is not used.

Once the spatial coordinates are discretized, the Hamiltonian

\[H=-\frac{\hbar^{2}}{2m}\sum_{j}\Delta_{j}+\sum_{j}V^{\mathrm{ext}}\left(\mathbf{r }_{j}\right)+\frac{1}{2}\sum_{j\neq k}V^{\mathrm{int}}\left(\mathbf{r}_{j}, \mathbf{r}_{k}\right) \tag{3}\]

can be written as a matrix, where \(m\) is the mass of the particles, \(V^{\mathrm{ext}}\) is the external potential, and \(V^{\mathrm{int}}\) is the inter-particle interaction. The matrices of the external potential and the interaction are diagonal and that of the kinetic energy is sparse. Hence, the expectation value of the Hamiltonian \(\left\langle H\right\rangle\) can be calculated by using the sparse-matrix technique. The ground-state wave function minimizes \(\left\langle H\right\rangle\); hence, \(\left\langle H\right\rangle\) is regarded as a loss function. Note that all the calculations are performed with double precision floating point numbers (float64). For simplicity, \(m=\hbar=1\) is assumed.

The procedure in the Tensorflow code is as follows:

1. Construct a model of the deep neural network;
2. Note that although a Tensorflow subroutine specification for the loss function technically requires two inputs--the training data (true_value) and the network output (predicts), the former is not referred in our training;
3. Fit the model (model.fit) where the initial value of predicts consists of positive random numbers;
4. The final wave function output_wf is given by using model.predict;
5. The ground-state energy is calculated using the wave function obtained by the last step.

The third step (model.fit) corresponds to determining the parameters inside the DNN; the fourth step (model.predict) corresponds to storing the wave function obtained in the previous step; the fifth step corresponds to calculating the ground-state energy using the wave function obtained in the third step. Note that predicts and the final wave function should be normalized whenever generated.

### One-dimensional one-particle systems

In this section, benchmark calculations of one-dimensional systems are shown. The dependence of the numbers of units and layers on calculation accuracy is also discussed. Since there exists only one particle, there is no interaction, \(V^{\mathrm{int}}\equiv 0\); and thus the Hamiltonian reads

\[H=-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V^{\mathrm{ext}}\left(x\right). \tag{4}\]

Since we focus only on bound states in this paper, it is enough to deal with the limited spatial region. In the calculation, \(\left|x\right|\leq x_{\mathrm{max}}\) is considered and the box is discretized within \(1024\) meshes, i.e., \(M=1024\). The Dirichlet boundary condition (\(\psi\left(\pm x_{\mathrm{max}}\right)=0\)) is used. For the second derivative, the three-point derivative is used for simplicity, while it can be straightforwardly improved for the accuracy [50]. Then, \(H\) is discretized as

\[H\simeq\tilde{H} =-\frac{1}{2h^{2}}\tilde{T}+\tilde{V}^{\mathrm{ext}}, \tag{5a}\] \[\tilde{T} =\begin{pmatrix}-2&1&0&\ldots&0&0&0\\ 1&-2&1&\ldots&0&0&0\\ 0&1&-2&\ldots&0&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\ 0&0&0&\ldots&-2&1&0\\ 0&0&0&\ldots&1&-2&1\\ 0&0&0&\ldots&0&1&-2\end{pmatrix},\] (5b) \[\tilde{V}^{\mathrm{ext}} =\begin{pmatrix}V_{1}^{\mathrm{ext}}&0&0&\ldots&0&0&0\\ 0&V_{2}^{\mathrm{ext}}&0&\ldots&0&0&0\\ 0&0&V_{3}^{\mathrm{ext}}&\ldots&0&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\ 0&0&0&\ldots&V_{M-3}^{\mathrm{ext}}&0&0\\ 0&0&0&\ldots&0&V_{M-2}^{\mathrm{ext}}&0\\ 0&0&0&\ldots&0&0&V_{M-1}^{\mathrm{ext}}\end{pmatrix} \tag{5c}\]

and the wave function is also discretized as a \((M-1)\)-dimensional vector

\[\psi\simeq\tilde{\psi}=\begin{pmatrix}\psi_{1}\\ \psi_{2}\\ \psi_{3}\\ \vdots\\ \psi_{M-3}\\ \psi_{M-2}\\ \psi_{M-1}\end{pmatrix}, \tag{6}\]

where \(\tilde{\psi}\) is assumed to be normalized, i.e., \(h\sqrt{\sum_{j}\tilde{\psi}_{j}^{2}}=1\), \(V_{j}^{\mathrm{ext}}=V^{\mathrm{ext}}\left(x_{j}\right)\), \(\psi_{j}=\psi\left(x_{j}\right)\), \(x_{j}=-x_{\mathrm{max}}+hj\), and \(h\) denotes the mesh size \(h=2x_{\mathrm{max}}/M\)[51]. This \(\tilde{\psi}\) is used for predicts and output_wf. Here, a tilde denotes discretized form. Then, \(\left\langle H\right\rangle\) can be calculated as

\[\left\langle H\right\rangle\simeq\left\langle\tilde{H}\right\rangle=\tilde{ \psi}^{\mathrm{T}}\tilde{H}\tilde{\psi}. \tag{7}\]

#### ii.2.1 Harmonic oscillator

First of all, the harmonic oscillator potential

\[V^{\mathrm{ext}}\left(x\right)=\frac{1}{2}\omega^{2}x^{2} \tag{8}\]

is tested. The ground-state wave function \(\psi_{\mathrm{gs}}\) and energy \(E_{\mathrm{gs}}\) are, respectively, known exactly as [28]

\[\psi_{\mathrm{gs}}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\exp\left(-\frac{\omega x^{ 2}}{2}\right), \tag{9a}\] \[E_{\mathrm{gs}} =\frac{\omega}{2}. \tag{9b}\]
In this calculation, \(x_{\text{max}}=5\) is used.

Table 1 shows the summary of calculations. In general, all the calculations give almost the correct energy [Eq. (9b)]. On the one hand, total optimization costs similar amount of time in all the calculation. On the other hand, different setup requires different number of epochs and time per epoch for optimizing the DNN. Small DNN tends to take shorter time for each epoch, while it requires longer epochs. It seems that 32 units per layer is too large, so it requires longer epoch and longer estimation time per epoch. It should be noted that the number of epochs differs in each time since the initial condition of the fitting procedure is generated by the random numbers. In addition, if one uses a different value of the learning rate, the number of epochs can be different.

Figure 2 shows relative errors of the loss function, \(\left\langle H\right\rangle\), to the exact ground-state energy \(E_{\text{gs}}\) as functions of the number of epochs. It can be seen that, although the loss function achieved the relative error of \(1.0\times 10^{-8}\), the final accuracy becomes about \(1.0\times 10^{-4}\). This may be due to the precision of the Tensorflow code.

Figure 3 shows calculated wave functions. The red thick lines correspond to the exact solution given in Eq. (9a), while thin lines corresponds to the results given in this work, where different colors correspond to different numbers of units and layers. The relative errors of the DNN wave function, \(\psi^{\text{DNN}}\), to the exact one, \(\psi^{\text{exact}}\),

\[\delta\psi\left(x\right)=\frac{\left|\psi^{\text{DNN}}\left(x\right)-\psi^{ \text{exact}}\left(x\right)\right|}{\psi^{\text{exact}}\left(x\right)} \tag{10}\]

are shown in Fig. 4. It can be seen that the DNN calculation, basically, reproduces the exact solution in our interest within the accuracy of \(10^{-4}\) or more. This deviation can be reduced if we use more tight convergence criterion [52]. In the tail region, the deviation \(\delta\psi\left(x\right)\) diverges, while this is because the denominator of Eq. (10), \(\psi^{\text{exact}}\left(x\right)\), reaches to zero. The deviation looks larger if \(\omega\) is smaller, which is related to the cutoff parameter for the spatial mesh \(x_{\text{min}}\). The exact value of \(\psi_{\text{gs}}\left(x_{\text{min}}\right)\) is \(2.8\times 10^{-6}\) for \(\omega=1\), while it is \(8.1\times 10^{-28}\) for \(\omega=5\) and it is much smaller for \(\omega=10\), while in the numerical calculation, they are approximated to zero. The value \(2.8\times 10^{-6}\) may be too large to assume to be zero.

It should be noted that rather small DNN is enough to reproduce the solution of the wave function. Owing to the simplicity, it is easy to analyse the weights and biases of the DNN. For instance, the DNN wave function for the single layer with four units includes only 13 parameters; the ground-state DNN wave function for \(\omega=1\) can be written as

\[\psi_{\text{gs}}\left(x\right) =\frac{1}{3.7451}\,\text{softplus}\left(a_{\text{gs}}\left(x \right)\right), \tag{11a}\] \[a_{\text{gs}}\left(x\right) =2.4069a_{1}\left(x\right)-1.8344a_{2}\left(x\right)-1.9778a_{3} \left(x\right)+2.3484a_{4}\left(x\right)-4.8998,\] (11b) \[a_{1}\left(x\right) =\text{softplus}\left(0.35953x+3.9226\right),\] (11c) \[a_{2}\left(x\right) =\text{softplus}\left(2.5821x+0.033213\right),\] (11d) \[a_{3}\left(x\right) =\text{softplus}\left(-0.65170x+2.9574\right),\] (11e) \[a_{4}\left(x\right) =\text{softplus}\left(0.15421x+2.2016\right), \tag{11f}\]

where the first coefficient of Eq. (11a) (\(1/3.7451\)) is not obtained by the DNN but instead by the normalization. Hence, smaller DNN is better not only due to the calculation cost but also for analysis of the structure of DNN.

Let us provide our interpretation of the obtained wave function [Eq. (11a)]. The rectified linear function (ReLU)

\[\text{ReLU}\left(x\right)=\begin{cases}0&\left(x<0\right),\\ x&\left(x\geq 0\right)\end{cases} \tag{12}\]

is a widely-used activation function, and the softplus function can be regarded as a smoothed version of the ReLU. Here, for interpreting Eq. (11a), we shall just replace the softplus function with the ReLU. The wave function obtained by the DNN [Eq. (11a)] and the obtained function before the output layer [Eq. (11b)] are shown in Fig. 5, where the normalization factor (\(1/3.7451\)) of \(\psi_{\text{gs}}\) is ignored. Equations (11a) and (11b) where the softplus function is replaced to the ReLU function are also plotted as \(\psi_{\text{gs}}^{\text{ReLU}}\) and \(a_{\text{gs}}^{\text{ReLU}}\), respectively. The DNN with the ReLU function can be understood as an approximation with a piecewise linear function. As shown in Fig. 5, the ground-state wave function in DNN is approximated by the following function:

\[\psi_{\text{gs}}\approx\begin{cases}-ax+b&\left(0\leq x\leq b/a\right),\\ ax+b&\left(-b/a\leq x\leq 0\right),\\ 0&\left(\text{otherwise}\right),\end{cases} \tag{13}\]

where \(a\) and \(b\) are positive numbers. The ReLU function at the output layer guarantees to make the wave function vanish for \(x<-b/a\) and \(x>b/a\), and thus, \(a_{\text{gs}}\) should be \(\pm ax+b\). This function can be represented by just two ReLU functions. Hence, even two units in hidden layer are enough to describe the brief structure of ground-state wave function, and with increase the number of units, the
ground-state wave function is reproduced easily. Since the ReLU function is not differentiable at \(x=0\), the ReLU wave function is not differentiable. Hence, the softplus is better to describe a differentiable function, while the ReLU function can also describe a differentiable function approximately if the number of units is large enough. In case of \(N\)-bodies systems, the similar function to Eq. (11a) can be represented by just \(2^{N}\) ReLU functions.

#### iii.2.2 Square wall potential

Next, the square wall potential

\[V^{\mathrm{ext}}\left(x\right)=\begin{cases}-V_{0}&(|x|<x_{0}),\\ 0&(\text{otherwise})\end{cases} \tag{14}\]

is tested (\(V_{0}>0\)). The analytical forms of the ground-state wave function \(\psi_{\mathrm{gs}}\) and energy \(E_{\mathrm{gs}}\) are unknown; thus, our values of the energy will be compared with the numerical calculation obtained by the orthodox method of Hamiltonian diagonalization. In this calculation, \(x_{\mathrm{max}}=20\) and \(x_{0}=1\) is used.

Table 2 shows the summary of calculations. In general, all the calculations give almost the correct energy. The calculation time per epoch and the number of epochs with respect to the number of layers and units is slightly longer than the case of the harmonic oscillator.

Figure 6 shows relative errors of the loss function, \(\langle H\rangle\), to the exact ground-state energy \(E_{\mathrm{gs}}\) as functions of the number of epochs. It can be seen that, although the loss function achieved the relative error of \(1.0\times 10^{-7}\), the final accuracy becomes about \(1.0\times 10^{-2}\). Note that in this calculation, the convergence criteria is needed to be set looser than the other case; otherwise, it could not reach convergence. This may be related to the shape of the potential: Asymptotic region of the square wall potential is zero, while the harmonic oscillator potential increases rapidly. It will be shown later that the double-wall potential, which is close to the latter situation, reaches convergence with the tight criterion.

Figure 7 shows calculated wave functions. The red thick lines correspond to the exact solution given by the exact diagonalization, where the same mesh size matrix form are used for comparison; thin lines corresponds to the results given in this work. It can be seen that the DNN calculation, basically, reproduces the solutions given by the exact diagonalization.

### One-dimensional many-particle systems

When one considers systems composed of many identical particles, the symmetrization for bosonic systems or the antisymmetrization for fermionic systems of the wave function must be considered. The ground state of the bosonic system is identical to that of the different particles; hence it has no extra difficulty as was done in Ref. [45], while the antisymmetrization is rather difficult. In this section, a simple method of (anti)symmetrization in the DNN wave function is provided, in which the symmetrization and the antisymmetrization can be performed with equal footing.

#### iii.3.1 Hamiltonian matrix

As was done in the last section, the discretized Hamiltonian \(\tilde{H}\) should be represented in a matrix form and the discretized wave function \(\tilde{\psi}\) should be represented in a vector form. Here, one-dimensional two-body systems are considered as an example, and their coordinates are denoted by \(x\) and \(y\). Each direction is discretized with \(M\) meshes, i.e., in total \(M\times M\) meshes. Then, the discretized wave function \(\tilde{\psi}\) is

\[\tilde{\psi}=\begin{pmatrix}&\psi_{11}\\ &\psi_{12}\\ &\vdots\\ &\psi_{1(M-1)}\\ &\psi_{21}\\ &\psi_{22}\\ &\vdots\\ &\psi_{2(M-1)}\\ &\vdots\\ &\psi_{(M-1)1}\\ &\psi_{(M-1)2}\\ &\vdots\\ &\psi_{(M-1)(M-1)}\end{pmatrix}, \tag{15}\]

where \(\psi_{jk}=\psi\left(x_{j},y_{k}\right)\), \(x_{j}=-x_{\mathrm{max}}+hj\), \(y_{k}=-y_{\mathrm{max}}+hk\), and \(x_{\mathrm{max}}=y_{\mathrm{max}}\). Accordingly, the discretized Hamiltonian \(\tilde{H}\) reads

\[\tilde{H}=-\frac{1}{2h^{2}}\tilde{T}_{1}-\frac{1}{2h^{2}}\tilde{T}_{2}+\tilde{ V}_{\mathrm{ext}}^{1}+\tilde{V}_{\mathrm{ext}}^{2}+\tilde{V}_{\mathrm{int}}, \tag{16}\]

where \(\tilde{T}_{1}\) and \(\tilde{T}_{2}\) are the kinetic energy matrices

\[\tilde{T}_{1} =T\otimes I_{2}, \tag{17a}\] \[\tilde{T}_{2} =I_{2}\otimes T, \tag{17b}\]

\(\tilde{V}_{\mathrm{ext}}^{1}\) and \(\tilde{V}_{\mathrm{ext}}^{2}\) are the external potential matrices

\[\tilde{V}_{1}^{\mathrm{ext}} =V^{\mathrm{ext}}\otimes I_{2}, \tag{18a}\] \[\tilde{V}_{2}^{\mathrm{ext}} =I_{2}\otimes V^{\mathrm{ext}}, \tag{18b}\]

and \(\tilde{V}_{\mathrm{int}}\) are the interaction matrix whose matrix elements are
\[\left(\tilde{V}_{\text{int}}\right)_{i+j(M-1),k+l(M-1)}=\begin{cases}\frac{1}{2} \left[V^{\text{int}}\left(x_{i},y_{j}\right)+V^{\text{int}}\left(y_{j},x_{i} \right)\right]=V^{\text{int}}\left(x_{i},y_{j}\right)&\text{(for $i=k$, $j=l$)},\\ 0&\text{(otherwise)},\end{cases} \tag{19}\]

and \(\otimes\) is the Kronecker product. For instance, the matrix elements of Eqs. (17) reads

\[\left(\tilde{T}_{1}\right)_{i+j(M-1),k+l(M-1)}=\begin{cases}-2&\text{(for $i=k$, $j=l$)},\\ 1&\text{(for $i=k\pm 1$, $j=l$)},\\ 0&\text{(otherwise)},\end{cases} \tag{20a}\] \[\left(\tilde{T}_{2}\right)_{i+j(M-1),k+l(M-1)}=\begin{cases}-2&\text{(for $i =k$, $j=l$)},\\ 1&\text{(for $i=k$, $j=l\pm 1$)},\\ 0&\text{(otherwise)}.\end{cases} \tag{20b}\]

\begin{table}
\begin{tabular}{l l l l l l l l} \(\omega\) & \multicolumn{2}{c}{\# of Unit} & \multicolumn{2}{c}{Energy} & \multicolumn{2}{c}{\# of Epochs} & \multicolumn{1}{c}{Time per Epoch} \\ \cline{2-7}  & \multicolumn{1}{c}{1st Layer} & \multicolumn{1}{c}{2nd Layer} & \multicolumn{1}{c}{Kinetic} & \multicolumn{1}{c}{Kinetic} & \multicolumn{1}{c}{Potential} & \multicolumn{1}{c}{Total} & \multicolumn{1}{c}{(\#s)} \\ \hline

[MISSING_PAGE_POST]

 & +2.499976 & +2.499745 & +4.999721 & 19833 & 996.732 \\ \end{tabular}
\end{table}
Table 1: Calculation summary of a one-body problem under the harmonic oscillator potential. Row with “—” in the column “# of unit for 2nd layer” corresponds to calculation performed only with one layer.
#### ii.1.2 Symmetrization and antisymmetrization

The discretized wave function \(\tilde{\psi}\) should be symmetrized or antisymmetrized. In general, for the arbitrary function \(f\left(x,y\right)\), \(f\left(x,y\right)+f\left(y,x\right)\) is a symmetrized function and \(f\left(x,y\right)-f\left(y,x\right)\) is an antisymmetrized function.

In order to perform (anti)symmetrization in TensorFlow, instead of the simple \(\tilde{\psi}\),

\[\tilde{\psi}_{\pm}=\left(\begin{array}{c}\psi_{11}\\ \psi_{12}\\ \vdots\\ \psi_{1(M-1)}\\ \psi_{21}\\ \psi_{22}\\ \vdots\\ \psi_{2(M-1)}\\ \vdots\\ \psi_{(M-1)(M-1)}\end{array}\right)\pm\left(\begin{array}{c}\psi_{11}\\ \psi_{21}\\ \vdots\\ \psi_{12}\\ \psi_{22}\\ \vdots\\ \psi_{(M-1)2}\\ \vdots\\ \psi_{1(M-1)}\\ \psi_{2(M-1)}\\ \vdots\\ \psi_{(M-1)(M-1)}\end{array}\right) \tag{22}\]

is assumed to be a trial wave function. In the Tensorflow code, instead of the original predicts and output_wf, \(\tilde{\psi}_{\pm}\) is used in the second (the calculation process of the loss function) and fifth (calculate the ground-state energy) steps in Sec. II.1. Note that this process can be easily done by using the following commands:

1. predicts_transpose = tf.reshape(predicts,

Figure 3: Wave function under the harmonic oscillator potential. The red thick line corresponds to the exact solution [Eq. (9a)], while thin lines correspond to results of DNN calculation. Different thin line corresponds to different number of units. We observe that all the simulated results overlap with the exact solution.

Figure 2: Relative error of \(\left\langle H\right\rangle\) to the exact ground-state energy \(E_{\text{gs}}\) for the harmonic oscillator potential as functions of the number of epochs.
[m, m]),
2. predicts_transpose = tf.transpose(predicts_transpose),
3. predicts_transpose = tf.reshape(predicts_transpose, [m**2, 1]),
4. predicts = tf.add(predicts, predicts_transpose) for bosonic systems or predicts = tf.subtract(predicts, predicts_transpose) for fermionic systems,
5. the final predicts is used to evaluate the loss function.

Here, m corresponds to the number of meshes \(M\). Note that this method can be straightforwardly extended to multi-body systems.

#### iii.2.3 Two-body systems

Two-body systems under the harmonic oscillator potential [Eq. (8)] is tested. If there is no interaction \(V^{\mathrm{int}}\equiv 0\), the ground-state wave function \(\psi_{\mathrm{gs}}\) and energy \(E_{\mathrm{gs}}\) are known exactly. If one considers bosonic systems, they read

\[\psi_{\mathrm{gs}}\left(x,y\right) =\sqrt{\frac{\omega}{\pi}}\exp\left[-\frac{\omega\left(x^{2}+y^{2 }\right)}{2}\right], \tag{23a}\] \[E_{\mathrm{gs}} =\omega,\] (23b) and if one considers the fermionic systems, they read \[\psi_{\mathrm{gs}}\left(x,y\right) =\frac{\omega}{\sqrt{\pi}}\left(x-y\right)\exp\left[-\frac{ \omega\left(x^{2}+y^{2}\right)}{2}\right],\] (24a) \[E_{\mathrm{gs}} =2\omega. \tag{24b}\]

Figures 8 and 9, respectively, show wave functions for bosonic and fermionic systems obtained in this work. The total energies and the calculation time are shown in Table 3. For comparison, the exact wave functions [Eq. (23a) or (24a)] are also shown. Here, \(x_{\mathrm{max}}=y_{\mathrm{max}}=5\) and \(M_{x}=M_{y}=256\) are used for the spatial mesh and two layers each of which contains 32 units are used for the DNN. For the interaction, the Gaussian-type interaction

\[V^{\mathrm{int}}\left(x,y\right)=\lambda\exp\left(-\left|x-y\right|\right) \tag{25}\]

is used, where \(\lambda\) is the strength of the interaction. The DNN results with \(\lambda=0\) show a good agreement with the exact results, demonstrating that the DNN technique works well.

Behaviour of wave functions for nonzero \(\lambda\) is consistent qualitatively with our expectation: if \(\lambda\) is negative, i.e., the interaction is attractive, the wave function tends to collapse, and if \(\lambda\) is positive, i.e., the interaction is repulsive, the wave function tends to be broad. Time cost per epoch is almost universal among all the calculation, while more epochs are required to reach convergence for

Figure 4: Relative error of DNN wave function to exact one. Different thin line corresponds to different number of units. In the region of large \(\left|x\right|\), the deviation \(\delta\psi\left(x\right)\) diverges, because the denominator of Eq. (10), \(\psi^{\mathrm{exact}}\left(x\right)\), reaches to zero. Hence, the figures plot only \(\left|\delta\varphi\left(x\right)\right|<10^{-1}\).

Figure 5: Wave function obtained by the DNN [Eq. (11a)] and the obtained function before the output layer [Eq. (11b)], where the normalization is ignored. Equations (11a) and (11b) where the softplus function is replaced with the ReLU function are also plotted as \(\psi_{\mathrm{gs}}^{\mathrm{ReLU}}\) and \(a_{\mathrm{gs}}^{\mathrm{ReLU}}\), respectively.
fermionic systems than for bosonic systems. This may be because all the values are positive for the initial condition, while there are negative values for fermionic ground-state wave functions. More epochs are required for the repulsive interaction (\(\lambda>0\)) than for the attractive interaction (\(\lambda<0\)). This may be because the topology of the wave function is more complicated and extended in the repulsive case than the attractive case.

Finally, we point out a strange behaviour of the obtained wave function of the two-body system of \(\omega=1\) without the interaction. Here, for simplicity, the two-layer DNN in which each unit is composed of 4 units is used. In the case of two layers, the function obtained by

\begin{table}
\begin{tabular}{c c c c c c c c} \(V_{0}\) & \multicolumn{2}{c}{\# of Unit} & \multicolumn{2}{c}{Energy} & \multicolumn{2}{c}{\# of Epochs} & Time per Epoch \\ \cline{2-7}  & \multicolumn{2}{c}{1st Layer} & \multicolumn{1}{c}{2nd Layer} & Kinetic & Potential & Total & (\(\mu\)s) \\ \hline

[MISSING_PAGE_POST]

 & \(+0.659308\) & \(-9.847178\)
optimized weights of DNN is

\[u_{\text{gs}}\left(x,y\right) =A\,\text{softplus}\left(\sum_{j=1}^{N_{\text{unit}}}w_{2j}u_{2j} \left(x,y\right)+b_{2}\right), \tag{26a}\] \[u_{2j}\left(x,y\right) =\text{softplus}\left(\sum_{k=1}^{N_{\text{unit}}}w_{1jk}u_{1k} \left(x,y\right)+b_{1j}\right),\] (26b) \[u_{1k}\left(x,y\right) =\text{softplus}\left(w_{0k0}x+w_{0k1}y+b_{0k}\right), \tag{26c}\]

where \(A\) is the normalization constant, \(N_{\text{unit}}\) is the number of units of each layer, \(k\) is a weight, and \(b\) is a bias. The left column of Fig. 10 shows \(u_{\text{gs}}\left(x,y\right)\). The upper and lower rows, respectively, correspond to the results with minimizing the bosonic or fermionic energy expectation value. It is shown that the obtained function \(u_{\text{gs}}\), which is referred to as the raw _wave function_, is not symmetric nor antisymmetric. After the symmetrization for bosonic systems \(\psi_{\text{boson}}\left(x,y\right)=\left[u_{\text{gs}}\left(x,y\right)+u_{ \text{gs}}\left(y,x\right)\right]/A_{\text{boson}}\) or the antisymmetrization for fermionic systems \(\psi_{\text{fermion}}\left(x,y\right)=\left[u_{\text{gs}}\left(x,y\right)-u_{ \text{gs}}\left(y,x\right)\right]/A_{\text{fermion}}\) is performed with the normalization constant \(A_{\text{boson}}\) or \(A_{\text{fermion}}\), \(\psi_{\text{boson}}\) or \(\psi_{\text{fermion}}\) can be regarded as the bosonic or fermionic ground-state wave function, respectively. The energy expectation value of the raw (\(u_{\text{gs}}\)), the symmetrized (\(\psi_{\text{boson}}\)), and the antisymmetrized (\(\psi_{\text{fermion}}\)) wave functions are shown in Table 4. A surprising fact is that even if the raw _wave function_ is obtained by minimizing the bosonic (fermionic) expectation value, fermion (bosonic) energy expectation value is close to the cor

Figure 6: Relative error of \(\langle H\rangle\) to the exact ground-state energy \(E_{\text{gs}}\) for the square wall potential as functions of the number of epochs.
rect fermion (bosonic) energy eigenvalue, and _vice versa_. The (anti)symmetrization corresponds to the projection of the raw _wave function_ to the boson (fermion) subspace, while the remaining part is not supposed to be optimized well. This unexpected coincidence may be due to the smallness of the number of parameters in the DNN architecture, and deserves further study.

#### iv.2.4 Three-body systems

Three-body systems under the harmonic oscillator potential [Eq. (8)] is tested. For simplicity, we consider a system without any interaction \(V^{\text{int}}\equiv 0\). Then, the ground-state wave function \(\psi_{\text{gs}}\) and energy \(E_{\text{gs}}\) are known exactly as

\[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27a) \[E_{\text{gs}} =\frac{3}{2}\omega,\] (27b) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27c) \[E_{\text{gs}} =\frac{3}{2}\omega,\] (27d) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27e) \[E_{\text{gs}} =\frac{3}{2}\omega,\] (27f) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27g) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27h) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27h) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27i) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27h) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27i) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27h) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27i) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27i) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27j) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega \left(x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega \left(x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right) =\left(\frac{\omega}{\pi}\right)^{3/4}\exp\left[-\frac{\omega\left( x^{2}+y^{2}+z^{2}\right)}{2}\right],\] (27k) for bosonic systems and \[\psi_{\text{gs}}\left(x,y,z\right)
\[\psi_{\rm gs}\left(x,y,z\right)=\left(\frac{\omega}{\pi}\right)^{3/4}\sqrt{\frac{ \omega}{6}}\left[\left(x-y\right)\left(1-2\omega z^{2}\right)+\left(y-z\right) \left(1-2\omega x^{2}\right)+\left(z-x\right)\left(1-2\omega y^{2}\right) \right]\exp\left[-\frac{\omega\left(x^{2}+y^{2}+z^{2}\right)}{2}\right], \tag{28a}\] \[E_{\rm gs}=\frac{9}{2}\omega \tag{28b}\]

for fermionic systems.

Figures 11 and 12, respectively, show wave functions for bosonic and fermionic systems obtained by this work. The total energies and the calculation time are shown in Table 5. For comparison, the exact wave functions [Eq. (27a) or (28a)] are also shown. Here, \(x_{\rm max}=y_{\rm max}=z_{\rm max}=5\) and \(M=64\) are used for the spatial mesh and two layers each of which contains 32 units are used for the DNN. The interaction is not considered.

The DNN calculations

\begin{table}
\begin{tabular}{l c c c}  & Raw & Symmetrized & Antisymmetrized \\ \hline Boson & 1.13190 & 1.00012 & 2.05195 \\ Fermion & 1.13628 & 1.08235 & 2.00056 \\ \end{tabular}
\end{table}
Table 4: Energy expectation value of the raw (\(u_{\rm gs}\) in Eq. (26a)), the symmetrized, and the antisymmetrized wave functions. The rows named “Boson” and “Fermion”, respectively, correspond to the results obtained by minizing the bosonic or fermionic energy expectation values.

\begin{table}
\begin{tabular}{l c c c c} Particles & \(\omega\) & \(\lambda\) & Energy & \# of Epochs & Time per Epoch (ms) \\ \hline Boson & 1.0 & \(-1.00\) & \(-89.869381\) & 8490 & 23.583 \\ Boson & 1.0 & \(-0.25\) & \(-19.848949\) & 3595 & 23.618 \\ Boson & 1.0 & \(+0.00\) & \(+0.999927\) & 27242 & 23.424 \\ Boson & 1.0 & \(+0.25\) & \(+3.298725\) & 20040 & 23.529 \\ Boson & 1.0 & \(+1.00\) & \(+3.835173\) & 21763 & 23.677 \\ \hline Boson & 5.0 & \(-1.00\) & \(-87.554311\) & 10194 & 23.573 \\ Boson & 5.0 & \(-0.25\) & \(-17.203647\) & 10893 & 23.712 \\ Boson & 5.0 & \(+0.00\) & \(+4.997829\) & 24477 & 23.648 \\ Boson & 5.0 & \(+0.25\) & \(+21.149827\) & 20721 & 23.797 \\ Boson & 5.0 & \(+1.00\) & \(+31.804917\) & 21718 & 23.591 \\ \hline Boson & 10.0 & \(-1.00\) & \(-84.129658\) & 18635 & 23.566 \\ Boson & 10.0 & \(-0.25\) & \(-13.118213\) & 24380 & 23.601 \\ Boson & 10.0 & \(+0.00\) & \(+9.991009\) & 27489 & 23.692 \\ Boson & 10.0 & \(+0.25\) & \(+32.287424\) & 23810 & 23.816 \\ Boson & 10.0 & \(+1.00\) & \(+72.688350\) & 19591 & 23.601 \\ \hline Fermion & 1.0 & \(-1.00\) & \(-71.409493\) & 18794 & 23.928 \\ Fermion & 1.0 & \(-0.25\) & \(-11.369632\) & 17999 & 23.818 \\ Fermion & 1.0 & \(+0.00\) & \(+1.999931\) & 19215 & 23.843 \\ Fermion & 1.0 & \(+0.25\) & \(+3.298786\) & 25187 & 23.804 \\ Fermion & 1.0 & \(+1.00\) & \(+3.839178\) & 106163 & 24.136 \\ \hline Fermion & 5.0 & \(-1.00\) & \(-68.409494\) & 18558 & 23.731 \\ Fermion & 5.0 & \(-0.25\) & \(-7.207718\) & 19618 & 23.787 \\ Fermion & 5.0 & \(+0.00\) & \(+9.995902\) & 75208 & 23.841 \\ Fermion & 5.0 & \(+0.25\) & \(+21.385884\) & 51691 & 23.607 \\ Fermion & 5.0 & \(+1.00\) & \(+31.804915\) & 27885 & 23.799 \\ \hline Fermion & 10.0 & \(-1.00\) & \(-63.026667\) & 26603 & 23.894 \\ Fermion & 10.0 & \(-0.25\) & \(+0.302312\) & 21783 & 23.805 \\ Fermion & 10.0 & \(+0.00\) & \(+19.975414\) & 10135 & 23.741 \\ Fermion & 10.0 & \(+0.25\) & \(+38.060907\) & 79973 & 24.054 \\ Fermion & 10.0 & \(+1.00\) & \(+73.245097\) & 59289 & 23.972 \\ \end{tabular}
\end{table}
Table 3: Calculation summary of a two-body problem under the harmonic oscillator potential.
state energies. The DNN wave functions are consistent with the exact solution. The number of epochs for three-body systems is comparable with those for two-body systems, where the number of units and layers are identical for these two cases. In contrast, the time per epoch for the three-body systems are about four times of that for the two-body systems. This is related to the number of spatial meshes: \(256\times 256=65536\) meshes are used for the two-body systems and \(64\times 64\times 64=262144\) meshes are used for the three-body systems; thus, the number of meshes for the three-body systems are four times more than those for the two-body systems. Hence, it can be concluded that the time per epoch is almost proportional to the number of spatial meshes. This is reasonable since we use numerical methods for sparse matrices, in which the number of the non-zero matrix elements is \(O\left(M^{Nd}\right)\).

Figure 8: Two-body wave function under the harmonic oscillator potential for bosonic systems. The exact wave function without the interaction is shown in the left-most column.

Figure 9: Same as Fig. 8 but for fermionic systems.
## III Excited-state calculation

In this section, based on the variational principle, a method to calculate low-lying excited states sequentially is explained. Assume that wave functions of the ground state and \(n\) excited states, \(\ket{\psi_{0}}\), \(\ket{\psi_{1}}\), \(\dots\), \(\ket{\psi_{n}}\), are derived, where \(\ket{\psi_{0}}=\ket{\psi_{\text{gs}}}\). We consider a problem of finding the \((n+1)\)-th excited state \(\ket{\psi_{n+1}}\) which satisfies the orthnormal condition

\[\left\langle\psi_{j}|\psi_{n+1}\right\rangle=\delta_{j,n+1}, \tag{29}\]

by using a trial wave function \(\ket{\psi}\). The \((n+1)\)-th wave function can be obtained with minimizing the expectation value

\[\left\langle H\right\rangle=\frac{\left\langle\psi|H|\psi\right\rangle}{ \left\langle\psi|\psi\right\rangle}, \tag{30}\]

where \(\ket{\psi}\) is assumed to be orthogonal to \(\ket{\psi_{j}}\) (\(j=0\), \(1\), \(\dots\), \(n\)). This can be implemented in Tensorflow with assuming that

\[\ket{\psi}-\sum_{j=0}^{n}\left\langle\psi_{j}|\psi\right\rangle\ket{\psi_{j}} \tag{31}\]

is a trial wave function, instead of the simple \(\ket{\psi}\). For one-body problem, \(x_{\text{max}}=5\) and \(M=1024\) are used for the spatial mesh and the single-layer DNN with eight units is adopted.

\begin{table}
\begin{tabular}{l c c c} Particles & Energy & \# of Epochs & Time per Epoch (ms) \\ \hline Boson & +1.497880 & 20183 & 101.216 \\ Fermion & +4.486830 & 22770 & 98.356 \\ \end{tabular}
\end{table}
Table 5: Calculation summary of a three-body problem under the harmonic oscillator potential. Calculation is performed with \(\omega=1.0\).

Figure 11: (Left) Three-body wave function under the harmonic oscillator potential without inter-particle interaction for bosonic systems. (Right) Slice of the three-body wave function at the plane \(x+y+z=0\).

Figure 10: DNN wave function of the raw (\(u_{\text{gs}}\) in Eq. (26a)), the symmetrized, and the antisymmetrized wave functions. The rows named “Boson” and “Fermion”, respectively, correspond to the results obtained by minizing the bosonic or fermionic energy expectation values.
### Harmonic oscillator potential

One-body one-dimensional harmonic oscillators are taken as examples. The exact wave functions for several low-lying excited states are [28]

\[\psi_{0}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\exp\left(-\frac{\omega x^{2 }}{2}\right), \tag{32a}\] \[\psi_{1}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\sqrt{2\omega}x\exp\left(- \frac{\omega x^{2}}{2}\right),\] (32b) \[\psi_{2}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\frac{2\omega x^{2}-1}{ \sqrt{2}}\exp\left(-\frac{\omega x^{2}}{2}\right),\] (32c) \[\psi_{3}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\sqrt{\frac{\omega}{3}} \left(2\omega x^{2}-3\right)x\exp\left(-\frac{\omega x^{2}}{2}\right),\] (32d) \[\psi_{4}\left(x\right) =\left(\frac{\omega}{\pi}\right)^{1/4}\frac{4\omega^{2}x^{4}-12 \omega x^{2}+3}{2\sqrt{6}}\exp\left(-\frac{\omega x^{2}}{2}\right), \tag{32e}\]

where \(\psi_{n}\) is the \(n\)-th excited state, and the energies are

\[E_{n}=\left(n+\frac{1}{2}\right)\omega. \tag{33}\]

Figure 13 shows the wave functions of the ground state and first, second, third, and fourth excited states. Table 6 shows the summary of calculations. Basically, not only the ground-state but also low-lying excited-states wave functions and energies are successfully calculated.

One can find that the DNN solutions are consistent with the exact solutions. Thus, it can be concluded that the method to calculate low-lying excited states proposed here works well.

The number of epochs are almost universal for all states calculated here. In contrast, the time per epoch for a higher excited state is slightly longer since calculation for orthogonal condition [Eq. (31)] is needed to be performed, while it takes just a few \(\mu\)s.

Let us explain why our simple DNN can describe even the excited states correctly. For simplicity, the single-layer DNN with the four unit is used. The obtained function for the ground-state (\(u_{\text{gs}}\)), the first (\(u_{\text{1st}}\)), and the second (\(u_{\text{2nd}}\)) excited states are shown in Fig. 14. The obtained function for the \(n\)-th excited state is

\[u_{n\text{-th}}\left(x\right)=\sum_{j=0}^{n}a_{j}\psi_{n\text{-th}}\left(x\right) \tag{34}\]

with \(\sum_{j=0}^{n}\left|a_{j}\right|^{2}=1\), where \(\psi_{n\text{-th}}\) is the \(n\)-th excited-state wave function. In the case of the first and second excited states,

\[u_{\text{1st}}\left(x\right) =0.981758\psi_{\text{0th}}\left(x\right)+0.190136\psi_{\text{1st }}\left(x\right), \tag{35a}\] \[u_{\text{2nd}}\left(x\right) =0.969707\psi_{\text{0th}}\left(x\right)-0.0445505\psi_{\text{1st }}\left(x\right)+0.240173\psi_{\text{2nd}}\left(x\right) \tag{35b}\]

### Double-well potential

In order to see the effect of degeneracy, we also test the double-well potential

\[V^{\text{ext}}\left(x\right)=\left(x^{2}-\alpha^{2}\right)^{2}. \tag{36}\]

If the central barrier is low enough, i.e., \(\alpha\) is small enough, each state is not degenerate. In contrast, if the central barrier is high, i.e., \(\alpha\) is large, low-lying excited states below the central barrier are twofold degenerate: one state \(\psi_{\text{L}}\) is localized into the left (\(x<0\)) region while the other state \(\psi_{\text{R}}\) is localized into the right (\(x>0\)) region, and \(\psi_{\text{L}}\left(x\right)=\psi_{\text{R}}\left(-x\right)\) holds. Using a linear combination of these two degenerate states, one can recognize each state is degenerate the following two states: \(\psi_{\pm}\left(x\right)=\left[\psi_{\text{L}}\left(x\right)\pm\psi_{\text{R} }\left(x\right)\right]/\sqrt{2}\), where \(\psi_{+}\) (\(\psi_{-}\)) is a positive (negative) parity state. According to the exact diagonalization, \(\alpha=1.0\) and \(1.25\) give non-degenerate ground and low-lying excited states and thus \(\psi_{j}\) is just a \(j\)-th excited state, while \(\alpha=2.0\) and \(3.0\) give ground and low-lying excited states, which are almost two-fold degenerer

\begin{table}
\begin{tabular}{l c c c} State & Energy & Epochs & Time per Epoch (\(\mu\)s) \\ \hline
0th & +0.499998 & 23419 & 516.238 \\
1st & +1.499991 & 25646 & 519.026 \\
2nd & +2.499986 & 23157 & 527.849 \\
3rd & +3.500193 & 37880 & 534.115 \\
4th & +4.500201 & 19101 & 542.224 \\ \end{tabular}
\end{table}
Table 6: Calculation summary of excited states for a one-body problem under the harmonic oscillator potential. Calculation is performed with \(\omega=1.0\).
ate: \(\psi_{0}\) and \(\psi_{1}\) correspond to the ground states and \(\psi_{2}\) and \(\psi_{3}\) correspond to the first excited states.

Figure 15 shows the wave functions of the ground state and first, second, third, and fourth excited states. Table 7 shows the summary of calculations. Basically, not only the ground-state but also low-lying excited-states wave functions and energies are successfully calculated, even for the degenerate states. It is not apparent which calculation gives, left-right bases (\(\psi_{\mathrm{L}}\) and \(\psi_{\mathrm{R}}\)), parity bases (\(\psi_{+}\) and \(\psi_{-}\)) or even general linear combinations. The DNN calculations for both \(\alpha=2\) and \(3\) obtained wave functions with the left-right bases, while it may depend on the initial condition. Note that the exact diagonalization for \(\alpha=2\) obtained wave functions with the parity bases, while wave functions with the left-right bases are plotted by using linear combinations in Fig. 15 to make a comparison with the DNN result easily.

### Two-body systems

Two-body one-dimensional harmonic oscillators are taken as the last examples. Here, \(x_{\mathrm{max}}=y_{\mathrm{max}}=5\) and \(M_{x}=M_{y}=256\) are used for the spatial mesh and two layers each of which contains \(32\) units are used for the DNN. Here, the inter-particle interaction is not considered. The exact wave functions for several low-lying excited states can be written as linear combinations of Eqs. (32):

\[\Psi_{0}\left(x,y\right) =\psi_{0}\left(x\right)\psi_{0}\left(y\right), \tag{37a}\] \[\Psi_{1}\left(x,y\right) =\frac{1}{\sqrt{2}}\left[\psi_{0}\left(x\right)\psi_{1}\left(y \right)+\psi_{1}\left(x\right)\psi_{0}\left(y\right)\right],\] (37b) \[\Psi_{2}\left(x,y\right) =\frac{1}{\sqrt{2}}\left[\psi_{0}\left(x\right)\psi_{2}\left(y \right)+\psi_{2}\left(x\right)\psi_{0}\left(y\right)\right],\] (37c) \[\Psi_{3}\left(x,y\right) =\psi_{1}\left(x\right)\psi_{1}\left(y\right), \tag{37d}\]

where the energy eigenvalue of \(\Psi_{0}\), \(\Psi_{1}\), \(\Psi_{2}\), and \(\Psi_{3}\) are equal to, respectively, \(1\), \(2\), \(3\) and \(3\) for bosonic systems and

\[\Psi_{0}\left(x,y\right) =\frac{1}{\sqrt{2}}\left|\begin{matrix}\psi_{0}\left(x\right)& \psi_{1}\left(x\right)\\ \psi_{0}\left(y\right)&\psi_{1}\left(y\right)\end{matrix}\right|, \tag{38a}\] \[\Psi_{1}\left(x,y\right) =\frac{1}{\sqrt{2}}\left|\begin{matrix}\psi_{0}\left(x\right)& \psi_{2}\left(x\right)\\ \psi_{0}\left(y\right)&\psi_{2}\left(y\right)\end{matrix}\right|,\] (38b) \[\Psi_{2}\left(x,y\right) =\frac{1}{\sqrt{2}}\left|\begin{matrix}\psi_{0}\left(x\right)& \psi_{3}\left(x\right)\\ \psi_{0}\left(y\right)&\psi_{3}\left(y\right)\end{matrix}\right|,\] (38c) \[\Psi_{3}\left(x,y\right) =\frac{1}{\sqrt{2}}\left|\begin{matrix}\psi_{1}\left(x\right)& \psi_{2}\left(x\right)\\ \psi_{1}\left(y\right)&\psi_{2}\left(y\right)\end{matrix}\right|, \tag{38d}\]

where the energy eigenvalue of \(\Psi_{0}\), \(\Psi_{1}\), \(\Psi_{2}\), and \(\Psi_{3}\) are equal to, respectively, \(2\), \(3\), \(4\) and \(4\) for fermionic systems. Note that the second excited states are twofold degenerate, \(\Psi_{2}\) and \(\Psi_{3}\), in both the bosonic and fermionic systems.

Figures 16 and 17, respectively, show the wave functions of the ground state and first and second excited

Figure 14: Optimized function \(u\) obtained by the DNN for the ground-state (0th), the first and the second excited states.

Figure 13: Wave functions of the ground and low-lying excited states under the harmonic oscillator potential. The top panel shows the exact wave functions and the bottom one shows the DNN wave functions. In order to make consistency for the phase factor, \(-\psi_{3}\left(x\right)\) is plotted for the exact wave function of the third excited state.
states. Table 8 shows the summary of calculations. Note that \(\left[\Psi_{2}\left(x,y\right)\pm\Psi_{3}\left(x,y\right)\right]/\sqrt{2}\) are plotted for the second excited states for exact solutions. Not only the ground-state but also low-lying excited-states wave functions and energies are successfully calculated even for two-body systems. In addition, as one-body problems, the numerical cost for a low-lying excited state is almost the same as that for the ground-state. Thus, this method to calculate low-lying excited states can work even for multi-body systems with a reasonable numerical cost.

## IV Summary

In this paper, we proposed a method to calculate the wave functions and energies of not only the ground state but also low-lying excited states of quantum multi-body systems using the deep neural network and the unsupervised machine learning technique. In order to calculate systems of many-particle systems of identical particles, a simple method of symmetrization for bosonic systems

\begin{table}
\begin{tabular}{l l l l l} \hline \hline \(\alpha\) & \(j\) & Energy & Epochs & Time per Epoch (\(\mu\)s) \\ \cline{3-5}  & & Exact diagonalization & Deep neural network & \\ \hline
1.0 & 0 & +0.869573 & +0.869706 & 29108 & 513.203 \\
1.0 & 1 & +1.661393 & +1.661685 & 35596 & 523.556 \\
1.0 & 2 & +3.543667 & +3.544327 & 55037 & 524.351 \\
1.0 & 3 & +5.665058 & +5.666010 & 20083 & 536.858 \\ \hline
1.25 & 0 & +1.417858 & +1.417886 & 30284 & 512.083 \\
1.25 & 1 & +1.725904 & +1.726013 & 36992 & 523.092 \\
1.25 & 2 & +3.717933 & +3.717949 & 30943 & 528.675 \\
1.25 & 3 & +5.424725 & +5.424850 & 24474 & 534.907 \\ \hline
2.0 & 0 & +2.762317 & +2.762333 & 23663 & 519.261 \\
2.0 & 1 & +2.762333 & +2.762343 & 24764 & 529.129 \\
2.0 & 2 & +7.988520 & +7.989654 & 19491 & 532.193 \\
2.0 & 3 & +7.990618 & +7.989601 & 18662 & 534.543 \\ \hline
3.0 & 0 & +4.214229 & +4.214253 & 20526 & 515.042 \\
3.0 & 1 & +4.214229 & +4.214284 & 20290 & 521.291 \\
3.0 & 2 & +12.526202 & +12.526214 & 51131 & 519.062 \\
3.0 & 3 & +12.526202 & +12.526278 & 12488 & 534.752 \\ \hline \hline \end{tabular}
\end{table}
Table 8: Calculation summary of excited states for a two-body problem under the harmonic oscillator potential. Calculation is performed with \(\omega=1.0\).

Figure 15: Wave functions of the ground and low-lying excited states under the double-well potential. The top panels show the exact wave functions and the bottom ones show the DNN wave functions.

\begin{table}
\begin{tabular}{l l l l l} \hline \hline \(\alpha\) & \(j\) & Energy & Epochs & Time per Epoch (\(\mu\)s) \\ \cline{3-5}  & & Exact diagonalization & Deep neural network & \\ \hline
1.0 & 0 & +0.869573 & +0.869706 & 29108 & 513.203 \\
1.0 & 1 & +1.661393 & +1.661685 & 35596 & 523.556 \\
1.0 & 2 & +3.543667 & +3.544327 & 55037 & 524.351 \\
1.0 & 3 & +5.665058 & +5.666010 & 20083 & 536.858 \\ \hline
1.25 & 0 & +1.417858 & +1.417886 & 30284 & 512.083 \\
1.25 & 1 & +1.725904 & +1.726013 & 36992 & 523.092 \\
1.25 & 2 & +3.717933 & +3.717949 & 30943 & 528.675 \\
1.25 & 3 & +5.424725 & +5.424850 & 24474 & 534.907 \\ \hline
2.0 & 0 & +2.762317 & +2.762333 & 23663 & 519.261 \\
2.0 & 1 & +2.762333 & +2.762343 & 24764 & 529.129 \\
2.0 & 2 & +7.988520 & +7.989654 & 19491 & 532.193 \\
2.0 & 3 & +7.990618 & +7.989601 & 18662 & 534.543 \\ \hline
3.0 & 0 & +4.214229 & +4.214253 & 20526 & 515.042 \\
3.0 & 1 & +4.214229 & +4.214284 & 20290 & 521.291 \\
3.0 & 2 & +12.526202 & +12.526214 & 51131 & 519.062 \\
3.0 & 3 & +12.526202 & +12.526278 & 12488 & 534.752 \\ \hline \hline \end{tabular}
\end{table}
Table 7: Calculation summary of a one-body problem under the double-well potential.
and antisymmetrization for fermionic systems were also proposed.

The obtained wave functions and energies are consistent with the exact solution. We found that the neural network is not necessarily large for one-body systems, which also enables us to analyze the internal structure of the deep neural network used. For instance, just only one hidden layer with four units is enough to describe the ground-state wave function of the harmonic oscillator. This can be understood by using the piecewise approximation with linear functions. We confirmed that our simple (anti)symmetrization method works for multi-body systems. The numerical cost per epoch for fermionic systems is almost the same as that for bosonic systems. The numerical cost is almost proportional to the number of spatial meshes since the sparse matrix representation is used. In addition, the numerical cost for a low-lying excited state is almost the same as that for the ground state.

The deep neural network has been applied to solve many-fermion systems where the ground-state wave function is assumed to be a Jastrow wave function [33; 34; 42]. The method proposed in this paper can be an alternative method to solve many-fermionic systems since the ansatz for the ground-state wave function is lenient, and the symmetrization and antisymmetrization are treated on an equal footing.

Since the numerical cost is not so large and our (anti)symmetrization is quite simple, this method can be an alternative method to calculate wave functions and

Figure 16: Two-body wave function for the ground and low-lying excited states under the harmonic oscillator potential for bosonic systems without the interaction. The exact wave function is shown in the top row.

Figure 17: Same as Fig. 16 but for fermionic systems.
energies of the ground and low-lying excited states, for instance, for the electronic structure of molecules and solids, for the nuclear structure of atomic nuclei including a tetra neutron [53; 54; 55], and for cold atoms [56].

At this moment, we only considered one-dimensional systems, while most problems interested are three-dimensional systems. In addition, spin components, or even isospin components for nuclear systems, are often important. The restricted Boltzmann machine has been applied to obtain the ground- and low-lying excited-states wave functions [57; 32; 58]. Since the input is discrete variables in the spin systems, the Boltzmann machine is suitable. such pioneering works may help to consider the spin (or isospin) components in this work. Such extensions are possible within our framework, and remain for future work.

As far as we know, all the calculations using the deep neural network for wave functions are static, while describing many phenomena including the interaction between matter and laser [59; 60], ion-cluster collision [61] heavy-ion collision [62] nuclear fission [63; 64], and fusion [65]. In order to describe such phenomena, time evolution from a state obtained by the deep neural network is also interesting, while it is left for a future study.

Finally, let us make a comment on the interpretation of the wave functions obtained in our work. As we have shown, thanks to the simplicity of the deep neural network, we could interpret the structure of the network easily. We found that replacing the softmax function with the ReLU activation provides a piecewise linear function which approximates the ground-state wave function. Since any wave function including those for excited states, which is naturally continuous, can be approximated by a piecewise-linear function, we intuitively conclude that the neural network representation can work for any physical quantum mechanical system in any dimensions. The physical meaning of the piecewise-linear functions is as follows. First of all, linear functions are solutions of the free Schrodinger equation with no potential term. So it is a good idea to start with linear functions in physical systems. Then the inclusion of the potential term in the Hamiltonian causes the curvature of the wave function. The curvature is determined by the interplay between the Laplacian and the potential term in the Hamiltonian. So, the kink structure of the wave function is dictated by the Hamiltonian. The kinks correspond to the ReLU activations, thus in effect, the nonlinearity in the Hamiltonian corresponds to the neural network structure. This reminds us of the work [66] in which the deep layers of the deep Boltzmann machine representing the ground-state wave functions of spin systems were interpreted as a Euclidean Hamiltonian evolution, or the work [67; 68; 69] in which the deep layers of the sparse neural network used for the AdS/CFT correspondence were interpreted as a bulk curved geometry. Further interplay between the sparsity of the interpretable neural networks and Hamiltonians of physical systems is to be discovered.

###### Acknowledgements.

 The authors acknowledge for the fruitful discussion with Haozhao Liang, Masaaki Kimura, and Hiroyuki Tajima. T. N. acknowledges the RIKEN Special Postdoctoral Researcher Program, the Science and Technology Hub Collaborative Research Program from RIKEN Cluster for Science, Technology and Innovation Hub (RC-STI), and the JSPS Grant-in-Aid for Research Activity Start-up under Grant No. JP22K20372. H. N. acknowledges the JSPS Grant-in-Aid for Scientific Research (C) under Grant No. JP19K03488. The work of K. H. was supported in part by JSPS KAKENHI Grant No. JP22H01217, JP22H05111 and JP22H05115.
