Density-dependent nucleon-nucleon interaction from three-nucleon forces

Density-dependent nucleon-nucleon interaction from three-nucleon forces

Alessandro Lovato    Omar Benhar    Stefano Fantoni    Alexey Yu. Illarionov    Kevin E. Schmidt SISSA and INFN, Sezione di Trieste. I-34014 Trieste, Italy
INFN, Sezione di Roma. I-00185 Roma, Italy
Dipartimento di Fisica, Università “La Sapienza”. I-00185 Roma, Italy
CNR-DEMOCRITOS National Supercomputing Center. I-34014 Trieste, Italy
Dipartimento di Fisica, Università di Trento. I-38123 Povo, Trento, Italy
Department of Physics, Arizona State University, Tempe, AZ 85287
July 12, 2019
Abstract

Microscopic calculations based on realistic nuclear hamiltonians, while yielding accurate results for the energies of the ground and low-lying excited states of nuclei with , fail to reproduce the empirical equilibrium properties of nuclear matter, that are known to be significantly affected by three-nucleon forces. We discuss a scheme suitable to construct a density-dependent two-nucleon potential, in which the effects of -particle interactions can be included by integrating out the degrees of freedom of -nucleons. Our approach, based on the formalism of correlated basis function and state-of-the-art models of the two- and three-nucleon potentials, leads to an effective interaction that can be easily employed in nuclear matter calculations, yielding results in good agreement with those obtained from the underlying three-body potential.

pacs:
21.30.Fe, 21.45.Ff, 21.65.-f

I Introduction

The results of ab initio microscopic calculations consistently suggest that realistic nuclear hamiltonians, including both two- and three-nucleon potentials, while providing a quantitative account of the energies of the ground and low-lying excited states of nuclei with steve_bob (); steve_RNC (), fail to explain the empirical equilibrium properties of nuclear matter. This problem can be largely ascribed to the uncertainties associated with the description of three-nucleon interactions, whose contribution turns out to be significant.

A signal of the limitations of the commonly employed three-nucleon potential models (e.g. the Urbana IX model of Ref. UIX ()) has been recently provided by the authors of Ref. AFDMC1 (), who carried out a study of symmetric nuclear matter within the Auxiliary Field Diffusion Monte Carlo (AFDMC) approach. Their results, obtained using a truncated version V8P () of the state-of-the-art nucleon-nucleon potential of Ref. av18 (), show that AFDMC simulations do not lead to an increase of the binding energy predicted by Fermi-Hyper-Netted-Chain (FHNC) and Brueckner-Hartree-Fock (BHF) calculations FHNC-BHF ().

Different three-nucleon potential models UIX (); ILL (), yielding similar results when applied to the calculation of nuclear properties, predict sizably different equations of state (EoS) of pure neutron matter at zero temperature and densities exceeding the nuclear matter saturation density, fm AFDMC2 (). In this region, the three-nucleon force contribution to the binding energy becomes very large, the ratio between the potential energies associated with two- and the three-body interactions being at density (see, e.g. Ref. akmal ()). The size of the three-body force contribution suggests that, at large , interactions involving four or more nucleons may also play an important role, and should be taken into account.

In view of the severe difficulties involved in the implementation of the existing models of three-nucleon interactions in many-body calculations, the explicit inclusion of four- and more-body potentials does not appear to be a viable option. In this paper we follow a different strategy, somewhat along the line of the Three-Nucleon-Interaction (TNI) model proposed by Lagaris and Pandharipande LP1 () and Friedman and Pandharipande FrP () in the 1980s.

The authors of Refs. LP1 (); FrP () suggested that the main effects of three- and many-nucleon forces can be taken into account through an effective, density-dependent two-nucleon potential. However, they adopted a purely phenomenological procedure, lacking a clearcut interpretation based on the the analysis of many-nucleon interactions at microscopic level.

The TNI potential consists of two density-dependent functions involving three free parameters, whose values were determined through a fit of the saturation density, binding energy per nucleon and compressibility of symmetric nuclear matter (SNM), obtained from FHNC variational calculations. The numerical values of the three model parameters resulting from recent calculations performed by using AFDMC simulations turn out to be only marginally different from those of the original TNI potential Gand2010 ().

The TNI potential has been successfully applied to obtain a variety of nuclear matter properties, such as the nucleon momentum distribution mom_dis (), the linear response response_1 (); response_2 (), and the Green’s function spec_func1 (); spec_func2 (), within the Correlated Basis Function (CBF) approach (for a review of CBF theory and its applications, see Ref.fabrocini_fantoni () and references therein).

