Tensor-optimized antisymmetrized molecular dynamics in nuclear physics

We develop a new formalism to treat nuclear many-body systems using bare nucleon-nucleon interaction. It has become evident that the tensor interaction plays important role in nuclear many-body systems due to the role of the pion in strongly interacting system. We take the antisymmetrized molecular dynamics (AMD) as a basic framework and add a tensor correlation operator acting on the AMD wave function using the concept of the tensor-optimized shell model (TOSM). We demonstrate a systematical and straightforward formulation utilizing the Gaussian integration and differentiation method and the antisymmetrization technique to calculate all the matrix elements of the many-body Hamiltonian. We can include the three-body interaction naturally and calculate the matrix elements systematically in the progressive order of the tensor correlation operator. We call the new formalism"tensor-optimized antisymmetrized molecular dynamics".


Introduction
The tensor interaction plays a very important role in nuclear physics. It originates from the pion exchange between nucleons, which is the most essential component in the nucleonnucleon interaction. The Green's function Monte Carlo (GFMC) simulation by the Argonne group for light nuclei demonstrated that the pion contribution on the nuclear binding energy is about 80% of the entire contribution of the two-body interaction [1]. The pion exchange interaction can be divided into tensor and central spin-spin components. If we take a phenomenological potential as the Argonne AV8 ′ potential [2,3], we see that the lightest nucleus, a deuteron, cannot be bound by the central interaction alone: the tensor interaction plays a dominant role in the binding of the deuteron [4]. Analyzing the energy contribution of the tensor interaction, we find that the transition matrix element between the s-wave and d-wave components provides the largest attraction for the deuteron binding.
This finding motivated the introduction of the tensor-optimized shell model (TOSM), where the TOSM wave function has a low-momentum shell-model state and also highmomentum 2-particle-2-hole (2p-2h) configurations, aiming at describing heavier nuclei and nuclear matter [5]. The 2p-2h configurations are excited by the tensor interaction from the important step is to take care of the antisymmetrization. For this difficult problem, we write all the interactions and correlations in the momentum space so that we are able to write all operators in separable forms in the particle coordinates [14]. We are then able to use the matrix technique to treat the antisymmetrization systematically. The third important ingredient is to treat all the necessary momentum integrations with a multiple of momenta in the Gaussian integrations. We introduce the source terms for the momenta and take derivatives of the integrated result of the fundamental Gaussian integrals. All calculations of the matrix elements are performed in the rectangular coordinates of particles using the above important ingredients for multi-body operators.
In light nuclei there appear cluster structures, which are difficult to describe in the shell + mean-field approximation. The cluster structure is better described in the AMD framework using effective interactions, where the 4 He wave function is written as (0s) 4 Gaussian functions. Although some trial was performed to include higher-spin states in the AMD framework for 4 He, it was not successful in treating the tensor interaction [15]. On the other hand, the TOSM was able to treat the tensor interaction in the shell-model basis by introducing high-momentum 2p-2h configurations. Use of the bare nucleon-nucleon interaction including the strong tensor interaction in the AMD framework is highly anticipated. This is achieved in the present TOAMD formulation, which is systematical and straightforward to apply to heavier-mass nuclei than those of the GFMC simulation. We would like to apply the TOAMD for the description of the coexistence of and competition between shell and cluster structures. One important application of the TOAMD is the alpha-condensed states in 8 Be, 12 C, and 16 O using the bare nucleon-nucleon interaction [16].
This paper is arranged as follows. In Sect. 2, we introduce the tensor-optimized antisymmetrized molecular dynamics (TOAMD), where the TOAMD wave function and the Hamiltonian are written explicitly for nuclear many-body systems. Here, we write all the matrix elements of the Hamiltonian in the AMD wave function so that we define all ingredients of the TOAMD theory. In Sect. 3, we calculate matrix elements of an interaction with one tensor correlation operator, which leads to multi-body operators up to four-body for the AMD wave function. In Sect. 4, we write all the necessary integrals and differential formulas for multiple-momentum integrations of Gaussian functions. In Sect. 5, we calculate two-body interactions with two tensor correlation operators as examples of using all the formula developed for the TOAMD theory. We describe a systematic method to write matrix elements for many-body systems. In Sect. 6, we introduce the short-range correlation operator as the sum of Gaussian functions. We explicitly give some of the matrix elements for a three-body interaction with the short range and tensor correlations. We present here a systematic method to calculate any complicated matrix elements in the TOAMD theory. In Sect. 7, we summarize the present paper. We further present two appendices. In Appendix A, we desribe the co-factor matrix theory to treat the antisymmetrization. In Appendix B, we explicitly give the Gaussian integrals for multi-body operators.

Tensor-optimized antisymmetrized molecular dynamics
We describe here the construction of the tensor-optimized antisymmetrized molecular dynamics (TOAMD). In this section, we give all the ingredients such as the TOAMD wave function, Hamiltonian, and the matrix elements of the Hamiltonian in the AMD wave function. 3/36

