New successive variational method of tensor-optimized antisymmetrized molecular dynamics for nuclear many-body systems

We recently proposed a new variational theory of"tensor-optimized antisymmetrized molecular dynamics"(TOAMD), which treats the strong interaction explicitly for finite nuclei [T. Myo et al., Prog. Theor. Exp. Phys. 2015, 073D02 (2015)]. In TOAMD, the correlation functions for the tensor force and the short-range repulsion and their multiple products are successively operated to the AMD state. The correlated Hamiltonian is expanded into many-body operators by using the cluster expansion and all the resulting operators are taken into account in the calculation without any truncation. We show detailed results for TOAMD with the nucleon-nucleon interaction AV8$^\prime$ for $s$-shell nuclei. The binding energy and the Hamiltonian components are successively converged to exact values of the few-body calculations. We also apply TOAMD to the Malfliet-Tjon central potential having a strong short-range repulsion. TOAMD can treat the short-range correlation and provided accurate energies of $s$-shell nuclei, reproducing the results of few-body calculations. It turns out that the numerical accuracy of TOAMD with double products of the correlation functions is beyond the variational Monte Carlo method with Jastrow's product-type correlation functions.

fully in the wave function without any truncation for the particle states. These 2p-2h excitations can describe the strong tensor correlation in nuclei via the coupling between 0p-0h and 2p-2h configurations. For the short-range correlation, it is in general difficult to express this correlation in the shell model based on the single-particle picture. We have combined TOSM with the central part of the unitary correlation operator method (UCOM) in which the central-type short-range correlation is explicitly treated [7,8]. The central-type shift operator is introduced in order to reduce the short-range amplitudes of the relative motion of a nucleon pair in nuclei. In UCOM, the transformed Hamiltonian is truncated up to twobody operators, while the exact transformation leads to many-body operators. The two-body approximation of UCOM is considered to be reasonable for the short-range correlation. In TOSM+UCOM, we explicitly treat the tensor and short-range correlations in the wave function. This method nicely describes the shell-model-like states with the correct order of energy levels in the p-shell nuclei [9][10][11][12].
On the other hand, nuclear clustering is an important aspect of nuclear structure, such as the triple-α Hoyle state in 12 C [13,14]. The clustering correlation is difficult to treat in the shell-model-type approach [11,15]. Recently, we have developed a new variational theory for the clustering description of nuclei using the N N interaction [16][17][18]. We employ antisymmetrized molecular dynamics (AMD) [19,20] as the basis state to describe the clustering states, and introduce two correlation functions of tensor-operator type and centraloperator type to treat the N N interaction. The correlation functions are multiplied by the AMD wave function and superposed with the original AMD wave function. We call this method "tensor-optimized antisymmetrized molecular dynamics" (TOAMD) [16]. The concept of TOAMD is similar to TOSM in treating the tensor and short-range correlations induced by the N N interaction. The scheme of TOAMD is extendable increasing the series of the multiple products of the two correlation functions by successive power expansion.
In TOAMD, there appear multiple products of the Hamiltonian and the correlation functions, which become the series of many-body operators in the cluster expansion. We treat all the resulting many-body operators, which makes TOAMD a variational method. In the previous work [17], we have described the s-shell nuclei with TOAMD within the double products of the correlation functions, and have shown that TOAMD nicely reproduces the results of the Green's function Monte Carlo (GFMC) using the AV8 ′ bare N N interaction [22]. The TOAMD provides a general formulation for nuclei with various mass numbers. The concept of TOAMD, which introduces the products of correlation functions successively, can be applied to the nuclear matter problem. Recently Toki and Hu have formulated tensoroptimized relativistic nuclear matter [21]. They proposed a variational framework for nuclear matter starting from the N N interaction, in which the correlation functions are multiplied successively by the plane-wave basis states. This framework is conceptually the same as TOAMD except for the basis functions.
We have also applied TOAMD to the description of the central-type correlation coming from the central N N interaction which has a short-range repulsion [18], in the same scheme as the previous work [17]. This subject gives scope to the application of TOAMD to other fields. We have focused on the central-type correlation and use the Malfliet-Tjon V (MT-V) central potential with a Yukawa-type tail and a strong short-range repulsion [23]. We have shown that the results of TOAMD for s-shell nuclei are good as compared with the few-body calculation. We have also compared the results of TOAMD with those using UCOM, which 2/25 showed that UCOM well describes the short-range correlations of nuclei quantitatively under the two-body approximation of the unitary transformation.
In this paper, we perform a detailed analysis of s-shell nuclei with the AV8 ′ bare N N interaction and the MT-V central interaction in TOAMD, the basic results of which are reported in Refs. [17,18]. We explain the detailed procedure of the TOAMD calculation for finite nuclei. We give the classification of the many-body operators in the cluster expansion of the correlated operators by using a diagrammatic representation. We also explain how to perform the energy variation including the multiple products of the correlation functions using Gaussian expansion in the TOAMD wave function. We report the properties of the TOAMD wave function, such as the spatial distribution of the correlation functions and the role of the many-body operators arising from the correlated Hamiltonian. The individual contributions of the many-body operators are discussed for each Hamiltonian component. In TOAMD, the wave function has a form of power series expansion with respect to the correlation functions, where each correlation function is determined independently in the energy minimization. This property of TOAMD is different from the ordinary variational approach with the Jastrow ansatz, in which the common correlation function is multiplied by every pair in nuclei. We discuss the detailed effect of the successive and independent optimization of the correlation functions on the solutions of TOAMD.
In Sect. 2, we explain the essential features of TOAMD. In particularly, we give the details of the cluster expansion leading to many-body operators. In Sect. 3, we present the results for s-shell nuclei. We show the results using the AV8 ′ bare interaction and the Malfliet-Tjon central interaction. A summary is given in Sect. 4.