The strategy based on the development of two–body density-dependent potentials has been later abandoned, because their application to the study of finite nuclei involves a great deal of complication, mainly stemming from the breakdown of translation invariance. While in uniform matter the density is constant and the expansion of the effective potential in powers of is straightforward, in nuclei different powers of the density correspond to different operators, whose treatment is highly non trivial.

However, the recent developments in numerical methods for light nuclei seem to indicate that the above difficulties may turn out to be much less severe then those implied in the modeling of explicit many–body forces and, even more, in their use in ab initio nuclear calculations.

In view of the observation, based on a variety of experimental evidence RMP1 (); Pand1997 (), that short range nucleon–nucleon (NN) correlations are a fundamental feature of nuclear structure, the description of nuclear dynamics in terms of interactions derived in coordinate space, like the Urbana-Argonne models, appears to be the most appropriate, for both conceptual and technical reasons.

First of all, correlations between nucleons are predominantly of spatial nature, in analogy with what happens in all known strongly correlated systems. In addition, one needs to clearly distinguish the effects due to the short–range repulsion from those due to relativity. Finally, quantum Monte Carlo methods have serious difficulties in dealing with highly non local interactions. For all the above reasons we stick to two–body density-dependent potentials of the Urbana-Argonne type.

Our approach is based on the tenet that -body potentials () can be replaced by an effective two-nucleon potential, obtained through an average over the degrees of freedom of particles. Hence, the effective potential can be written as a sum of contributions ordered according to powers of density, the -th order term being associated with -nucleon forces.

Obviously, such an approach requires that the average be carried out using a formalism suitable to account for the full complexity of nuclear dynamics. Our results show that, in doing such reduction, of great importance is the proper inclusion of both dynamical and statistical NN correlations, whose effects on many nuclear observables have been found to be large RMP1 (); Pand1997 ().

In this work, we use CBF and the Fantoni-Rosati (FR) cluster expansion formalism FR_cluster (); fabrocini_fantoni (); ariasdesaavedra_bisconti_co_fabrocini () to perform the calculation of the terms linear in density of the effective potential, arising from the irreducible three-nucleon interactions modeled by the UIX potential.

It should be noticed that our approach significantly improves on the TNI model, as the resulting potential is obtained from a realistic microscopic three-nucleon force, which provides an accurate description of the properties of light nuclei.

While being the first step on a long road, the results discussed in this paper are valuable in their own right, as the effective potential can be easily implemented in the AFDMC computational scheme to obtain the EoS of SNM. Similar calculations using the UIX potential are not yet possible, due to the complexities arising from the commutator term. In addition, the density-dependent potential can be used to include the effects of three-nucleon interactions in the calculation of the nucleon-nucleon scattering cross section in the nuclear medium. The knowledge of this quantity is required to obtain a number of nuclear matter properties of astrophysical interest, ranging from the transport coefficients to the neutrino emission rates BV (); BFFV ().

In Section II we discuss the main features of the existing theoretical models of the three-nucleon force, while Section III is devoted to a brief review of the many-body approach based on CBF and the cluster expansion technique. In Section IV we describe the derivation of the density-dependent interaction, pointing out the role of dynamical and statistical correlation effects. In Section V we compare the energy per particle of nuclear matter obtained from the effective potential to that resulting from highly refined calculations, carried out using the Argonne V6P () and V8P () nucleon-nucleon potentials and the Urbana IX three-nucleon potential UIX (). Finally, in Section VI we summarize our findings and state the conclusions.

Ii Three nucleon forces

Nuclear many-body theory (NMBT) is based on the assumption that nuclei can be described in terms of point like nucleons of mass , whose dynamics are dictated by the hamiltonian

(1)

Before describing the three nucleon potential , let us discuss the two-body potential model, that will be used throughout the paper. It is given by

(2)

where

(3)

In the above equation, and , where and are Pauli matrices acting on the spin or isospin of the -th, while

(4)

with , is the tensor operator, is the relative angular momentum

(5)

and is the total spin of the pair

(6)

Such potentials have exactly the same form as the first eight components of the state-of-the-art Argonne potential av18 (). We will be using the so called Argonne and Argonne potentials, which are not simple truncations of the Argonne potential, but rather reprojections V6P ().

The Argonne potential is obtained by refitting the scattering data in such a way that all and partial waves as well as the wave and its coupling to are reproduced equally well as in Argonne . In all light nuclei and nuclear matter calculations the results obtained with the are very close to those obtained with the full , and the difference can be safely treated perturbatively.