Wave function of TOAMD
In this subsection, we introduce the TOAMD wave function for finite nuclei. In the concept of the tensor-optimized shell model (TOSM), it is important to prepare a basic wave function to represent a correct density profile with low-momentum components, and add high-momentum components to be excited by the strong tensor interaction [5]. We take the following form for the TOAMD wave function: Here, |AMD is an AMD wave function for mass number A: where p i denotes various quantum numbers of a single nucleon and |p i = |ψ pi χ pi ξ pi . The antisymmetrizer A makes sure that the exchange of particle coordinates among all particles have an opposite sign from the original wave function (Slater determinant). We can include multiple tensor correlation operators F D · · · F D |AMD for the description of the multi-cluster states and short-range correlation to be discussed later. The spatial wave function of single nucleon ψ pi ( r i ) is written in terms of the shifted Gaussian function: The particle coordinate is written as r i for all nucleons i = 1, . . . , A. The size parameter ν and the position parameter D pi are the variational parameters to specify the spatial wave function. The position vectors D pi are in general complex variables, but we write them as real variables in this paper to simplify the notation. For numerical calculations it is important to take the AMD wave function in order to calculate matrix elements analytically using the Gaussian integral formula. As the nuclear system becomes heavy, we ought to include more Gaussian functions and additionally perform angular momentum and parity projections. The spin wave function χ pi (s i ) is written as: The spin wave function χ p (s) is written as a linear combination of spin-up | ↑ and spindown | ↓ wave functions, where β p is a complex variational parameter. As for the isospin part ξ pi (t i ), we take pure proton and neutron states: Here, |proton and |neutron are pure proton and neutron states, respectively. The tensor correlation operator is expressed as: with the variational function written by the sum of Gaussian functions: 4/36 where C µ and a µ are variational parameters. Here the tensor operator is written as: The spin operators σ i and σ j are for particles i and j. The isospin τ i · τ j operator is added for the tensor correlation operator, where the tensor correlation is very strong because its origin is one-pion exchange. Although we do not show it explicitly, we also consider the isospinindependent tensor correlation in the calculation. The relative coordinate r ij is the difference between the positions r i and r j of two particles i and j, andr ij is a unit vector with its direction. It is essential to include the tensor-correlated wave function F D |AMD explicitly so that the strong tensor interaction, exciting |AMD to the F D |AMD state, provides a large attractive contribution to the total energy. All the parameters in the wave function (1) are variational parameters. They are fixed by the energy minimization of the many-body Hamiltonian: We have to calculate all the necessary matrix elements for the two and three-body interactions using the TOAMD wave function. This form of the wave function (1) was studied for the 4 He nucleus by Nagata et al. to study the role of the tensor interaction [17]. The overlap integral of the AMD wave function is written as: The single-particle matrix elements are written as: The spatial matrix element is: As for the spin part, we use the following notation for the matrix element: For the isospin part, we use the following notation: where the states p and q are both the proton states or neutron states forM pq = 1, and they are different states forM pq = 0. Altogether the overlap matrix element is written as: We will give the overlap matrix element of the tensor-correlated wave function F D |AMD later, where we describe the method of handling multiple tensor correlation operators. In the next subsection, we introduce the Hamiltonian and later calculate the matrix element of the Hamiltonian for the AMD wave function AMD|H|AMD . 5/36

Nucleon-nucleon interaction and three-body interaction
We take the many-body Hamiltonian as a summation of the kinetic, two-body and threebody interactions: Here, T is the many-body kinetic energy: Here, we use natural units: = c = 1. The first term is the sum of the individual kinetic energies and T c.m. is the center-of-mass (c.m.) kinetic energy. Hence, T denotes the kinetic energy of the intrinsic motion of nucleons. Wiringa et al. constructed phenomenological two-body nucleon-nucleon (N N ) interactions without and with ∆(1232) degrees of freedom, the Argonne v14 (AV14) and v28 (AV28) models, respectively [2,3]. In principle we are able to treat any interactions, but here we write all the necessary ingredients to treat the AV14 potential. We briefly review the content of the N N interactions of the AV14 potential. We consider the ∆(1232) degrees of freedom in terms of a three-body interaction instead of treating ∆(1232) explicitly [18]. The two-body interaction AV14 is written as the sum of many operators: There are 14 operators in the AV14 potential. We write the two-body interaction as 1 2 i =j instead of i<j in order to treat the antisymmetrization systematically. The Argonne potential has three radial components: a long-range one-pion exchange part v π , and phenomenological intermediate-range and short-range core parts v I and v S : where operators O p represent operators of spin, isospin, tensor, spin-orbit, squared angular momentum, and squared spin-orbit interactions. In our calculations, we expand all the radial dependence in terms of the sum of Gaussian functions: As for the tensor operator, we multiply r 2 ij in the radial dependence. We have the three-body interaction in order to describe quantitatively the nuclear system: There are two components for the three-body interaction in the Urbana series of threenucleon potentials [1], where one term originates from two-pion exchange through delta excitation and the other term originates from the relativistic effect. We write 1 2 i =j =k instead of 6/36 cyclic i<j<k U (ijk) = i<j<k (U (ijk) + U (jki) + U (kij)) for systematic manipulations. The two-pion exchange term is written as: Here, all the operators with the radial dependences are those of the one-pion exchange: Here, ξ Y and ξ T represent the short-range cut-off factors of these interactions [1]. We expand Y (r) and T (r) in the sum of Gaussian functions. As for the relativistic effect term, the Argonne group chose the following expression: This three-body interaction does not depend on spin and isospin.

