Introduction

Quantum systems, from atomic nuclei to solids, are composed of many particles; thus, finding the states of many-body systems is an essential problem. In advanced quantum many-body methods, such as density functional theory [1-3], tensor network [4, 5], and variational Monte Carlo method [6-8], the ground-state energy can be treated in a variational form, turning finding the ground state of the system into a variational problem. Among the methods for solving the variational problems, unsupervised machine learning (ML) is gaining attention owing to its powerful performance for optimization, making it an ideal tool for solving many-body problems in quantum systems.

In unsupervised ML for variational problems, the ML structure can be treated as a trial function. In the case of quantum many-body systems, the energy expectation value with respect to the system Hamiltonian H is treated as the loss function, and the output of the ML is usually regarded as the wave function.

For unsupervised ML applications in spin systems, restricted Boltzmann machines (RBMs) were first employed as a wave function ansatz to address both static and time-evolution problems [9]. Subsequently, improved algorithms were developed to enhance accuracy [10], and deterministic time-evolution schemes were introduced [11]. Additionally, the Jastrow-Slater wave function ansatz was adapted for unsupervised ML [12]. The excited-state calculation was presented with two different ansätze [13, 14].

In the field of nuclear physics, ML has emerged as a crucial tool with extensive and impactful applications [15, 16]. These applications span the prediction of both nuclear ground- and excited-state properties, including nuclear mass [17–21], charge radius [22, 23], excited states [24, 25], α decay [26, 27], β decay [28, 29], charge density [30–32], density functional [33], nuclear level density [34, 35], ground-state magnetic moments [36], photoabsorption cross section [37], single-Λ hypernuclei [38], and the distribution of ground-state spin in the two-body random ensemble [39]. In statistical mechanics, the transfer matrix for a spin-glass model was calculated using an unsupervised deep neural network [40].

For many-body systems of bosons or fermions depending on continuous coordinates, incorporating the required (anti) symmetrization into the neural network wave function presents a significant challenge. In Ref. [41], a single-layer neural network was used to calculate deuterons in the momentum space with a pre-training process. In Ref. [42], further analysis was conducted on the structure of neural networks and numerical uncertainty. For bosonic systems, Ref. [43] introduced a neural network wave function for studying the Calogero-Sutherland model and Efimov bound states, while Ref. [44] proposed incorporating a pooling layer to enforce symmetrization. For fermionic systems, such as electronic structure of atoms and molecules and nuclear structure with an ab initio Hamiltonian, the Jastrow-Slater ansatz is often introduced to consider the fermion antisymmetrization and the interparticle correlation [45–51]. The hidden-fermion technique [52], which is an extension of the Jastrow–Slater ansatz, was proposed and applied with an unsupervised ML. This technique was immediately applied to nuclear systems to search for the ground-state energy of ¹⁶O [53], and both energies and other ground-state properties of A ≤ 20 nuclei [54]. Another proposed method for calculating the electronic structure is to map the fermionic system into a spin system [55]. Recent progress is summarized in Ref. [56].

Among these methods, in Ref. [44], bosonic symmetrization is considered using the pooling layer, whereas it does not work for fermionic systems. A Jastrow-Slater wave function is applied in many works to solve fermionic systems [45–51]; however, it does not work for bosonic systems.

In Ref. [57], another method using an unsupervised deep neural network (DNN) was proposed, which can solve both bosonic and fermionic systems in coordinate space. The calculation steps are the same for both bosonic and fermionic systems, differing only in the sign in the trial wave function to account for (anti)symmetrization. According to the universal approximation theorem [58, 59], a multilayer neural network can approximate any function with the desired accuracy by changing the size of the neural network. Therefore, a neural network can be used as a trial wave function [9–11, 41–43, 56, 57]. By utilizing the energy expectation as the loss function for minimization, this method successfully calculated the ground state of one-dimensional (1D) one- to three-body systems. However, the spin and isospin degrees of freedom, which are crucial for discussing the properties of magnetic materials and atomic nuclei, have not yet been considered in this method.¹

In this study, we extend the method proposed in Ref. [57] to include spin and isospin degrees of freedom. We will focus on two-body systems without external potential as the first step. Owing to translational symmetry, only the relative motion was considered. This setup corresponds to the onebody calculation in Ref. [57]. This method is validated by calculating the deuteron, the simplest realistic many-body system. Its Hamiltonian contains a noncentral tensor force, making the inclusion of spin and isospin degrees of freedom essential. Because a deuteron has only one bound state, we do not focus on the excited states in the present discussion.