The Argonne is not just a truncation of , as the radial functions associated with the first six operators are adjusted to preserve the deuteron binding energy. Our interest in this potential is mostly due to the fact that AFDMC simulations of nuclei and nuclear matter can be performed most accurately with –type of two–body interactions. Work to include the spin–orbit terms in AFDMC calculations is in progress. On the other hand we need to check the accuracy of our proposed density-dependent reduction with both FHNC and AFDMC many–body methods before proceeding to the construction of a realistic two–body density-dependent model potential and comparing with experimental data.

It is well known that using a nuclear Hamiltonian including only two-nucleon interactions leads to the underbinding of light nuclei and overestimating the equilibrium density of nuclear matter. Hence, the contribution of three-nucleon interactions must necessarily be taken into account, by adding to the Hamiltonian the corresponding potential, e.g. the widely used Urbana IX (UIX) UIX ().

The potential of Ref. UIX () consists of two terms. The attractive two-pion () exchange interaction turns out to be helpful in fixing the problem of light nuclei, but makes the nuclear matter energy worse. The purely phenomenological repulsive term prevents nuclear matter from being overbound at large density.

Figure 1: Feynman Diagram associated with the Fujita Miyazawa three-nucleon potential term.

The term was first introduced by Fujita and Miyazawa Fujita_Miyazawa () to describe the process whereby two pions are exchanged among nucleons and a resonance is excited in the intermediate state, as shown in the Feynman diagram of Fig. 1. It can be conveniently written in the form

(7)

where

(8)

and

(9)

The radial functions associated with the spin and tensor components read

(10)
(11)

while the are short-range cutoff functions defined by

(12)

In the UIX model, the cutoff parameter is kept fixed at fm, the same value as in the cutoff functions appearing in the one-pion exchange term of the Argonne two-body potential. On the other hand, is varied to fit the observed binding energies of H and He. The three-nucleon interaction depends on the choice of the NN potential; for example, using the Argonne model one gets .

The repulsive term is spin-isospin independent and can be written in the simple form

(13)

with defined in Eq. (11). The strength , adjusted to reproduce the empirical nuclear matter saturation density, is with .

The two parameters and have different values for and . We disregard such differences in the present analysis, mostly aimed at testing the quality of our density-dependent reduction of the UIX three–body potential, rather than reproducing empirical data.

Iii formalism

iii.1 Correlated basis theories

One of the most prominent features of the nucleon-nucleon (NN) interaction is the presence of a repulsive core, giving rise to strong correlations that cannot be taken into account within the independent particle picture.

This problem has long been recognized, and was clearly pointed out by Blatt and Weisskopf over fifty years ago. In their classic Nuclear Physics book, first published in 1952, they warn the reader that “the limitation of any independent particle model lies in its inability to encompass the correlation between the positions and spins of the various particles in the system” BW ().

Let us consider uniform nuclear matter, defined as a translationally invariant system of protons and neutrons, in which the electromagnetic interaction is turned off. In the absence of interactions, such a system can be described as a Fermi gas at zero temperature, and its ground state wave function reduces to the antisymmetrized product (Slater determinant) of orbitals associated with the single particle states belonging to the Fermi sea:

(14)

with

(15)

and

(16)

In the above equations, is the normalization volume, and are Pauli spinors, describing the nucleon spin and isospin and . Here is the Fermi momentum, while and denote the density and the degeneracy of the momentum eigenstates [2, 4 pure neutron matter (PNM) and symmetric nuclear matter (SNM), respectively]. The generalized coordinate represents both the position and the spin-isospin variables of the -th nucleon, while denotes the set of quantum numbers specifying the single particle state.

The antisymmetrization operator can be written in the form

(17)

where

(18)

is the two-particle exchange operator, defined by the relation

(19)

Note that, as shown by Eq.(18), the exchange operators act on both the radial and spin-isospin components of the nucleon wave function.

Due to the strong repulsive core, the matrix elements of between eigenstates of the non interacting system

(20)

turn out to be very large, or even divergent if the core of the NN potential is infinite. As a consequence, perturbative calculations carried out using the bare NN potential and the Fermi gas basis states are unavoidably plagued by lack of convergence.

To circumvent this problem, one can follow two different strategies, leading to either G-matrix or CBF perturbation theory. Within the former approach, the bare potential is replaced by a well behaved effective interactions, the so called G-matrix, which is obtained by summing up the series of particle–particle ladder diagrams. In the second approach, nonperturbative effects are handled through a change of basis functions.