One-, two-, and three-body matrix elements
We describe here how to calculate one-, two-, and three-body matrix elements. We want to show all the necessary ingredients for matrix elements of the AMD wave function. We give the matrix elements of the kinetic energy and those for the central, tensor, spin-orbit, squared angular momentum, and squared spin-orbit interactions, and for the three-body interactions.

Kinetic energy
The matrix element of the individual kinetic energy for the AMD wave function is written as: Here, C(p : q) = (B −1 ) qp |B| is the determinant of a co-factor matrix for the pq element of the overlap matrix B and it is obtained using the properties of the determinant. All the details of the co-factor matrix C are given in Appendix A. Here, (B −1 ) qp denotes the qp component of the inverse matrix of the overlap matrix B. The single-particle matrix element is calculated in rectangular coordinates for the AMD single-particle states: Here, M pq is the spin matrix element (13) andM pq the isospin matrix element (14).

7/36
As for the center-of-mass kinetic energy, the center-of-mass wave function of the AMD wave function is: where R is the center-of-mass coordinate,ν = Aν, andD = 1 The center-of-mass kinetic energy is −1 2(Am) ∇ 2 R andD = 0 is assumed. Hence, the matrix element is: This is a single-particle kinetic energy. Since the interactions and correlations are written in terms of the relative (intrinsic) coordinates and momenta, the center-of-mass state is unaffected. Therefore, the intrinsic energy should be obtained after calculating everything using the individual kinetic energy and subtracting the center-of-mass kinetic energy (28).

Central interaction
The central interaction is expanded in the Gaussian functions: In order to perform the antisymmetrization in a systematic manner (Appendix A), we write the interaction in the momentum representation so that the operator has a separable form.
The momentum integration is written in a shorthand notation: We can write the matrix element as: By writing the interaction in the momentum representation, we are able to handle the antisymmetrization systematically. We give the antisymmetrization in terms of a co-factor matrix for the pq : st element of the B matrix: We can calculate the single-particle matrix element easily using the Gaussian integration formula and write the final expression as: Here, the spatial matrix element is: 8/36 All the coefficients A, B, C in the Gaussian integral (34) are given in Appendix B. Since we have to handle many Gaussian integrals, we have introduced a notation for the Gaussian integrals of I (type) (A, B, C). We introduce the superscript as (type) for the integral, where type denotes which nucleons interact with each other. In the present example (34), particles 1 and 2 interact with each other by a two-body interaction (type = 12). All the Gaussian integrals I (type) (A, B, C) up to three-body operators are given explicitly in Appendix B.
In the case of spin-spin interaction σ 1 · σ 2 , we have to replace the spin matrix elements We calculate these spin matrices using the Pauli spin operators and the spin-spinor wave functions. In the case of isospin-isospin interaction τ 1 · τ 2 , we have to replace the isospin matrix elementsM psM qt with xM ps xM qt x , wherē

Tensor interaction
The matrix element of the tensor interaction, where we see all the essential features of the tensor operator, is interesting. The tensor interaction in the momentum space is written as: The matrix element of the tensor interaction is written as: Using various matrix elements we can write it as: Here, the momentum integration is defined as: where the coefficients A, B, C are given in Appendix B. The subscript 1x1y is related to the momentum integrations for momenta k 1x k 1y outside of the exponent (38) for multiple momentum integrations to be discussed in the next section. If we have isospin-isospin interaction, we should replaceM psM qt with xM ps xM qt x . Hereafter, we do not write this 9/36 statement for other interactions. Later we will introduce a matrix form for treating multiple tensor operator cases and describe a systematic method to calculate momentum integrations.

Spin-orbit interaction
The spin-orbit interaction is written as: where the relative orbital angular momentum l ij is written as l ij = r ij × 1 2 ( p i − p j ). In the momentum space it is written as: We need the matrix element of p x = −i∇ x : The matrix element of the spin-orbit interaction is written as: Using various matrix elements we can write it as:

Squared angular momentum interaction
The squared angular momentum interaction is written as: 10/36 In the momentum space it is written as: We write the matrix element of the double-differential operator p x p y = −∇ x ∇ y : The matrix element of the squared angular momentum interaction using the single (43) and double derivatives (48) is written as:

Squared spin-orbit interaction
The squared spin-orbit interaction is written as: In the momentum space it is written as: 11/36 The matrix element of the squared spin-orbit interaction is written as:

Short-range three-body interaction
The three-body interaction U contains two components, one due to two-pion exchange exciting the ∆ state and one due to the relativistic effect. For calculation of the matrix element it is simpler to give the relativistic-effect one first. The short-range three-body interaction due to the relativistic effect is written as: Hence, the matrix element is Here, C(pqr : stu) is the determinant of the co-factor matrix for the pqr : stu matrix of the overlap matrix B as explained in Appendix A. The explicit form of the Gaussian integral I (12:23) is written in Appendix B, where (type = 12 : 23) indicates that the first two-body operator acts on particles 1 and 2 and the second two-body operator on particles 2 and 3.

Two-pion three-body interaction
There are two terms for the two-pion three-body interaction, First we decompose {X ij , X jk } in the spin and tensor operators: These spatial-spin operators are multiplied by the isospin operators. 12/36 We start with the spin-spin three-body part: Hence, the matrix element is: We write the tensor-tensor three-body part: Hence, the matrix element is: We write the spin-tensor three-body part: Here, because of the change in the particle coordinates, the third and fourth terms in Eq. (55) are identical and we multiply by 2 for the spin-tensor three-body term. Hence, the matrix 13/36 element is: As for the commutator terms: the additions for the spin part of the above matrix elements ((58), (60), (62)) are replaced by subtractions and the additions for the isospin part are also replaced by subtractions.

Two-body interaction with one tensor correlation operator
We give here the transition matrix element from the AMD wave function |AMD to the AMD wave function with tensor correlation F D |AMD . This matrix element includes the overlap matrix of the tensor-correlated state of the TOAMD state.
Here, V is given as a summation over the particle coordinates in Eq. (18) and F D is also given as a summation over the particle coordinates in Eq. (6). Hence, there are various many-body operators: Here, the symmetry factors S in front of each term in Q n such as 1 2 , 1, and 1 4 are obtained by taking into account the symmetry of interchange of the particle coordinates. S is tabulated in Table 1 of Sect. 6. Hence, there appear two-body, three-body and four-body operators, which are written as Q 2 , Q 3 , and Q 4 . We shall discuss the matrix elements of these operators one by one.

Two-body term
We write the case of the tensor interaction explicitly. As for the central interaction, the expressions are similar: Here, the first bracket corresponds to the tensor interaction and the second bracket the tensor correlation. Each tensor operator is expressed in the momentum space using the following 14/36 expression: From here we introduce a shorthand notation for the coefficient: The matrix element of the tensor interaction is written as: Since we can write the single-particle matrix element, in Gaussian form, we are left with integrations over momentum for the entire matrix element. We describe the systematic method of the Gaussian integrals in the next section and in Appendix B, we give the final expression for the tensor matrix element: Here, the function I ( (12) 2 ) 1x1y2z2u (A, B, C) is the result of the momentum integration, which will appear in the next section. The explicit form for A, B, C is given in Appendix B.

Three-body term
The Q 3 operator is written as: The Q 3 operator corresponds to the case of type = 12 : 23 for the momentum integration in the specification of Appendix B. Hence, we can write the matrix element as: 15/36

Four-body term
The Q 4 operator is written as: The momentum integration of this four-body operator should be that of type = 12 : 34 in Appendix B. Hence, we are able to write the final result as: We are able to write the momentum integral I (12:34) 1x1y2z2u (A, B, C) in a separable form and simplify the calculation of the matrix element. If we take an interaction other than the tensor interaction, we should change the above formula slightly depending on the type of operators.