Definition
We explain here the essential features of TOAMD. The details of TOAMD are given in Ref. [16]. We employ the AMD wave function as a reference state for TOAMD. The AMD wave function Φ AMD is given as the Slater determinant consisting of Gaussian wave packets of nucleons with mass number A as follows: The single-nucleon wave function φ( r) has a Gaussian wave packet with a range parameter ν and a centroid position D, a spin part χ σ and an isospin part χ τ . In the present study of s-shell nuclei, χ σ is fixed with the up or down component and χ τ is the proton or neutron. The range parameter ν is common for all nucleons. In this condition we can factorize the center-of-mass (c.m.) wave function from Φ AMD . In the shell-model limit with the condition of D i = 0 (i = 1, . . . , A) for all nucleons, the range ν has relation to ω as 2 ν/m = ω/2 where m is the nucleon mass. In TOAMD, we consider two kinds of correlations induced by the tensor force and shortrange repulsion, both of which are difficult to express in Φ AMD . Following the concept in Refs. [24,25], we introduce the pair-type correlation functions, F D for the tensor force and 3/25 F S for the short-range repulsion. We multiply them by the AMD wave function individually and superpose these components with the AMD wave function. The correlation functions F D and F S are determined variationally. This concept of TOAMD is similar to TOSM [5][6][7]. Here we define the two-body correlation functions as with a relative coordinate r ij = r i − r j . The pair functions f t D (r) and f t,s S (r) are the variational parameters and are explained later. The labels t and s stand for the isospin and spin channels of the two nucleons, respectively. The functions F D and F S change the relative motion of a nucleon pair in Φ AMD and do not excite the c.m. motion of Φ AMD . The function F D produces the D-wave transition owing to the tensor operator S 12 given as The functions F D and F S are scalar operators and do not change the angular-momentum state of Φ AMD . In general, two functions F D and F S are not commutable. Physically, the functions F D and F S can excite two nucleons in nuclei to the high-momentum state, which corresponds to the 2p-2h excitations in TOSM [5]. We multiply these correlation functions by the AMD wave function. In the present study of TOAMD, we include up to the double products of the correlation functions. We define the TOAMD wave function with the single correlation functions as We call this form of the TOAMD wave function "single TOAMD". As an extension of Eq. (6), we increase the order of the correlation function in TOAMD by adding the double products consisting of F D and F S following the previous study [16]. As with the single case, we define the TOAMD wave function with the double products of the correlation functions below, which we call "double TOAMD" as This form is based on the power series expansion in terms of the correlation functions F D and F S . It is noted that F D and F S in each term in Eq. (7) are independent and variationally determined. This indicates that we have five kinds of F D and five kinds of F S in the double TOAMD wave function. For simplicity, we denote the common symbols of F D and F S in Φ double TOAMD . This wave function of TOAMD has a general form with respect to mass number A and is commonly used for all nuclei and also for nuclear matter [21]. As an extension of Eq. (7), we can successively increase the order of power expansion with F D and F S to triple products such as F D F D F S when we want to increase the variational accuracy of the results. 4/25