Note that the system we consider is the same as that in Refs. [41, 42]. However, our method is technically different from theirs. In the process of optimization, Ref. [41] used a supervised method to pre-train DNN parameters, facilitating faster convergence and ensuring accurate behavior at the origin [42, 60]. Our method does not require pre-training, primarily because of the choice of the DNN ansatz for the wave functions. This choice influences the behavior at the origin (discussed in the Appendix), along with other technical differences, such as mesh point distributions and optimizers. It was pointed out in Ref. [42] that the structure of a neural network affects the convergence, where a deeper neural network yields a lower convergence rate. In contrast, our DNN model always converges with different numbers of layers by ensuring a relatively smaller number of hidden nodes. In the results, overfitting is not observed in our DNN wave functions, whereas it appeared in Ref. [42]. In Ref. [41], it was also pointed out that the extension to arbitrary spin and isospin was a future perspective; our work derived a general formula for arbitrary spin and isospin and calculated the deuteron with 18 possible partial waves by applying a nonfully connected DNN structure.

The remaining of this paper is organized as follows. In Sect. 2, we discuss the theory to include the spin and isospin degrees of freedom, and also the center of mass to the previous method with partial wave expansions. The corresponding DNN model is proposed in Sect. 3. In Sect. 4, we discuss the calculation results of the deuteron. Finally, we provide a summary and future perspectives in Sect. 5.

Method

Neural Network Models

3.2

The design of the non-fully-connected neural network

The design of the present non-fully-connected DNN, which generates partial wave functions with corresponding coefficients djϕj, is shown in Fig. 1. The input is the relative distance ri. The outputs were the partial wave functions. After determining the cutoff value of the orbital angular momentum L and the spin quantum number S, the possible angular momentum sets {LSJ} of partial waves are then determined by the LS coupling. For each single output, the nodes in the hidden layers are fully connected, while the hidden layer nodes of different outputs are not connected to each other. The “softplus” function(28)is used as the activation function.

Fig. 1Full-screen View

Schematic of a non-fully connected deep neural network. The shared input is the relative distance between the particles. The output are the partial wave functions with coefficients djϕj, where j (j = 1, 2, ..., q) represents the combination of the angular momenta, L, S, and J. There is no connection between the hidden layer nodes of different outputs. The computational accuracy can be enhanced by increasing the number of layers and nodes in each layer

Throughout the training process, the relative distances of the particles are discretized into M uniformly distributed lattice points, which are passed to the DNN as a dataset. The output of the DNN is a dataset containing all partial wave functions. The energy expectation value is regarded as the loss function, which is calculated using Eq. (27). The parameters of DNN are then updated by the Adam optimizer [71] and the ML is implemented by using the Tensorflow [72].

Calculation results for deuteron

The feasibility of the method was validated by applying it to deuterons. To include a more general case, at the beginning, we consider the orbital angular momentum quantum number L from zero to four and spin quantum number S from zero to one in Eq. (24) with 18 possible states. The 18 partial wave functions are computed in a 20 fm box with 1500 mesh points. A non-fully connected DNN is used, which is composed of two layers, each containing 16 units. Moreover, we also tested a fully connected version of such a DNN to explore the impact of connections between nodes on the training results. Both models are trained to testify whether only the ³S₁ and ³D₁ states are obtained, which is a known fact regarding the properties of deuterons [73]. Table 1 lists the percentage of each state in the wave function. The results show that only the ³S₁ and ³D₁ states make significant contributions in the ground state, while the contributions of the other states are sufficiently small to be neglected, which is consistent with the known fact. It was also found that both the fully connected and non-fully connected DNN can generate the ground state with the same level of accuracy, showing that the connection between nodes does not make a significant difference in precision. We chose the non-fully connected one for subsequent calculations. The relative error of the ground-state energy is 0.05% of the benchmark value obtained by the Green’s function Monte Carlo calculation [61].

Since it is proved that only the ³S₁ and ³D₁ states contribute to the ground-state energy, we can reduce the number of outputs from 18 to 2 to reduce the size of the DNN. Accordingly, the calculation cost is also reduced; hence, hereinafter we apply a more complicated hidden layer structure that is composed of three layers each containing 16 units with 30 fm box size and 2000 mesh points. The energies obtained using the AV18, AV8′, and AV4′ potentials are listed in Table 2. In order to see the effect of the electromagnetic interaction between a proton and a neutron, calculations using the AV8′potential with and without the electromagnetic interaction are also shown. The accuracy is both within 2 keV (0.1%) for AV18 and AV8′ calculations compared to the benchmark [61]. The AV4′ potential comprises only four terms; hence, it yields a final error of approximately 1% [70]. In Fig. 2, the deuteron wave functions are shown based on AV18 and AV8′ potentials, which are compared with deuteron wave functions of the AV18 calculation in Ref. [61]. It is observed that the DNN results of both AV18 and AV8′ calculations reproduce the AV18 benchmark [61] quite nicely. The relative errors and are shown in Fig. 3, where u and w are the partial wave functions of ³S₁ and ³D₁ states, respectively. Because the contribution of the ³D₁ state in the loss function is much smaller than that of ³S₁, its relative error is accordingly larger.