Momentum integration and systematic differentiation
As we have seen in the calculation of matrix elements, we have to perform various types of momentum integrations (k 1 , . . . , k l ) with Gaussian functions in the following form: We have already used the cases where the number of the momentum integration is 1 in Eq. (39) and 2 in Eq. (71). Here, a vector k, a matrix A and other quantities are defined as: Here, A ij is a fraction of 1/ν and depends on the type of momentum integrals. and By construction, the matrix A is a real symmetric matrix and can be written as A = A 1/2 A 1/2 . Again, the subscripts as p, q · · · of the D p ′ s depend on the type of momentum integrals. Here, we have included the source terms whose coefficients are written as b 1 , b 2 , . . . , b l . 16/36 The source terms provide the term e i bi ki in the integrand. These source terms are used for the calculation of integrals of the form: Here, i, j, . . . , k (≤ l) represent the momenta and b stands for the set of b i . The momentum directions are denoted by x, y, . . . , z. We shall obtain the final results by setting the coefficients of the source terms b to zero and write: These results of the momentum integrations are written in terms of the AMD wave functions with ν and various D p ′ s and the interaction ranges a µ .
We calculate I ixiy···kz one by one. In the following, n denotes the number of momenta in the integration multiplied by the Gaussian functions. n = 0 We first take the integration of the basic integral: This multiple Gaussian integration is verified for a symmetric matrix A with the existence of a square root matrix A 1/2 . The front factor 1 (2π) 3l comes from the definition of the momentum integration in Eq. (31). det |A| in the denominator appears from the Jacobian of the change of the integration variables and det |A −1/2 | = 1 √ det |A| is used.
Here, b ix is included in B and the derivative is: where Here, a very interesting relation is: 17/36 where 1 appears at the i th position in the above vector of dimension l. n = 2 where Here, we have a symmetry in that the results are unchanged by changing the order of ix and jy: Here, we have used the fact that the derivative of D ij δ xy is zero: Since the derivative is done successively, we find several interesting rules found by deriving these expressions, which are useful for derivation of higher-order derivative terms: • The derivative terms are written in terms of only D ij δ xy and iE ix . Because of this fact, we simply write D ij δ xy → D αβ and iE ix → E α to express the derivative formula. • These α and β denote ix etc. At the same time, they can also mean the successive order of derivatives: α < β < · · · .
It is then interesting to write the properties of the differentiation using the numbering notation I αβ··· (A, B, C : b) using the above results up to n = 3. I αβ··· consists of the sum of D k αβ (E n−2k γ ) terms with k = 1 · · · [ n 2 ]. I αβ··· is symmetric for any exchange of the order of α, β, . . . due to the interchangeable property of the differentiation. Observing the terms with subscript α and the procedure of the above manipulations, the first subscript α appears once in all the terms keeping its position at the beginning for each term. 18/36 Hence, we can formally write the integrals as: where [ n 2 ] is n/2 for even n and (n − 1)/2 for odd n. Q k q denotes configurations of all the derivative terms αβγ · · · to appear in the D k and E n−2k terms. N n k is the number of all the terms, where all αβγ · · · are partitioned in the D k and E n−2k terms.
In order to write possible partitions of αβγ · · · , we consider the term D k E n−2k and write rules for the configurations. We write the configurations as: Here, αβγ · · · are written in the form of a i b j c k · · · . Given the observations above for the exchange property of αβ · · · and the order of the appearance of α, we can write the following rules: • Rule 1: With these rules, we can write all the configurations without double counting of them. It is interesting to calculate the number of terms for each partition N n k . There are altogether n! ways to order the derivatives αβγ · · · . There are three rules to avoid double counting of partitions in the n! ways. Now, the order of c i in E n−2k is fixed to a unique one from rule 3. There are (n − 2k)! ways to order the derivatives, but the order is fixed to one by the rule 3. Hence, we have to divide by (n − 2k)! out of the entire possibility n!. We can have k pairs, we have to divide by 2 k from rule 1. In addition, there are k D ′ s; we have to order these D ′ s using rule 2 and we have to divide by k!. Hence the number of configurations for each term is: For n = 3 and k = 1 it is N 3 1 = 3, and for k = 0 it is N 3 0 = 1. These numbers of terms agree with the above results.
Using the above rules, we write the n = 4 and n = 5 cases explicitly, keeping in mind the orders of α · · · γ in each configuration.  19/36 This form is written following the three rules above. Using the formula, we can calculate the number of configurations: N 4 2 = 3, N 4 1 = 6, N 4 0 = 1. These numbers agree with the expression of I 1234 (A, B, C : b) in Eq. (96). n = 5 Using the formula, we get: N 5 2 = 15, N 5 1 = 10, N 5 0 = 1. These numbers agree with the numbers of each configuration in the above expression. We have derived the above expressions in various ways. Rules 1, 2, and 3 for the construction of all the partitions have been verified in the mathematical induction method. With the above rules we are able to write explicitly the derivative formula for any number of momentum integrations and derivatives.
After getting all the terms, we should set b = 0 and write the integrals in the same notation: In addition, we shall introduce the superscript (type) for the integrals to specify which nucleons interact with each other; that is discussed in Appendix B.

Two-body interaction with two tensor correlation operators
We write explicitly the case of the tensor interaction with two tensor correlations: where Here, R 2 , . . . , R 6 are two-to six-body operators. We have written the multiple of two twobody interactions (correlations) as the sum of two-, three-and four-body operators: in Eq. (65), where We calculate the multiple of three two-body interactions (correlations) R 2 , . . . , R 6 baesd on the multiple of two two-body interactions (correlations) Q 2 , Q 3 , and Q 4 . We start with the 20/36 two-body operator R 2 basing on Q 2 : We obtain the three-body operator R 3 based on Q 2 : We have the three-body operator R 3 based on Q 3 : We obtain the four-body operator R 4 based on Q 2 : We obtain the four-body operator R 4 based on Q 3 : We obtain the four-body operator R 4 based on Q 4 : We obtain the five-body operator R 5 based on Q 3 : We obtain the five-body operator R 5 based on Q 4 : We obtain the six-body operator R 6 based on Q 4 : The symmetry factor in front of the summation is tabulated in Table 1 of Sect. 6. 21/36