Cluster expansion of the correlated operators
We use the Hamiltonian with a two-body bare N N interaction V for mass number A as Here, t i and T c.m. are the kinetic energies of each nucleon and the center-of-mass, respectively. We employ a bare N N interaction v ij AV8 ′ [1] consisting of central v C ij , tensor v T ij , and spinorbit v LS ij terms, which is used in the benchmark calculation of 4 He [22]. The total energy E in TOAMD is given as: where F stands for F D and F S . The operatorsH andÑ in the last equation are the correlated Hamiltonian and norm operator, respectively. We calculate the matrix elements of the correlated operators with the AMD wave function. The operatorsH andÑ consist of the various products of correlation functions, such as F † HF and F † F . These operators are individually expanded into a series of many-body operators in terms of the cluster expansion, the detailed procedure of which is given in Ref. [16]. In the case of the kinetic energy T , F † T F is expanded into many-body operators from two to five, with various combinations of particle index. For the two-body interaction V , F † V F is expanded into many-body operators from two-to six-body ones. In the same way, F † F † V F F gives up to ten-body operators.
In the calculation of Eq. (10), we use all the resulting many-body operators in the cluster expansion without any truncation. This is of importance to keep TOAMD as a variational theory. The procedure of the calculation of the matrix elements is performed systematically for any order of the multiple products of the correlation functions and also for any mass number. Later we explain the cluster expansion for the single TOAMD case.
In the Hamiltonian given in Eq.
We discuss typical cases of the cluster expansion of the correlated HamiltonianH and norm operatorÑ . We show all the diagrams to be calculated in the single TOAMD in Figs Fig. 1 Diagrams of the cluster expansion of F † F . Vertical lines indicate the particles numbering from the left side as 1, 2, 3, 4. Short horizontal lines indicate the two-body correlation function F . Each diagram represents the ladder (left), linked (middle), and unlinked (right) types. They are two-body, three-body, and four-body operators, respectively. We also write the configurations with symmetry factor below each diagram. and 3. We follow the rule of cluster expansion explained in Ref. [16], in which the two-body correlation function F is symbolically expressed as where the factor 1/2 is a symmetry factor to set the number of interaction pairs in F . Hereafter we write the condition of "i = j" as "i, j" in the summation of the particle index for simplicity. The notation of the configuration with square brackets [12] indicates the term of A i,j f ij , where the particle index i and j correspond to 1 and 2 in the configuration, respectively. We notice that the symmetry [12] = [21] from the particle exchange. In addition, we put the symmetry factor 1/2 in front of the configuration to express the function F in Eq. (12).
Following the rule of this expression for F , we expand the double product F † F , which appears in the correlated norm operatorÑ , into three terms as We again put the symmetry factor in front of the configuration with the square brackets. When the symmetry factor is unity, we do not show it. In general, the symmetry factor has a value of (1/2) N , where N is the number of symmetry with respect to the particle exchange in the configuration. When there is no symmetry, N = 0 and the symmetry factor becomes unity. In Eq. (15), the colon in the configurations represents the partition between the different correlation functions. In this case, we obtain three kinds of configurations: a two-body ladder operator, 1   34]. The four-body operator is the product of the two independent operators, which give asymmetry factor of 1/2 individually and 1/4 in total. These configurations can be visualized in terms of the diagrammatic representation shown in Fig. 1. Each configuration has a maximum number n in the square brackets, which leads to the n-body operator. This kind of cluster expansion is systematically performed in any power of the multiple products of the correlation functions [16]. The kinetic energy operator T is treated as one-body term except for the c.m. term T cm , which has already been given in Eq. (11). The configuration of the one-body operator T = A i t i is expressed as [1] with a single number in the square brackets. We take the case of F † T F . This correlated operator is expanded into nine diagrams from [12 : 1 : 12] to 1 4 [12 : 3 : 45], where the single numbers in the middle represent the one-body kinetic operator. The corresponding diagrams are displayed in Fig. 2.
For the two-body interaction V , we take the case of F † V F , which is expanded into 16 diagrams from 1 2 [12 : 12 : 12] to 1 8 [12 : 34 : 56]. The corresponding diagrams are displayed in Fig. 3. The same scheme of the cluster expansion limited to the single correlation function is presented using shell-model basis states [26,27].
In the double TOAMD given in Eq. (7), we consider the double products of the correlation functions in the wave function. Among the various kinds of correlated operators, we show the cases of the correlated operators of norm F † F † F F , one-body kinetic energy F † F † T F F and two-body interaction F † F † V F F as follows:  Fig. 3 Diagrams of the cluster expansion of F † V F . The dotted lines indicate the two-body interaction V . Table 1 Numbers of diagrams of many-body operators in the cluster expansion of 2  3  4  5  6  7  8  9  10   norm  1  13  46  47  25  6  1  --T  1  40  183  259  163  55  10  1  -V  1  40  295  587  516  235  65  10  1 The numbers of diagrams of many-body operators arising from the correlated norm, T and V are listed in Table 1. It is necessary to take all the resulting operators in order to retain the variational principle for TOAMD. In general, multiple products of many correlation functions produce a large number of many-body operators in the cluster expansion. Among these operators, the higher-body terms require larger calculation costs to obtain their matrix elements numerically, which often occurs for larger mass nuclei.