Fig. 2Full-screen View

Deuteron wave functions with the AV18 and AV8′ potentials, where u and w denote the S -wave and D -wave components, respectively. The benchmark is based on the AV18 potential [61] taken from Ref. [74]. The hidden layers of each output are composed of three layers, each of which has 16 units. Both results the AV18 and AV8′ potentials show good agreement with the benchmark

Fig. 3Full-screen View

Relative error of DNN deuteron wave functions to the AV18 benchmark results [61, 74]. The hidden layers of each output are composed of three layers each of which has 16 units

The DNN energies with respect to different numbers of mesh points in the AV18 potential are listed in Table 3. It was found that 500 mesh points were sufficiently large to produce results with less than 1% relative error, with the fluctuation of the DNN training results becoming an important factor in the precision.

Table 4 shows the performance of the program by varying the number of layers and nodes with the AV18 potential [74-76]. In the case of a single hidden layer containing 16 units, the relative error in energy is 2.857%. With an increase in the number of units to 32 or the addition of a hidden layer, the relative errors do not exceed 0.1%. A more complicated network structure does not result in higher precision, as the training results fluctuate slightly every time. This finding differs from that of Ref. [42], where adding a second layer induces a sharp decrease in model accuracy, whereas increasing the number of hidden nodes partially compensates for this loss. This shows that our DNN model is more stable and can be effectively extended to model complex physical systems, where expressive power necessitates the use of deeper network architectures.

In Fig. 4, the relative errors of the loss function with respect to the benchmark energy are shown as functions of epochs in the AV18 calculation for two optimizers. When the Adam optimizer was used, spikes appeared in the loss values, although the loss still decreased between them. In contrast, when we use the stochastic gradient descent (SGD) optimizer [77], these spikes do not appear; instead, it requires more steps to reach convergence.

Fig. 4Full-screen View

Relative error of with respect to the benchmark AV18 ground-state energy E_gs = -2.22458 MeV [75, 76] as functions of the number of epochs. The hidden layers of each output were composed of three layers, each of which has 16 units. The DNN with the Adam optimizer converges around 20000 epochs, while the one with the SGD optimizer fails to converge within 10⁵ epochs

Summary

In Ref. [57], an unsupervised machine learning technique was developed to calculate the ground state using a deep neural network. In this study, we extended the method by introducing spin and isospin degrees of freedom through the introduction of partial wave expansions generated by a non-fully connected deep neural network.

The method was verified by calculating the simplest two-body nuclear system—deuteron. At first, the partial waves with the orbital angular momentum quantum number L from zero to four and the possible total spin S from zero to one were calculated. The results were consistent with the fact that only the ³S₁ and ³D₁ states contribute to the deuteron ground state. In the following calculations, the obtained wave functions and energies showed consistency with the benchmark [61, 70]. We found that the deep neural network does not need to be large. In the present case, two hidden layers, each containing eight units, were sufficient to generate faithful representations of the wave function.

We believe that these improvements to the deep neural network approach hold promise for further studies that could extend the methodology of this work. For example, one can extend from two-body calculations to N -body calculations. For three-body systems, we believe the ground state can be obtained by solving the Faddeev equation using hyperspherical coordinates [78], in which the corresponding partial wave expansion includes two angular terms with spin and isospin components. Based on the present work and our deep neural network approach to solve the Dirac equation [79], another promising direction is to incorporate the spin and isospin degrees of freedom into the Dirac framework. We anticipate that the present deep neural network architecture can be extended to generate Dirac partial-wave functions.

DNN output analysis

We show the wave functions corresponding to the two strategies for selecting the DNN outputs: 1) the output corresponds to ξLSJ and 2) the output corresponding to ϕLSJ and ξLSJ is calculated from the output. In our test, we use 1700 mesh points in a 20 fm box to generate ³S₁ and ³D₁ states. We find that the DNN can generate more precise results of ϕLSJ than ξLSJ close to the origin, as shown in Fig. 5. Therefore, ϕLSJ was utilized as the DNN output in our calculation.

Fig. 5Full-screen View

Wave functions ϕLSJ obtained by 1) the output corresponds to ξLSJ (Blue) and 2) the output corresponding to ϕLSJ and ξLSJ is calculated from the output (Red). The benchmark calculation [61] is also shown (Green)