Correlated basis theories of Fermi liquidsfabrocini_fantoni (); ariasdesaavedra_bisconti_co_fabrocini (); clark (); OCBfantoni () are a natural extension of variational approaches in which the trial ground state wave function is written in the form

(21)

In the above equation, is a correlation operator, whose structure reflects the complexity of the nucleon-nucleon potential wiringa_pandha_1 ():

(22)

with

(23)

Note that the symmetrization operator is needed to fulfill the requirement of antisymmetrization of the state , since, in general, . The correlated basis (CB) is defined as

(24)

where is a generic n particle – n hole Fermi gas state. The CB states are normalized but not orthogonal to each other. They have been used within non orthogonal perturbation theory clark (); krotscheck () to study various properties of quantum liquids. An exhaustive analysis of the convergence properties of this perturbation scheme has never been carried out, but the truncation of the series at a given perturbative order is known to lead to nonorthogonality spuriosities, whose effects are not always negligible. A much safer and efficient procedure, in which one first orthogonalizes the CB states by using a combination of Schmidt and Lwdin transformations and then uses standard perturbation theory, has been proposed by Fantoni and Pandharipande OCBfantoni ().

The radial functions , appearing in the definition of the correlation operator are determined by the minimization of the energy expectation value

(25)

which provides an upper bound to the true ground state energy . In principle, that can be done by solving the Euler equations resulting from the functional minimization of with respect to the correlation functions , in analogy with what has been done in Jastrow theory of liquid He. However, the presence of the spin-isospin dependent correlation operators and their non–commutativity makes the application of this procedure to nucleonic systems almost prohibitive. In this case the functional minimization can be carried in a more straightforward fashion on the lowest order cluster contribution to , paying the price of introducing proper constraints and the associated variational parameters, as discussed below.

iii.2 Cluster expansion

In CBF theories the calculation of is carried out by i) expanding the r.h.s. of eq. (25) in powers of proper expanding functions that vanish in uncorrelated matter and ii) summing up the main series of the resulting cluster terms by solving a set of coupled integral equations. The FR cluster expansion FR_cluster () has been derived to accomplish the first of these two steps for the case of Jastrow correlated models. It has been obtained through a generalization of the concepts underlying the Mayer expansion scheme, originally developed to describe classical liquids mayer (), to the case of quantum Bose and Fermi systems. In this case the expanding quantity is given by

(26)

where is the only correlation of the Jastrow model. Notice that in the calculation of the kinetic energy operators also act on the correlation functions, giving rise to additional expanding quantities. For the sake of simplicity, and since here we are interested in calculating the expectation values of two- and three-body potentials, we will not address this issue. It has been proved FR_cluster () that the FR cluster expansion is linked, and therefore does not suffer the appearance of infinities in the thermodynamic limit. The FR techniques have been subsequently extended and extensively used to deal with spin–isospin dependent correlation operators, like those of Eq. (22) wiringa_pandha_1 (); wiringa_pandha_2 (). In this case, besides the expanding function of eq. (26) one has to also consider the following ones:

(27)

The cluster terms are most conveniently represented by diagrams wiringa_pandha_1 (); fabrocini_fantoni (). The diagrammatic representation of the above expanding functions is given by the bonds displayed in Fig. 2, in which is represented by a dashed line, by a single wavy line, and by a doubly wavy line .

Figure 2: Different kinds of correlation bonds.

The is the largest of all the expanding quantities and the cluster terms (and similarly the corresponding cluster diagrams) involving these functions (or bonds) have to be summed up as massively as possible. This can be accomplished by solving the FHNC equations of Ref. fantoni_rosati (), which already in their basic form sum up all diagrams at all orders, with the exception of the so called elementary diagrams.

The cluster diagrams involving operatorial bonds, like those representing the functions given in Eq. (27), cannot be summed up as massively as the scalar diagrams of FHNC type. This is due to the additional complexity associated with the non commutativity of the spin–isospin dependent correlation operators. The most powerful summation scheme which has been derived so far is the so called Single Operator Chain of Refs. (wiringa_pandha_1 (); wiringa_pandha_2 ()), generally denoted as FHNC/SOC approximation.

For the sake of clarity, in the following we summarize the main features of the FR cluster expansion and the FHNC/SOC approximation, extensively reviewed in wiringa_pandha_1 (); fabrocini_fantoni ().

Let us consider the expectation value of the NN potential. Exploiting the symmetry properties of the wave function it can be written in the form

(28)