Energy variation
The present TOAMD wave function has two kinds of variational functions, the AMD wave function Φ AMD and the correlation functions F D and F S . We determine these functions using the Ritz variational principle with respect to the total energy in TOAMD as δE = 0 8/25 in Eq. (10). For Φ AMD , the centroid positions of the Gaussian wave packet of each nucleon (2) are variationally obtained using the cooling method [19]. The radial forms of F D and F S are optimized in four spin-isospin channels to minimize the total energy E. We adopt the Gaussian expansion method to express the pair functions f t D (r) in Eq. (3) and f t,s S (r) in Eq. (4) as follows: where a t n , a t,s n , C t n and C t,s n are variational parameters. We take the number of Gaussian functions N G = 7, in which we get converging solutions. For the Gaussian ranges a t n , a t,s n , we search for their optimized values in a wide range to cover the spatial correlation. The coefficients C t n and C t,s n are linear parameters in the single correlation terms of TOAMD. They are determined variationally by diagonalizing the Hamiltonian matrix. For the double products of the correlation function such as the F D F S term in Eq. (7), the products of the two Gaussian functions in Eq. (16) become the basis functions and, correspondingly, the products of C t n and C t,s n are the variational parameters. We express the TOAMD wave function in the form of a linear combination using the coefficients of the Gaussian expansion in the correlation functions: where the labels α and β are the set of the Gaussian index n and the quantum numbers of spin s and isospin t for two nucleons in the correlation function. The Hamiltonian and norm matrix elements are H α,β and N α,β , respectively. For the AMD wave function Φ AMD , we give it the labels α = 0 and β = 0, and the coefficientC 0 in the linear combination. The corresponding operators are given asH 0,0 = H andÑ 0,0 = 1. For the basis states with a single correlation function,C α indicates C t n and C t,s n , given in Eq. (16). For the double products of the correlation functions, these components are treated as single basis states with the two kinds of correlation functions. The corresponding expansion coefficients are denoted byC α , which can be C t Here α includes information on two kinds of correlation functions. It is noted that the coefficientC α determined in the calculation cannot be inversely decomposed into C t n and C t,s n . Finally, we solve the following eigenvalue problem to determine the total energy E and all the coefficientsC α in Eq. (17).
We briefly explain the calculation procedure of the Hamiltonian matrix elements in TOAMD, which are reduced to the matrix elements of the correlated Hamiltonian using the AMD wave function in Eqs. (10) and (18). We express the N N interaction V as a sum of Gaussians, similarly to the correlation function F . The correlated operatorsH andÑ involve the products of F and V . After the cluster expansion ofH andÑ to the many-body 9/25 operators, the resulting many-body operators include various combinations of the interparticle coordinates in the exponent of the Gaussian. This structure of the particle coordinates in the many-body operators makes it difficult to analytically evaluate the matrix elements in general. In the present approach, we use the Fourier transformation of the Gaussians in F and V [16,28]. This transformation decomposes the square of the interparticle coordinates r 2 ij in the exponent into the product of the plane waves, each having single particle coordinate r i and r j . In the momentum space, the matrix elements of the many-body operators result in the products of the single-particle matrix elements of the plane waves. Using the single-particle matrix elements in AMD, we perform multiple integration of the associated momenta in the last step and obtain the matrix elements of TOAMD. We explain the above procedure in Appendix A. Several typical cases are given in Ref. [16].

Central interaction
In TOAMD, the central-type correlation function F S is introduced to express the central-type correlation including the short-range repulsion. We investigate the applicability of TOAMD to the central interaction case only using the F S and F S F S terms in the double TOAMD wave function in Eq (7). We choose the Malfliet-Tjon V (MT-V) of N N interaction [23,29], which gives a strong short-range repulsion and a Yukawa-type tail. The MT-V potential is explicitly defined as in units of MeV and r in units of fm. The essential results of TOAMD for MT-V are given in Ref. [18] for s-shell nuclei in comparison with UCOM. Hence, we show the additional results in this paper.
In the AMD wave function, we optimize the range parameter ν, which is 0.11 fm −2 for 3 H and 0.25 fm −2 for 4 He in the results of the energy minimization. We also obtain that D i = 0 for all nucleons in two nuclei in the single TOAMD calculation. This indicates that the swave configuration is favored as the AMD wave function. We keep this condition throughout the analysis with MT-V. In Fig. 4, we show the results of TOAMD for two nuclei. We start from the results of AMD and further add the correlation terms successively. Here we simply denote the correlation function F S as "S". The symbol +S indicates the wave function of (1 + F S ) × Φ AMD for single TOAMD, and the symbol +SS indicates the wave function of (1 + F S + F S F S ) × Φ AMD for double TOAMD. We can see a nice convergence of energies in double TOAMD to the few-body calculations for two nuclei. In Fig. 4, at each correlationaddition step, the magnitudes of the kinetic and interaction energies increase together and their cancellation becomes large. The increase of the kinetic energies indicates that the short-range correlations arising from F S increase at each step.
In Table 2, we summarize the results of double TOAMD compared with other theories. It is found that the energies of 3 H and 4 He obtained in TOAMD nicely reproduce the fewbody results. In addition, the energies in TOAMD become lower than those of variational Monte Carlo (VMC) [30] for two nuclei. This result indicates that the accuracy of TOAMD is beyond that of VMC from the variational point of view. In VMC, they employ a Jastrow-type correlation function in which the two-body correlation function is multiplied for every pair  of nuclei. This approach is widely used for treating many-body theory in various fields. In this calculation, they essentially assume the common radial form of the correlation function for any pair. On the other hand, in TOAMD, we take the power series expansion of the correlation functions, and each correlation function at each order is treated independently and determined variationally with respect to the total wave function. This indicates that we can optimize the correlation functions of every term fully in TOAMD. We shall discuss later the advantage of the independent treatment of the correlation functions using the bare N N interaction AV8 ′ .