Matrix element of a multiple of three operators
We start with a two-body operator of the category type = (12) 3 : The momentum integral I We calculate the three-body operator R 3 (Q 2 ) of type = 12 : (23) 2 : (2) µ1C (2) µ2C We come to the first term of the three-body operator R 3 (Q 3 ) of the category type = (12) 2 : 23: The other two terms for R 3 (Q 3 ) are categorized as type = 23 : 12 : 23 and 13 : 12 : 23. The superscripts of the integral I have to be changed according to the type and the spin and isospin matrix elements should be changed slightly for these terms. We write the four-body operator R 4 (Q 2 ) of the category type = 12 : (34) 2 : There are many more matrix elements, which are obtained in the same way. We skip writing these matrix elements here. Other interactions are written in a similar way by changing the 22/36 characters of the operators, as discussed in Sect. 2. The systematic way to calculate matrix elements will be described in Sect. 6.

Short-range correlation
We have to introduce further the short-range correlation in the many-body wave function. The discussion of the short-range correlation has been delayed up to this section, since the main difficulty with the nuclear many-body problem is the treatment of the tensor correlation. We have developed various methods to handle the tensor correlation operator. We shall use the same concept as the tensor correlation for the short-range correlation. Several methods for the short-range correlation have been developed in the past. One popular method is the Jastrow correlation operator method [19] and another is the unitary correlation operator method (UCOM) [20]. In the Jastrow method, the correlation operator is written as a product of correlation functions. In the UCOM, a Hermite correlation operator is placed on an exponential so that the correlation operator is unitary. In all these methods, the matrix elements are obtained by introducing an approximation to take the resulting operators up to few-body operators. As discussed for the case of the tensor correlation operator, we are able to add correlation operators systematically one after another to see the convergence of the solutions. Hence, it is important to know that the present formulation is able to calculate all the matrix elements systematically and straightforwardly.

Wave function with short-range correlation
We introduce the short-range correlation operator F S in the TOAMD wave function as in the case of the tensor correlation operator.
where |Ψ was introduced as the TOAMD wave function in Eq. (1). This arrangement indicates that the full wave function is the sum of the following four components: The first term provides low-momentum components representing the shape of nucleus, the second term provides high-momentum components due to the short-range correlation, and the third term provides intermediate-hight-momentum components due to the tensor correlation. The last term is an interference term for the short-range and tensor correlations.
Here, we expand the short-range correlation operator in the sum of Gaussian functions: These expansion parameters are considered as variational parameters of the many-body wave function. The short-range correlation is strong in the non-spin, non-isospin channel and we show only this case. However, we will use spin-and isospin-dependent short-range correlations in the calculation. We can then obtain the many-body Schrödinger equation H|Ψ = E|Ψ and the eigenvalue is: Although the operator structure of the short-range correlation is simple, there appear many-body operators. Hence, as an example we discuss the case of the three-body interaction 23/36 with the short-range and tensor correlations in the next subsection. We then discuss a general method to calculate all the necessary matrix elements in the subsequent subsection.