Numerator and denominator of the above equation are expanded in powers of the functions defined above. The expansion of , with for the numerator and the denominator, leads to series of terms, say , where the labels N and D stand for numerator and denominator, respectively, each characterized by the number of correlated nucleons, i.e. those appearing in the argument of the expanding functions. Integrating such terms over the variables of the remaining uncorrelated nucleons amounts to multiplying by the –body Fermi distribution operator . Consider for instance one of the cluster terms of the numerator, whose structure is given by

(29)

where and stands for integration over the coordinate and tracing over the spin and isospin variables of the –th nucleon, and

(30)

Note that can be moved to the left of to obtain the second line of Eq. (29) as, on account of the property

(31)

one needs to antisymmetrize only.

Summing over the states belonging to the Fermi sea for each independently leads to an expression which does not depend on the number of particle . There is no violation of the Pauli principle because of the antisymmetrization of . More specifically, each term of the r.h.s. of Eq. (31) coming from the antisymmetrization of is Pauli violating, but the total sum it is not. On the other side, the independence of has very useful consequences. One of them is that the numerator of can be easily recognized as the product of the denominator times the sum of linked cluster terms. In addition, the FR cluster expansion is exact for any number of particles, not just in the thermodynamic limit like, for example, the Mayer expansion. This property has been exploited in FHNC calculations of finite nuclear systems like nuclei ariasdesaavedra_bisconti_co_fabrocini () or nucleon confined in periodical box PBFHNC ().

The operatorial –body Fermi distribution function includes a direct term corresponding to in Eq. (17) and a number of exchange terms generated according to the algebra of the exchange operators .

The basic statistical (exchange) correlation is described by the Fermi gas one–body density matrix

(32)

where the Slater function is given by

(33)

Diagrammatically, the exchange correlation , referred to as “exchange bond”, is represented by an oriented solid line. The –algebra implies that the exchange bonds form closed loops which never touch each other. If is made of scalar correlations only, then all nucleons in a given exchange loop must be in the same spin–isospin state and the Fermi distribution operators of Eq. (30) reduces to the standard Fermi gas distribution functions. For example, the two–body Fermi distribution function is given by

(34)

As an example, consider the two-body cluster contribution. From Eq. (29) it can be written as

(35)

The sum over the states belonging to Fermi sea implies a sum over the spin-isospin quantum numbers, which amounts to computing the trace of the spin and isospin operators appearing in . The trace is normalized to unity, as summation over the momenta leads to the appearance of a factor in both the direct and exchange term. The final result is

(36)

iii.2.1 Diagrammatic rules

The diagrams consist of dots (vertices) connected by different kinds of correlation lines. Open dots represent the active (or interacting) particles ( and ), while black dots are associated with passive particles, i.e. those in the medium. Integration over the coordinates of the passive particles leads to the appearance of a factor .

Active correlations must be treated differently from the passive ones, as the components of the two–body potential may be singular, thus leading to divergent integrals. In the diagrammatic expansion of , the quantity is represented by a thick solid line, denoted as “interaction line” and depicted in Fig. 3.

Exchange lines form closed loops, oriented clockwise or counterclockwise, the simplest of which is the two–body loop yielding a contribution . In addition to the extra factors coming from the algebra arising from the spin–isospin structure of the corresponding cluster term, an –vertex loop carries a factor , where is associated with each exchange operator and is due to the presence of spin–isospin species of the loop, combined with the existence of different orientation and to the minus sign coming from the permutations. The two–body loop is an exception to this rule, because there is only one such loop. Therefore, the corresponding factor is rather than .

The correlation bonds of Figs. (23) cannot be superimposed to each other. They can only be superimposed to exchange bonds.

The allowed diagrams are all linked, as a result of the linked cluster property discussed above.

Figure 3: Graphical representation of an interaction line.

A typical diagram of the FR cluster expansion is sketched in Fig. 4. Its contribution to the potential energy expectation value is given by

(37)

where the trace is carried out in the spin-isospin spaces of particles 1, 2 and 3. The factor comes from the relation , due to translation invariance, and is a symmetry factor. The four orderings appearing on the r.h.s. of Eq. (37) correspond to the two possible positions of , on either the left- or right-hand side of the operator .

Figure 4: Example of diagram appearing in the cluster expansion of .

iii.2.2 FHNC/SOC approximation

All the linked cluster diagrams or sub–diagrams built with scalar passive bonds only, with the only exception of the so called elementary diagrams, can be summed up in closed form by solving the FHNC equations fantoni_rosati (); fabrocini_fantoni (). The contributions associated with the elementary diagrams can be formally included in the FHNC equations, but there is no exact procedure to sum all of them. They can only be taken into account approximatively, by explicit calculation of the n–point () basic structures (FHNC/n approximation). However, it is well known that in nuclear matter calculations the FHNC approximation provides very accurate results.