AV8 ′ potential
We discuss the results of TOAMD with the AV8 ′ potential for 3 H and 4 He. This part is the extended analysis reported in Ref. [17]. We report two kinds of the results based on single and double TOAMD. In Table 3, we list the range parameter ν of Φ AMD in the TOAMD wave function determined at the single and double levels. Similar to the case of MT-V central potential, we obtain D i = 0 for all nucleons in two nuclei, indicating the s-wave configuration of AMD for two nuclei [17].
In Table 4, we show the results of TOAMD, successively adding the correlation term stepby-step, where we use the ν-value optimized in double TOAMD. Here we use the labels D and S, indicating the correlation functions F D and F S , respectively. The symbol +D is the 11/25 The symbol +DD is the full calculation in double TOAMD defined in Eq. (7). The components of F S F D and F D F S haven an almost identical effect on the solutions; hence, the their results are combined and denoted by +SD+DS. We confirm that, starting from AMD, which gives positive energies, the energies of 3 H and 4 He are converged successively to the results of GFMC. This behavior is also shown in Figs. 12 and 13, explained later. The energy difference between GFMC and double TOAMD is 80 keV for 3 H and 1.2 MeV for 4 He.

Single TOAMD.
We first discuss the results of TOAMD with the single correlation function given in Eq. (6). This analysis provides the fundamental properties of the TOAMD wave function. After that, we proceed to the analysis of double TOAMD.
The energies of 3 H and 4 He are −5.34 MeV and −15.68 MeV, respectively, in single TOAMD with the optimized ranges ν listed in Table 3. It is interesting to plot the correlation functions f S (r) for the short-range part and f D (r) for the tensor part in Eq. (16), which are determined variationally using Gaussian expansion. We multiply the relative wave function of two nucleons φ rel ( r) in AMD with the correlation functions, which is physically meaningful. The relative wave function φ rel ( r) is defined in the form of a Gaussian wave packet as Here, D 1 , D 2 , and D r are the centroid of the two nucleon positions and thee difference between them, respectively; they are set to zero in the results of variation in TOAMD. This result makes φ rel spherical. In addition only the even channels are considered because of the s-wave configuration of AMD. The range parameter ν r is for the relative wave function of each nucleus. In Figs This indicates that f S (r) has two characters; one is the short-range correlation to avoid the repulsion in the interaction and the other is the intermediate correlation to obtain the attraction from the interaction. The results of 3 H and 4 He exhibit a similar trend. The reduction of short-range amplitude in SE is larger than that in the case of TE. This is related to a feature of the AV8 ′ potential: short-range repulsion of SE is stronger than that of TE [2]. The size of the intermediate amplitude in TE is larger than that in SE. This property is related to the presence of coupling between the S-and D-wave states via the tensor force in TE to obtain the total energy.
For the tensor correlation in Fig. 6, the TE channel is shown, where the term of r 2 is multiplied in relation to the definition of F D in Eq.  f S (r) as shown in Fig. 5. This agreement of the peak positions makes the tensor coupling increase in the TE channel. From the behavior of the correlation functions, the ranges of f S (r) and f D (r) are not short. This property increases the contributions of the many-body terms of the correlated operators in the cluster expansion beyond the two-body approximation. This point will be discussed later with the Hamiltonian components.
In Fig. 7, we plot the ν-dependence of the total energy and the Hamiltonian components for 3 H. Here, small ν represents the spatially extended AMD wave function and large ν represents the spatially compact case. It is noted that the ν-dependence corresponds to the ω-dependence in the shell-model prescription, because of the s-wave configuration of the AMD wave function. We clearly confirm the saturation behavior of the total energy E. Among the Hamiltonian components, the contribution of the tensor force is larger than that of the central force in the AV8 ′ potential for any range of ν. As the range of ν becomes large, the contributions of the interaction energy and also the kinetic energy increase, which causes a large cancellation between them to obtain the total energy.
We discuss the explicit role of the many-body operators in the cluster expansion of the correlated Hamiltonian. In Fig. 8, we decompose the total energy and the Hamiltonian components into many-body terms, while we sum up all the cluster components in the correlated norm part appearing in the denominator in Eq. (10). The symbols [2] and [3] represent the two-body and three-body cluster terms of the correlated operators, respectively. The symbol [1+2] appearing in the total and kinetic energies represents the sum of the one-body kinetic energy T , in which the c.m. term is subtracted and the two-body term. This [1+2] term corresponds to the calculation under the two-body approximation of the transformed Hamiltonian such as UCOM. For the total energy E, it is found that the [1+2] term cannot provide the saturation behavior around the energy minimum point at ν = 0.14 fm −2 . On the other hand, the three-body term [3] shows the saturation at small ν value of 0.10 fm −2 and becomes repulsive at large ν. The sum of all the cluster terms gives the saturation behavior 14/25  The symbols [2] and [3] indicate the two-body and three-body terms in the cluster expansion, respectively. For the kinetic energy, its half-value (K/2) is shown. The term of [1+2] in the total and kinetic energies represents a sum of the one-body kinetic energy and the two-body term.
of the total energy as a function of ν. This decomposition indicates the important role of the three-body term in the correlated Hamiltonian.
We show the contribution of each Hamiltonian component of 3 H in Fig. 8. For the kinetic energy (K), central (C) and tensor (T) forces, we decompose their correlated operators into each cluster term. For the kinetic energy, the terms [1+2] and [3] are shown. It is found that the [1+2] term makes a dominant contribution and the thee-body term [3] makes a the small contribution in scale. For central and tensor forces, the dominant contributions come from the two-body term [2], which is common, and the three-body term [3] gives the small contribution in scale. These results are similar to those of the kinetic energy. From the decomposition, all those up to the two-body term make the main contribution in the correlated Hamiltonian, but these are largely canceled out between the kinetic and interaction energies. The three-body terms of each Hamiltonian component commonly give the small value in magnitude. In total, owing to the large cancellation between the kinetic and interaction energies, the three-body term [3] can be compared to the [1+2] term in the total energy as shown in Fig. 8.    We proceed to the analysis of 4 He. In Fig. 9, we plot the ν-dependence of the Hamiltonian components of 4 He. We clearly confirm the saturation behavior for the total energy, where the optimized range ν = 0.17 fm −2 is larger than the value of 3 H of 0.14 fm −2 in Fig. 7, because 4 He is spatially more compact than 3 H. Among the Hamiltonian components, the 16/25 contribution of the tensor force is larger than that of the central force in the AV8 ′ potential, similar to the results of 3 H.
In Fig. 10, we decompose the total energy E into the many-body terms of the correlated Hamiltonian in the cluster expansion as [1+2], [3], and the newly four-body term [4]. It is found that the [1+2] term cannot present the saturation behavior at around the optimized range ν in the total energy. The three-body [3] and four-body [4] terms provide the contributions, which are not so small. Owing to the presence of higher-body terms beyond the two-body one, the total energy has a proper energy minimum, as any variational method should provide.
We see the cluster contributions of each Hamiltonian component of 4 He in Fig. 10 up to the four-body cluster term. For the kinetic energy, the [1+2] term makes the main contribution, and the three-body term [3] makes a smaller contribution in scale than that of the [1+2] term. The four-body term [4] makes a much smaller contribution in scale than the values of [1+2] and [3]. For the central and tensor forces, the dominant contributions commonly come from the two-body term [2], the three-body term [3] makes a small contribution in scale, and the four-body term [4] makes a much smaller contribution than the others. From these results, all those up to the two-body term make the main contribution in the correlated Hamiltonian. The three-and four-body terms of each Hamiltonian component commonly give small values. The higher-body terms tend to make smaller contributions in each component of the correlated Hamiltonian. On the other hand, the contributions from the kinetic and interaction energies are largely canceled out, which results in the small total energy. In total, the sum of the three-body and four-body terms in the total energy E is comparable to the sum of the two-body term, as shown in Fig. 10.
In the summary of single TOAMD, we confirm the importance of many-body cluster terms in the correlated Hamiltonian. This result indicates the inevitable role of many-body operators in the correlated Hamiltonian to describe the energy minimum required from the variational point of view. A similar result is obtained with the Brueckner-Bethe-Goldstone approach for nuclear matter [32], in which the three-body correlations induced by the Gmatrix are necessary to explain the energy saturation in nuclear matter. In TOAMD, the higher-body terms tend to make a smaller contribution in the cluster expansion of the correlated operators, but we cannot ignore these terms because of the large cancellation between the kinetic and interaction energies. The contribution of the higher-body terms is related to the spatial range of the correlation functions F S and F D , which is not short in TOAMD, as shown in Figs. 5 and 6. The present analysis also indicates that the two-body approximation should be treated carefully in relation to the range of the correlation functions in many-body theory.