Three-body interaction with the short-range correlation for the tensor component
We explicitly write a complicated matrix element where the number of momentum integrals is 6 : 2 from the short-range correlation, 2 from the tensor correlation, and 2 from the three-body interaction. We consider the case of the repulsive three-body interaction with the short-range correlation: Here, we have introduced an approximation that the short-range correlations act only on the same nucleon pairs, i and j in the three-body interaction. This is because the shortrange correlation F S is large only at very short distances and the probability of more than three nucleons coming to the region of the short-range correlation is negligibly small. If we multiply the two tensor correlation operators by the three-body operator (122) of the three-body interaction with two short range correlations F S U R F S , we get many terms: Here, one of the R 3 terms for the type = (12) 3 : 23 : (12) 2 , where three two-body operators act on particles 1 and 2, one two-body operator acts on particles 2 and 3, and two two-body operators act on particles 1 and 2: (2) µ6 k1 k2 k3 k4 k5 k6 e −k 2 1 /4aµ1 e −k 2 2 /4aµ2 e −k 2 3 /4aµ3 e −k 2 4 /4aµ4 e −k 2 5 /4aµ5 e −k 2 6 /4aµ6 e i k1( ri− rj) e i k2( ri− rj) e i k3( ri− rj) e i k4( rj− rk) e −i k5( ri− rj) e −i k6( ri− rj) 24/36 Here we write the case, where two tensor operators and two short-range operators act on i, j pairs and the three-body interaction works for i, j, k nucleons: We can further write other terms for the three-body operators in a similar way to that above. We have in addition the four-, five-, six-, and seven-body operators. All these matrix elements are written in a similar way to the above expression, with major changes in the momentum integrals and small changes in the spin matrix elements. Depending on the number of nucleons involved for the multi-body operators, the coefficients of the co-factor matrix change as C(pq · · · r : p ′ q ′ · · · r ′ ).
Here, S is a symmetry factor for a many-body operator. We list S for various configurations up to n = 3 in Table 1. The symmetry factors are the same for configurations obtained by the interchange of particle numbers. pq · · · rp ′ q ′ · · · r ′ are the quantum numbers of AMD states, where each quantum number p ′ s runs from 1 to A. µ1µ2 · · · µn are the expansion parameters 25/36 of interactions and correlations in terms of Gaussian functions and the number of µ ′ s is n representing the number of momenta. The expansion factorsC (m) µ are written as: where m = 1 for the spin-orbit interaction, m = 2 for the tensor, squared angular momentum and squared spin-orbit interactions, and m = 0 otherwise. The coordinates x, y, . . . , z run from x to z of the rectangular coordinates. The momentum integral I (type) X(x,y,...,z) (A, B, C) depends on the type of configurations for the coefficients A, B, C with differentiation of momentum given by a function X(x, y, . . . , z). F (x, y, . . . , z) represents the type of interactions and correlations to be specified by the interaction. M pq···rp ′ q ′ ···r ′ Z(x,y,...,z) is the spin matrix element of all the AMD states andM pq···rp ′ q ′ ···r ′ U (x,y,...,z) the isospin matrix element. The co-factor C(pq · · · r : p ′ q ′ · · · r ′ ) takes care of the antisymmetrization of particles in the many-body operators.
We describe the procedure of writing matrix elements for various operators. We first fix which matrix elements to calculate the configuration: type. Once we fix type, we can get the symmetry factor S for this configuration. We order all the operators from left to right and assign momenta k 1 , . . . , k n . We then write all the operators explicitly keeping the order of the operators. For each operator for particles i and j, we write the following factors:

Central interaction
We write the case of the spin-spin and isospin-isospin interaction: with m = 0. If there are no spin-and/or no isospin-dependent operators, we simply drop these spin and isospin operators.

Squared angular momentum interaction
26/36 with m = 2. For calculation of the matrix elements, we need matrix elements of p x = −i∇ x and p x p y = (−i∇ x )(−i∇ y ) for the calculation of angular momenta, which are given in Eqs.

Squared spin-orbit interaction
with m = 2. We can then write the matrix element explicitly. As an example, we show the case of type = 12 : (34) 3 : 56 for the tensor operator: The momentum integral is given in Appendix B. All matrix elements can be written in a similar systematic manner. As for the spin-orbit interaction with the short-range and tensor correlations, we give explicitly the case of type = 12 : (34) 3 : 56.
The spin-orbit interaction has a derivative term. Therefore the expressions are slightly different for different configurations. We write one similar configuration case: type = 12 : 27/36 (34) 3 : 35.

Summary
We have developed a powerful many-body theory to describe finite nuclei, which is calld "tensor-optimized antisymmetrized molecular dynamics" (TOAMD). The TOAMD theory is based on AMD, in which the concept of the TOSM is incorporated, in order to treat the strong tensor interaction in the nucleon-nucleon interaction. The tensor interaction is treated by the tensor correlation operator acting on the AMD wave function. Since the tensor interaction is of long and intermediate range, we have to explicitly treat many-body operators due to the tensor correlation operators and the two-and three-body interactions in the Hamiltonian for the AMD state. For the TOAMD theory, we have to treat up to 6-body operator terms for the two-body interaction and 7-body operator terms for the three-body interaction for the AMD state.
In order to treat multi-body operators for many-body nuclear systems efficiently, we should use all the powerful mathematics. The AMD wave function consists of shifted Gaussian functions with spin and isospin wave functions and all the interactions are expanded in Gaussian functions. In addition, we have to take into account the antisymmetrization of all the nucleons. For this purpose, we take the Fourier transforms of all the interactions so that any multi-body operators are written in separable forms in particle coordinates. We are then able to calculate multi-body matrix elements with antisymmetrization using the multi body co-factor matrix of the norm matrix. We should then take the multi-momentum integrations, where we have developed a systematic integral and differentiation method with source terms. The final results are written as the sum of many terms for the norm and energy matrices, which do not involve any numerical integrations. We have to minimize the total energy with respect to the variational parameters, which are the shift coordinates and spin weights in the AMD wave function, the tensor correlation function, and additionally, the short-range correlation function.
The TOAMD theory provides the total energy as function of the variational parameters, where the total energy can be calculated systematically in a straightforward manner. Since the matrix elements can be calculated in systematic methods, we are able to improve the calculated results by adding more and more complicated correlations including both the tensor and short-range correlations. For heavy nuclei, we have to perform the angular momentum projection and take the sum of the Slater determinants for better description 28/36 of many-body systems. The TOAMD theory is a powerful and economical method to treat nuclear many-body system. The formulation is transparent and we are able to calculate twoand three-body interactions with any order of tensor and short-range correlations. In the TOAMD theory, we should be able to calculate nuclei with many nucleons using the present powerful computers.
In the very near future, we shall publish numerical results of s-shell nuclei, and successively, p-shell and heavier nuclei in the TOAMD theory using the bare nucleon-nucleon interaction.
Here, C(r 1 r 2 : l 1 l 2 ) is the determinant of a co-factor matrix of B, where the r 1 and r 2 rows and l 1 and l 2 columns are removed from the A × A matrix. Again C(r 1 r 2 : l 1 l 2 ) is a function of the single particle overlap p i |q j and does not depend on the one-body matrix element.