On the other hand, diagrams having one or more passive operatorial bonds are calculated at leading order only. Such an approximation is justified by the observation that operatorial correlations are much weaker than the scalar ones. Based on this feature, one would be tempted to conclude that the leading order amounts to dressing the interaction line with all possible FHNC two–body distribution functions. This is not true as, besides the short range behavior, the intermediate range behavior of NN correlations also plays an important role that needs to be taken into account. In particular, tensor correlations, and to some extent also exchange correlations, have a much longer range than the central ones.

In order to handle this problem, summing the class of chain diagrams turns out to be to be of great importance. To see this, let us consider an extreme example of a long range correlation, namely a function that heals to zero as , implying that its Fourier transform behaves as in the long wavelength limit. Chain diagrams of –bonds are calculated by means of the convolution integral of the various in the chain. In Fourier space convolution integrals are given by products of . One can easily verify that, in the long wavelength limit, any chain diagram is more singular than . On the contrary, the sum of all the chain diagrams has exactly the same degree of singularity. Hence, summing up the series of chain diagrams takes care of long range correlations 111This feature is critical to the calculation of the long wavelength limit of the static structure function and the phonon excitations..

The above issue is taken care of by summing up the Single Operator Chains (SOC) in the corresponding FHNC/SOC approximation wiringa_pandha_1 (); wiringa_pandha_2 (). SOC are chain diagrams in which any single passive bond of the chain has a single operator of the type or , with , or FHNC–dressed versions of them. Note that if a single bond of the chain is of the scalar type then the spin trace of the corresponding cluster term vanishes, as the Pauli matrices are traceless. Then the SOC is the leading order, and at the same time it includes the main features of the long range behavior of tensor and exchange correlations.

The calculation of SOC, as that of FHNC chains, is based upon the convolution integral of the functions corresponding to two consecutive bonds. Unlike FHNC chains, however, the SOC have operatorial bonds. Therefore, the basic algorithm is the convolution of two operatorial correlations having one common point. Let us consider two such correlation operators, say and . Their convolution gives rise to a correlation operator of the same algebraic form :

(38)

where the functions depend on the internal angles of the triangle . The above equation includes also the convolution of the scalar correlations, which is already taken into account by the FHNC chain equations. Hence, . If one of the bonds is scalar and the second is operatorial the convolution vanishes, i.e. . The explicit expressions of can be found in Refs. wiringa_pandha_1 (); fabrocini_fantoni ().

The ordering of the operators within an SOC is immaterial, because the commutator is linear in and , and Pauli matrices are traceless. The only orderings that matter are those of passive bonds connected to the interacting points or . The reason is that the interaction line may have up to four operators. Therefore, or may be reached by up to five operators and one has to take into account the different orderings. The underlying spin algebra is lengthy but straightforward, and it is given in Ref. wiringa_pandha_1 (). As an example, consider the cluster diagram of Fig. 4 and the corresponding cluster term of Eq. (37). The two orderings and give rise to the same trace, which in the case of turns out to be . The full expression of can be easily extracted from the energy term displayed in Eq. (7.7) of Ref. wiringa_pandha_1 (), and written in terms of the matrices and and the vector , defined as follows

(39)

where the sign applies if

(40)

while the sign applies if

(41)

The –matrices are associated with normal orderings, like , whereas the –matrices apply to alternate orderings, like .

A second important contribution which is included in FHNC/SOC approximation is the leading order of the vertex corrections. They sum up the contributions of sets of subdiagrams which are joined to the basic diagrammatic structure in a single point. Therefore, a vertex correction dresses the vertex of all the possible reducible subdiagrams joined to it. The FHNC equations for the full summations of these one–point diagrams are given in Ref. fabrocini_fantoni (). In the FHNC/SOC approximation they are taken into account only at the leading order, i.e. including single operator rings (SOR), which are nothing but loops of SOC. Vertex corrections play an important role for the fulfillment of the sum rules.

iii.2.3 Two–body and three–body distribution functions

The expectation value (28), can be conveniently rewritten in the form

(42)

where

(43)

are the operatorial components of the two–body distribution function.

The FHNC diagrams are divided in 4 separate classes, characterized by the type of bonds ending at the vertices associated with particles and . The different types of vertices are denoted “d” for direct, i.e. involving no exchange lines, “e” for exchange, i.e. the vertex of an exchange loop, and “c” for cyclic, i.e. the vertex of an exchange line. Using this notation we can write,

(44)

The two–body distribution functions satisfy the following sum rules

(45)

Note that the above relations also hold true for the distribution functions of the Fermi gas model, as well as for those obtained retaining the correlations only. Another sum rule is given by the expectation value of the kinetic energy, which can be written in three equivalent forms, known as Pandharipande–Bethe (PB), Jackson–Feenberg (JF) and Clark–Westhaus (CW). In an exact calculation they would all give the same results. Numerical differences between them gauge the degree of accuracy of the approximations employed in the calculation.

The generalization of Eq. (28) to the case of a three-body potential, e.g. the UIX model, reads

(46)

As for the case of the two–body distribution functions , it is useful to define three–body distribution functions , reflecting the operatorial structure of given in Eqs. (8) and (13). Let us write as a sum of spin-isospin three–body operators multiplied by scalar functions, depending on the relative distances only

(47)

From Eqs. (8) and (13) it follows that the sum of the above equation involves 19 operators. The expectation value of can be written as

(48)

with

(49)

In Ref.carlson_pandha_wiringa (), the above expectation value has been computed in FHNC/SOC. The cluster expansion and the corresponding diagrammatic rules of the cluster diagrams are very similar to those outlined in the case of the two-body potential, with the only difference of three external points, instead of two, all linked by interaction lines.

Figure 5: Diagrammatic representation of Eq. (51): the two-body density-dependent potential includes the effects of both the bare three-body potential and the correlation and exchange lines. While dresses the line joining particles and , the dressing being depicted by a line with a big bubble in the middle, dresses the lines , , and .

Iv Derivation of the effective potential

Our work is aimed at obtaining a two-body density-dependent potential that mimics the three-body potential. Hence, our starting point is the requirement that the expectation values of and of be the same:

(50)

implying in turn (compare to Eqs.(42) and (48))

(51)

A diagrammatic representation of the above equation, which should be regarded as the definition of the , is shown in Fig. 5. The graph on the left-hand side represents the three-body potential times the three-body correlation function, integrated over the coordinates of particle . Correlation and exchange lines are schematically depicted with a line having a bubble in the middle, while the thick solid lines represent the three-body potential. The diagram in the right-hand side represents the density-dependent two-body potential, dressed with the two-body distribution function. Obviously, has to include not only the three-body potential, but also the effects of correlation and exchange lines.

The left-hand side of Eq.(51) has been evaluated in carlson_pandha_wiringa () within the FHNC/SOC scheme. Here we discuss the derivation of the explicit expression of the two-body density-dependent potential appearing in the right-hand side of the equation. The procedure consists of three different step, each corresponding to a different dressing of the diagrams involved in the calculation

For each of these steps the final result is a density-dependent two-body potential of the form

(52)

where, depending on the step, the can be expressed in terms of the functions appearing in the definition of the UIX potential, the correlation functions and of the Slater functions.

Step I. Bare approximation

As a first step in the derivation of the density-dependent potential one integrates the three-body potential over the coordinate of the third particle

(53)

Diagrammatically the above equation implies that neither interaction nor exchange lines linking particle with particles and are included. Only the two-body distribution function is taken into account in the calculation of the expectation value of

(54)

Note that only the scalar repulsive term and one permutation of the anticommutator term of the three-body potential provide non vanishing contributions, once the trace in the spin–isospin space of the third particle is performed.

As shown in Fig 10, the contribution of the density-dependent potential to the energy per particle of SNM and PNM is more repulsive than the one obtained from the genuine three-body potential UIX. Thus, the scalar repulsive term is dominant when the three-body potential is integrated over particle .

iv.1 Step II. Inclusion of statistical correlations

As a second step we have considered the exchange lines that are present both in and . Their treatment is somewhat complex, and needs to be analyzed in detail.

Consider, for example, the diagram associated with the exchange loop involving particles , and , depicted in Fig. 6. Its inclusion in the calculation of the density-dependent two-body potential would lead to double counting of exchange lines connecting particles 1 and 2, due to the presence of the exchange operator in . This problem can be circumvented by noting that the antisymmetrization operator acting on particles , and can be written in the form

(55)

in which the exchange operators contributing to the density-dependent potential only appear in the first term of the right-hand side.

On the other hand, the second term in the right-hand side of Eq. (55) only involves the exchange operators , whose contribution is included in and must not be taken into account in the calculation of .

Two features of the above procedure need to be clarified. First, it has to be pointed out that it is exact only within the SOC approximation that allows one to avoid the calculation of commutators between the exchange operators and and the correlation operators acting on particles and . The second issue is related to the treatment of the radial part of the exchange operators. Although it is certainly true that one can isolate the trace over the spin-isospin degrees of freedom of particle 3, arising from and , extracting the Slater functions from these operators is only possible in the absence of functions depending on the position of particle 3 pandha_scambi ().

Figure 6: Three particle exchange loop.

However, this restriction does not apply to the case under consideration, as both the potential and the correlations depend on and . As a consequence, retaining only the and exchange operators involves an approximation in the treatment of the the Slater functions, whose validity has been tested by carrying out a numerical calculation.

By singling out the radial dependence of the exchange operators, and by computing the inverse of the operator , where denotes the spin-isospin part of , it is possible to find a “Slater Exact” density-dependent potential whose calculation does not involve any approximations concerning the Slater functions. It can be rewritten in the form

(56)
Figure 7: Contributions of the density-dependent potential to the energy per particle (see Eqs. (56) and (57)), arising from the scalar term of UIX (upper panel) and from the anticommutator term (lower panel).

Note that in the above equation we have omitted all correlations functions, whose presence is irrelevant to the purpose of our discussion. The density-dependent potential obtained from Eq.(56) must be compared to the one resulting from the approximation discussed above, which (again neglecting correlations) leads to the expression

(57)

where “S. A.” stands for Slater Approximation. We have computed and for SNM within the FHNC/SOC scheme, for both the scalar and the anticommutator terms of the UIX potential.

The results, plotted in Fig. 7, clearly show that Eq.(57) provides an excellent approximation to the exact result for the exchanges of Eq. (56). Hence it has been possible to use Eq. (57) also to compute the contribution coming from the commutator of the UIX potential, avoiding the difficulties that would have arisen from an exact calculation of the exchanges.

The second step in the construction of the density-dependent potential is then

(58)

which is a generalization of the bare potential of Eq. (53).

Figure 10 shows that taking exchanges into account slightly improves the approximation of the density-dependent potential. However the differences remain large because correlations have not been taken into account.

iv.2 Step III. Inclusion of dynamical correlations

The third step in the construction of the density-dependent potential amounts to bringing correlations into the game. We have found that the most relevant diagrams are those of Fig. 8.

Figure 8: Diagrams contributing to the density-dependent potential. The dashed lines with diamonds represent the first order approximation to , discussed in the text. Only diagrams with at most one operator attached to a given point are taken into account.

Note that, in order to simplify the pictures, all interaction lines are omitted. However, it is understood that the three-body potential is acting on particles , and . Correlation and exchange lines involving these particles are depicted as if they were passive interaction lines. Moreover, in order to include higher order cluster terms, we have replaced the scalar correlation line with the Next to Leading Order (NLO) approximation to the bosonic two-body correlation function:

(59)

The full bosonic or might be used instead of the NLO approximation. However, including higher order terms would have broken our cluster expansion. The correction to of Eq. (59), whose diagrammatic representation is displayed in Fig. 9, can indeed be considered to be of the same order as the operatorial correlations.

Figure 8 shows that the vertices corresponding to particles 1 and 2 are not connected by either correlation or exchange lines. All connections allowed by the diagrammatic rules are taken into account multiplying the density-dependent potential by the two-body distribution function, according to the definition of Eq.(51).

We have already discussed the exchange lines issue, coming to the conclusion that only the exchanges and have to be taken into account. This is represented by the second diagram, where the factor is due to the symmetry of the three-body potential, that takes into account both and .

Figure 9: NLO approximation to the bosonic two-body correlation function.

The explicit expression of obtained including the diagrams of Fig. 8 can be cast in the form

(60)

where denotes the sum of non central correlations

(61)

Note that, in principle, an additional term involving the anticommutator between the potential and the correlation function should appear in the second line of the above equation. However, due to the structure of the potential it turns out that

(62)

The calculation of the right-hand side of of Eq. (60) requires the evaluation of the traces of commutators and anticommutators of spin-isospin operators, as well as the use of suitable angular functions needed to carry out the integration over .

Figure 10: Contributions of the density-dependent potential to the energy per particle of SNM (upper panel) and PNM (lower panel), compared to the expectation value of the genuine three-body potential UIX: .

As for the previous steps, we have computed the contribution of the density-dependent potential to the energy per particle. The results of Fig. 10 demonstrate that the density-dependent potential including correlations is able to reproduce the results obtained using genuine three-body UIX to remarkable accuracy.

To simplify the notation, at this point it is convenient to identify

(63)

Note that the above potential exhibits important differences when acting in PNM and in SNM. For example, in SNM for , while in PNM for