Double TOAMD.
We discuss the results of double TOAMD. We include the double products of the correlation functions in TOAMD given in Eq. (7). We keep the nucleon positions D i = 0 for all nucleons, which are obtained in the results of single TOAMD.
In Fig. 11, we show the ν-dependence of the energy of 3 H in double TOAMD by adding the terms successively. It is confirmed that the energy surface of 3 H becomes deeper and flatter by adding the correlation terms in the wave function. This ν-independence corresponds to the ω-independence. The flat property of the energy curve represents the flexibility of the correlation functions, which are optimized as much as possible at any range of ν to include 17/25   the correlations in the TOAMD wave function. As a result the correlation functions gives a stable solution of TOAMD with respect to the range parameter ν.
In Fig. 12, we show the total energy E of 3 H obtained by adding the single and double correlation functions one by one with the solid circles. In each calculation, we set the value of ν to 0.095 fm −2 , as shown in Table 3, but the correlation functions are optimized in each calculation. The final energy with F D F D (+DD) is −7.68 MeV, as shown in Table 4. The Hamiltonian components are shown in Table 5. We can see the good agreement between TOAMD and other theories for each component. The matter radius is obtained as 1.746 fm in TOAMD. Figure 12 also explains the convergence of the energy, which is the same as that shown in Table 4. We confirm the converging behavior of the total-energy curve toward the values of GFMC. Similarly, we discuss the contributions of the kinetic energy (K), central (C), tensor (T), and LS (LS) forces in Fig. 12. The contributions of the tensor and LS forces appear after adding the tensor correlation (+D), which leads to the mixing of the D-wave component. We can see a nice convergence of each Hamiltonian component toward the results of few-body calculations (FB). From these results, TOAMD is able to treat the N N interaction explicitly owing to the two kinds of correlation functions F D and F S and is regarded as a successive variational method for nuclear many-body systems. We discuss the case of 4 He. The solid circles in Fig. 13 show the total energy obtained by successively adding the correlation functions in the TOAMD wave function in Eq. (7). The behavior is very similar to that of 3 H. We see a good convergence with the correlation functions. Finally we obtain −24.74 MeV of the energy in double TOAMD. The matter radius is obtained as 1.497 fm in TOAMD, reproducing the GFMC value of 1.490 fm [22]. In Fig. 14, we show the D-state probability of 4 He. This quantity is defined as the component of the wave function in which the total orbital angular momentum and total intrinsic spin are both two and they are coupled to be zero. In double TOAMD, the D-state probability is 13.04 which is close to the value of 13.91 obtained in the Faddeev-Yakubovsky (FY) calculation [22]. We expect that the small difference will be recovered by adding the higherorder terms in TOAMD such as the triple products of the correlation functions including F D , similarly to the Hamiltonian components.
In TOAMD, we emphasize that the correlation functions F D and F S are optimized independently in each term of Eq. (7). This is different from the Jastrow approach in which they assume the same functional form of the correlation function for every pair. It is interesting to examine the effect of the independent treatment of the correlation functions on the solutions of TOAMD. For this purpose we perform the following calculation: First, We also perform the calculation with a different condition from above, in which, after the full calculation of double TOAMD, we extract F S and F D from the single correlation function terms. These Including the double TOAMD components for both the MT-V central potential and the AV8 ′ bare potential, the binding energies and the Hamiltonian components of 3 H and 4 He are successively converged to values very close to those in the few-body calculations. It is found that the accuracy of the results obtained in TOAMD depends on the interactions. For the difference between the TOAMD and few-body calculations, the amounts of the difference using AV8 ′ , shown in Table 5, are larger than those using MT-V, shown in Table 2. This result suggests that more correlations should be included successively in TOAMD in the case of bare interactions, which is possible by increasing the multiple products of the correlation functions including F D . We shall consider the triple case, such as F D F D F S and F D F S F S in subsequent work; this is expected to increase the numerical accuracy of TOAMD using the bare interaction. The other way to extend TOAMD is the multi-basis representation of the AMD wave function from the single basis case.

Summary
We have performed a detailed analysis of s-shell nuclei with a new variational theory of "tensor-optimized antisymmetrized molecular dynamics" (TOAMD) [16][17][18]. In TOAMD, we introduce two kinds of correlation functions, F D and F S for the tensor and short-range correlations, respectively to treat the strong interaction directly in the description of nuclei. We employ AMD as a reference state and multiply the correlation functions by the AMD basis state. The scheme of TOAMD is extendable by increasing the series of the multiple products of the correlation functions by successive power expansion. It is noted that each correlation function included in each term of the multiple products is treated independently and their functional form can be different. This property of TOAMD is different from the ordinary Jastrow method, in which the correlation functions are multiplied by every particle pair in a common form. The present formulation of TOAMD is applicable to all nuclei with various mass numbers.
In TOAMD, the product of the correlation function and the Hamiltonian is taken into account as the correlated Hamiltonian and is treated in terms of the cluster expansion. The correlated Hamiltonian becomes a series of many-body operators. We explicitly manipulate these operators in the calculation of the matrix elements in TOAMD without any truncation of the higher-body operators beyond the two-body one. This is an important point to retain the variational principle for the TOAMD wave function. In this paper, we employ diagrammatic representation to distinguish each diagram obtained in the cluster expansion of the correlated operators.
We have shown the results of s-shell nuclei with TOAMD including up to the double products of the correlation functions (double TOAMD). We employ two kinds of N N interactions of the MT-V central potential with strong short-range repulsion and the AV8 ′ bare potential with tensor and LS forces. The energies and Hamiltonian components obtained in TOAMD nicely reproduce the few-body results of 3 H and 4 He such as GFMC in both interactions. These results show the power of TOAMD and the efficiency of two correlation functions F D and F S to treat the N N interaction. 21/25 We show the spatial distributions of the correlation functions F D and F S , which are found to be not short-range. In particular F S has two roles in the reduction of the short-range amplitude of the nucleon pair and the attraction at the intermediate range, where the latter is related to F D via the S-D coupling in the triplet-even channel. We also discuss the role of many-body terms of the correlated Hamiltonian in the cluster expansion, in particular, more than the two-body term. It is found that the contributions of three-body and four-body terms are smaller than that of the two-body term, but not negligible in each Hamiltonian component of the kinetic and interaction energies. Owing to the large cancellation between the kinetic and interaction energies, the contribution of the higher-body terms survives and becomes not small in the total energy. In order to obtain the proper energy saturation variationally, it is necessary to include higher-body terms in the correlated Hamiltonian.
We investigated the advantage of the independent treatment of the correlation functions in TOAMD. We perform constraint calculations in which the correlation functions are fixed to be common for every term of TOAMD. This leads to an energy loss of few MeV for s-shell nuclei. This result indicates that the correlation functions should be different at each order of the power series expansion in TOAMD. From the variational point of view, this treatment is better than the case using common correlation function for all orders. In fact, for the central MT-V potential, it is found that the numerical accuracy of TOAMD with double products of the correlation functions is beyond the variational Monte Carlo calculation using the product-type common correlation functions.
In future work, we shall increase the multiple products of the correlation functions to the triple case, which is expected to increase the numerical accuracy of TOAMD. The other way to extend TOAMD is the superposition of many AMD basis states. Based on the success for the s-shell nuclei in the present analysis, we shall apply TOAMD to the p-shell nuclei. We also have a plan to treat the three-nucleon interaction explicitly, e.g., the Fujita-Miyazawa type, whose treatment is essentially similar to the case of the many-body operators of the correlated Hamiltonian in TOAMD.

A. Matrix elements in TOAMD
We explain briefly how to calculate the matrix elements of the operatorÔ with the TOAMD wave function Φ TOAMD ; these are equivalent to the matrix elements of the correlated operator O with the AMD wave function Φ AMD as O =Ô + F †Ô +ÔF + F †Ô F + · · · . (A2) The correlated operatorÕ consists of the original operatorÔ, the multiple products ofÔ and the two-body correlation function F . Each term including F in Eq. (A2) is expanded into a series of n-body operators in terms of the cluster expansion with n = 2, · · · , A. We consider the case of the specific n-body operatorÕ [n] , which has the following form:Õ i1i2...in with symmetry factor S [n] for particle exchange and n-particle indices of i 1 , · · · , i n . The matrix elements ofÕ [n] with the AMD wave function are given as Φ AMD |Õ [n] |Φ AMD = S [n] A i 1 ,i 2 ,...,in j 1 ,j 2 ,...,jn where the quantity B ij is the single-nucleon overlap matrix element in AMD. We shall consider the interaction V as the operatorÔ. We express V and F in terms of the sum of Gaussian functions in coordinate space. In order to calculate the kernel of the matrix elements φ i1 · · · φ in |Õ [n] 1 2···n |φ j1 · · · φ jn in Eq. (A3), we use the Fourier transformation of the Gaussian for the interparticle coordinate r ij = r i − r j with range a in V and F [16,28]  d k e − k 2 /4a · e i k· ri e −i k· rj · k 2 S 12 (k) , (A7) The tensor operator S 12 is also transformed into the momentum space. Equation (A6) is used in the calculation of the LS term. The merit of this transformation is that the coordinate r ij in the Gaussian function becomes separable for each particle coordinate r i and r j with the plane-wave form in momentum space. We explain the case of the central interaction for V and the central correlation for F to calculate the kernel matrix elements. In the correlated interaction,Ṽ , the specific nbody operatorṼ [n] becomes the products of several Gaussian functions and has various connections with the interparticle coordinates, such as the diagrams shown in Fig. 3. We transform each Gaussian function inṼ [n] into the corresponding momentum space one by one with the Fourier transformation. We denote the transformed operator in momentum space asṼ [n] k which is a function of multi-momenta having Gaussian functions and includes the products of plane waves for each particle coordinate, as shown in Eq. (A5). The structure ofṼ [n] k for particle coordinates is schematically written as 1 ( r 1 ) · g [n] 2 ( r 2 ) × · · · × g [n] n ( r n ).
Each operator g i ( r i ) has a form of e i Ki· ri depending on the momentum K i , which is the sum of the momenta for which the underlying Gaussian functions in Eq. (A5) commonly involve the coordinate r i . The operator g i ( r i ) can have spin and isospin operators, which 23/25 are ignored at present. Hence the kernel matrix elements of the n-body operator are given by the products of the single-particle matrix elements of the plane wave in AMD as φ i1 · · · φ in |Ṽ [n] k,1 2···n |φ j1 · · · φ jn ∝ φ i1 |g [n] 1 ( r 1 )|φ j1 × · · · × φ in |g [n] n ( r n )|φ jn , where φ i |e i K· r |φ j = φ i |φ j · e i K·( D * i + Dj)/2− K 2 /8ν .
Finally we perform the multiple integration for all momenta analytically, in which the integrand includes the Gaussian functions for each momentum. For the kinetic energy and LS force, we take the explicit derivatives by the corresponding coordinate [16]. We explicitly give the form of g i ( r i ) in the case of F † F appearing in the correlated norm operator. The operator F † F gives two momenta of k 1 and k 2 from each F and is expanded into three diagrams of 1 2 [12 : 12], [12 : 13], and 1 4 [12 : 34] in the cluster expansion in Fig. 1. For the configuration 1 2 [12 : 12] as the two-body operator with ladder-type, g [2] 1 ( r 1 ) = e i( k1+ k2)· r1 , g [2] 2 ( r 2 ) = e −i( k1+ k2)· r2 .