Determinant of the co-factor matrix
We derive here the determinant C(r 1 r 2 · · · r k : l 1 l 2 · · · l k ) of the co-factor matrix of the A × A matrix B. We have the following identity for the A × A matrix B with the matrix elements a ij : Here, B −1 is the inverse matrix of B. Since this is an important formula to get the determinant of the co-factor matrix, we verify this explicitly using the definition of the determinant: Hence, we can write the determinant of B as: a k1l1 a k2l2 · · · a krlr C(k 1 · · · k r : l 1 · · · l r ) .
In the above formula, k 1 , k 2 , . . . , k r are any row numbers of the original matrix B. Here, C(k 1 · · · k r : l 1 · · · l r ) is the determinant of a co-factor matrix of B, where k 1 , k 2 ,. . . and k r 30/36 rows and l 1 , l 2 ,. . . and l r columns are removed from the matrix B. The coefficient C(k 1 · · · k r : l 1 · · · l r ) is written using a ij of the original matrix B. Explicitly, the coefficient C is given as C(k 1 · · · k r : l 1 · · · l r ) = Compared with the matrices of one-, two-and multi-body operators ((A1), (A2), and (A3)), the coefficients C(k 1 · · · k r : l 1 · · · l r ) are the same as that given in Eq. (A7). In the case r = A, the co-factor C becomes ǫ(P (k 1 · · · k r : l 1 · · · l r )), which is the phase of indicated permutation. Hence, for r = A, the rhs of Eq. (A3) is simply the determinant of the full matrix of p i |O|q j .

B. Gaussian integrals
We give here all the necessary Gaussian integrals: · · · kn f p |e iK1(k1,··· ,kn)r |g p ′ · · · f q |e iKm(k1,··· ,kn)r |g q ′ = 1 (2π) 3n Here, type denotes which operators act between which particles. The functions K i are functions of k i , whose explicit forms depend on the type of multi-body operators. In order to understand the meaning of A, B, C, we write the single-particle matrix element: We write here all the possible integrals up to the three two-body operators. There are still several Gaussian integrals necessary for calculations of matrix elements, but they can be obtained in a similar way to those presented in this appendix.

type = 12
In this appendix, we write vectors ( D p ) ′ s simply without the vector notation as D p . We write here type = 12, which indicates that a two-body operator act on particles 1 and 2.
Here type = (12) 2 means that two two-body operators act on particles 1 and 2. type = 12 : 23 Here, type = 12 : 23 means that one interaction acts on particles 1 and 2 and the other on particles 2 and 3. We omit the explanation of type in the following. type = 12 : 34 We give a more complicated case as an example so that the general rule for A, B, C can be understood instead of writing all the cases. In this example, type = 12 : (34) 3 : 56 indicates that one two-body operator acts on particles 1 and 2, three two-body operators on particles 3 and 4, and one two-body operator on particles 5 and 6.
The diagonal terms in A are fixed simply as in all the other cases, while the non-diagonal terms depend on the type of configuration. They can be obtained by calculating the k 2 35/36 terms: k 2 1 + k 2 1 + (k 2 + k 3 + k 4 ) 2 + (k 2 + k 3 + k 4 ) 2 + k 2 5 + k 2 5 = 2(k 2 1 + k 2 2 + k 2 3 + k 2 4 + k 2 5 ) + 4(k 2 k 3 + k 2 k 4 + k 3 k 4 ) . (B20) For B, the k 1 term is taken by the pp ′ : qq ′ state, the k 2 , k 3 , k 4 terms are taken by the rr ′ : ss ′ state and the k 5 term by the tt ′ : uu ′ state. For C, it is simply written by using all the 12 states. We can give the Gaussian integrals systematically, as shown here, but for the case in which the interactions and correlations are separable, the Gaussian integrals can be written as a product of various Gaussian integrals. For an example of this case, we can write: