Saturation of nuclear matter and roles of many-body forces: nuclear matter in neutron stars probed by nucleus—nucleus scattering

Right © The Author(s) 2016. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/4.0/ ), which permits unrestricted reuse,distribution, and reproduction in any medium, provided the original work is properly cited.Funded by SCOAP3 DOI 10.1093/ptep/ptw072


Introduction
Yoichiro Nambu left great footprints in physics. His greatest achievement will be his genius idea of spontaneous symmetry breaking (SSB) inspired by the BCS theory for superconductivity in condensed matter. The famous Nambu-Jona-Lasinio (NJL) papers [1,2] marked a great milestone in particle and nuclear physics. The idea proposed in the NJL papers has led to the idea of dynamical breaking of chiral symmetry in quantum chromodynamics (QCD) of the standard model, and the original NJL model is now understood as a low energy effective theory of QCD of a highly nonperturbative nature.
Since Nambu's contributions to various fields in physics are reviewed by other authors in the review articles for the present memorial symposium, this paper will be devoted to some specific aspects of his contribution, particularly to the repulsive nature of the nuclear force and the saturation property of nuclear matter, and related topics relevant to the saturation property of nuclear matter probed by nucleus-nucleus scattering will be presented.
To be more specific, we introduce a recently proposed unique scenario of extracting important information about the saturation property of nuclear matter and the stiffness of neutron stars through the analysis of nucleus-nucleus elastic scattering in laboratories. Key ingredients of the scenario are the short range repulsive three-body forces that act among nucleons and baryons in nuclear and neutron star matter. The idea is that the unknown strength of the three-body force can be fixed quite precisely through analyses of elastic scattering experiments in laboratories by the latest microscopic theory for interactions between finite nuclear systems.
The neutron stars J1614-2230 [3] and J0348-0432 [4] have had a great impact on the maximum mass problem; their observed masses are (1.97 ± 0.04)M and (2.01 ± 0.04)M , respectively, much larger than was observed before (say 1.4 ∼ 1.5)M . These large masses have added a severe condition for the stiffness of the equation of state (EOS) of neutron star matter that serves as a great challenge to nuclear physics.
In order to clarify the implication of the scenario, we briefly review a microscopic theory of interaction between finite nuclei, called the double-folding model (DFM), and Brueckner's G-matrix theory; these are the keys to connect infinite nuclear matter, including the baryonic matter in neutron stars, and nucleus-nucleus scattering experiments in laboratories. To this end, we first present a brief history of the theoretical study of effective nucleon-nucleon (NN) interaction in the nuclear medium derived from Brueckner's G-matrix theory and its application in constructing nucleusnucleus interactions in terms of various types of DFM.

Basic properties of finite nuclei
Before discussing the latest status of the effective interactions and their application to nucleusnucleus scattering and the saturation of nuclear matter, we first describe the basic properties of atomic nuclei.

Nuclear force and repulsive core
The origin of the "nuclear force" (or NN interaction) that sustains protons and neutrons in forming an atomic nucleus with an extremely small radius of less than 10 −14 m = 10 fm was first proposed by Hideki Yukawa [5] as the exchange of a pion (π meson) between nucleons (the one-pion exchange OPE). The OPE is known to generate the outermost (long range) part of the nuclear force because of the light mass of the exchanged pion, m π 140 MeV, with the so-called Yukawa tail ∼ e −m π r /r of the nuclear force. The nuclear force has a stronger attraction in the middle range (0.5 ∼ 1.0 fm) region that is now understood by the exchange of heavier mesons such as σ , ω, and ρ mesons, as well as by the exchange of multiple pions.
In addition to the attraction in the middle and long range parts of the nuclear force, a series of experimental measurements of nucleon-nucleon scattering at high energies indicates the existence of a strong repulsion at short distances of less than about 0.5 fm. The simplest version of the repulsive part of the nuclear force is the so-called hard core of infinite strength, which was first introduced phenomenologically by Jastrow in 1950 [6] to explain the nucleon-nucleon scattering experiments at high energies. To date, various types of repulsive core, either of the hard-core or soft-core type, have been implemented in realistic nuclear force models widely used in nuclear structure and reaction studies, although the origin of the short range repulsive core remains to be determined.
One pioneering work on a possible origin of the short range core of the nuclear force was by Nambu, who proposed a vector meson called the omega (ω) meson in 1957 [7] as a possible source of the repulsive core, in addition to the source of the attraction in the middle range of the nuclear force. The pioneering work by Nambu on the origin of the repulsive core was followed by various types of quark model [8,9], but it is only recently that a milestone in understanding the origin of the nuclear force, including the repulsive core, was marked in the epoch-making work by Ishii and his collaborators [10] with first-principle calculations based on QCD, the standard model for strong interaction. They derived the radial shape and strength of the nuclear force in lattice QCD calculations 2/25 PTEP 2016, 06A106 Y. Sakuragi that reproduce well the basic properties of the nuclear force, not only at the qualitative level but also quantitatively, derived from the meson exchange mechanism supplemented by a phenomenological repulsive core at short distances.

Nuclear force and saturation of nuclear matter
The nuclear radius observed in various experiments is found to be proportional to the cube root of the mass number A(= Z + N ), R r 0 × A 1/3 , which implies that nuclear volume V = 4πR 3 /3 is proportional to A and, hence, the average nucleon density in the nucleus ρ = A/V = 3/(4πr 3 0 ) is approximately constant, irrespective of the mass number A of the nucleus. The constant value ρ 0 0.17 fm −3 is called the saturation density. The mass of an atomic nucleus, M (Z, N ), with atomic number Z and neutron number N , defines the binding energy of the nucleus by B(Z, N ) = Zm p + Nm n − M (Z, N ) c 2 , m p and m n being the mass of a proton and a neutron respectively. The experimental measurement of nuclear masses shows that the binding energy per nucleonĒ = B(Z, N )/A has an almost constant value ofĒ 8 ± 0.5 MeV irrespective of the number of nucleons in nuclei; this is called the saturation of binding energy in atomic nuclei. This indicates that the total binding energy of a nucleus is approximately proportional to its mass number A (and hence to its The saturation of the nucleon density and that of the binding energy per nucleon in atomic nuclei are closely related to the nature of the nuclear force, particularly of its short range nature and the existence of the short range repulsive core. From the following intuitive consideration, it is easily understood that saturation of the binding energy per nucleon is a direct consequence of the saturation of nucleon density in nuclei and the short ranged nature of the nuclear force. First of all, one should note that the maximum range of the nuclear force of about 1.5 ∼ 2.0 fm is approximately the same order as the mean distance between nucleons in atomic nuclei evaluated from the saturation density of ρ 0 0.17 fm −3 . Because of the constant nucleon density in a nucleus (except in the region near the nuclear surface), the number of nucleons surrounding a given nucleon in uniform nuclear matter is also constant. This indicates that the binding energy for a nucleon by its surrounding nucleons in uniform nuclear matter must be almost constant because the nuclear force is of short range nature and does not act beyond the neighboring nucleons. In other words, the total binding energy of nucleons in a nucleus should be proportional to the number of nucleons (mass number A), which is also proportional to the nuclear volume because V ∝ R 3 ∝ A. This term of the binding energy is then called the volume term.
One might expect that the proportional constant a V of the volume term would be identified with the measured constant value of the binding energy per nucleon, a V ≈Ē 8 ± 0.5 MeV. However, this is not true at all and the measured constant value 8 ± 0.5 MeV is known to be a consequence of large corrections due to the finiteness of nuclei (or the effect of surface tension), the Coulomb energy, and the symmetry energy originating from the difference in neutron and proton numbers in a nucleus, all of which are the correction for reducing the binding energy. Subtracting those correction terms, the genuine volume term of the binding energy per nucleon is found to be a V 16 MeV in symmetric nuclear matter of infinite extension in space, namely having no surface. The binding energy per nucleonĒ 16 MeV is called the saturation energy of symmetric nuclear matter.
The saturation of the nucleon density and that of the binding energy per nucleon in atomic nuclei suggest the existence of a fundamental mechanism sustaining atomic nuclei in their size and binding and preventing the nuclei from shrinkage or collapse due to the strong attraction of the nuclear force. It is the short range repulsion of the nuclear force that is one of the candidates for the origin of the 3/25 PTEP 2016, 06A106 Y. Sakuragi saturation mentioned above. However, the short range repulsion itself of the nuclear force between two interacting nucleons, namely the two-body repulsion, is not sufficient to explain the saturation property observed in nuclear experiments as well as explaining the maximum mass of neutron stars M ∼ = 2M obtained from the latest astronomical observations [3,4]. A major subject to be discussed in the present paper will be that the effect of many-body forces of repulsive nature plays a decisive role in realizing the saturation of nuclear matter as well as of neutron star matter.

Effective interaction in the nuclear medium
Understanding nuclear structure and reactions microscopically starting from the basic nucleonnucleon (NN) interaction has been the fundamental subject of nuclear physics. The basic properties of NN interaction in free space, say its shape, strength, and spin-isospin dependence, are well understood phenomenologically from the observed properties of the deuteron, which is the only bound state of the two-nucleon system, and from NN scattering experiments at a wide range of scattering energies with various spin orientations of colliding nucleons. A number of theoretical models for NN interaction in free space have also been proposed, ranging from purely phenomenological models to those based on various meson theories and quark models. A recent derivation of the nuclear force in terms of first-principle lattice QCD calculation [10] based on the standard model for the strong interaction QCD marks a milestone in understanding the nuclear force, as mentioned in Sect. 1.
However, it should be emphasized that the interaction between nucleons in nuclei is not equivalent to that in free space and is strongly modified by the surrounding nucleons. This is called the medium effect or, in other words, the many-body effect. The NN interaction modified in the nuclear medium is called the effective NN interaction.

Brueckner's G-matrix and saturation properties of nuclear matter
Brueckner's G-matrix theory (see [11,12] and references therein) is one of the most successful theories for handling the short range correlation (SRC) of the bare NN interaction and generating an effective NN interaction in the nuclear medium. In the G-matrix theory, a pair of nucleons is embedded in infinite nuclear matter with a uniform nucleon density and one solves the "two-body" Schrödinger equation, called the Bethe-Goldstone equation, with the Pauli blocking operatorQ that prevents the two interacting nucleons from sitting in states already occupied by other nucleons of the surrounding nuclear matter. In other words, the interacting nucleon pair are virtually excited to unoccupied highly excited states in the intermediate states in the series of ladder diagrams under the effect of Pauli blocking by the occupied nucleons in nuclear matter. This implies that the motion of the nucleon pair is strongly restricted, though indirectly, by the surrounding nucleons through the Pauli blocking operator Q, particularly in high density nuclear matter. In this sense, the medium effect is interpreted as the many-body effect. Let us consider a pair of nucleons in momentum space, one being a bound nucleon with momentum p (in units of ) in nuclear matter with Fermi momentum k F and the other being an incoming nucleon with momentum k, where |p| < k F and |k| > k F in the initial state before the interaction. It is noted that the Fermi momentum k F is proportional to the cube root of the uniform density ρ of the nuclear matter, k F ∝ ρ 1/3 . Relative and center-of-mass momenta are given as k r = (k−p)/2 and K c = k+p, respectively. With an NN interaction V , the G-matrix equation giving the scattering of the pair in 4 where e(q) is a single-particle (s.p.) energy in an intermediate state with momentum q, and Q(q 1 , q 2 ) is the above-mentioned Pauli blocking operator defined by The starting energy ω is given as a sum of the energy E(k) of the propagating nucleon and the s.p. energy e(p) = 2 2m p 2 + U (p, e(p)) of a bound nucleon [13]. The G-matrix calculations are often performed with the so-called continuous choice for intermediate nucleon spectra, playing an essential role especially for imaginary parts of G-matrices. This choice means that s.p. energies are calculated self-consistently not only for bound states (q ≤ k F ) but also for intermediate states (q > k F ). The scattering boundary condition iε in the denominator leads to complex G-matrices, summation of which gives the complex s.p. potential U (q, e(q)). The plausible way is to use this self-consistent complex potential in the G-matrix equation (2). Avoiding numerical complexities in such a procedure, however, its real part U R (q) = Re U (q, e(q)) is usually used in the self-consistency process [14][15][16][17][18].
The G-matrix also depends on the NN interaction model in free space V adopted in the calculation, because of the difference in the properties of central, spin-orbit, and tensor forces and their spinisospin dependence, as well as in the off-energy-shell properties, although all of them reproduce more or less the basic properties of two-nucleon systems, both in the bound state (deuteron) and scattering states observed in experiments. (One should note that experimental data of the two-body NN system give us only information about the on-energy shell because of the energy conservation of the two-body system.)

Saturation curves and stiffness of symmetric nuclear matter
The single-particle energy in nuclear matter is calculated in the Brueckner-Hartree-Fock (BHF) framework with the use of the G-matrix interaction. This is nothing but the binding energy of a nucleon in nuclear matter, which corresponds to the volume term of the binding energy per nucleon E = E/A in a finite nucleus. The single-particle energy depends on the density ρ of the nuclear matter, and the curve ofĒ(ρ) as a function of ρ is called the saturation curve of the symmetric nuclear matter. (An illustrative example is shown in Fig. 1.) The saturation curve should have a minimum valueĒ 0 at a certain value of the density ρ 0 ; the point (Ē 0 , ρ 0 ) in the two-dimensional The history of the study of effective interactions and the saturation curve for symmetric nuclear matter tells us that it is impossible to reproduce the saturation point by any effective interactions evaluated from the non-relativistic G-matrix calculations starting from various kinds of NN interaction models as long as only the two-body NN interaction is taken into account in the calculation. Each model gives a saturation point apart from the empirical one, and it is known that the evaluated saturation points lie along an almost straight line called the Coester line [19] that passes by the empirical saturation point at a non-negligible distance. Thus, it is one of the fundamental subjects in nuclear physics to reproduce the saturation point in symmetric nuclear matter. A key is the role of the many-body forces in dense nuclear matter. As will be shown later, the introduction of three-body and four-body forces that act among three or four nucleons simultaneously in high density nuclear matter can resolve the problem. Within the non-relativistic framework, it is inevitable to introduce such many-body forces to reproduce the correct saturation point.
In addition to the saturation point, the shape of the saturation curve is also important for understanding the dynamical properties of nuclear matter with the variation of its density. It is also an important subject to investigate the static and dynamic response of nuclear matter at high densities well above the saturation density (ρ ρ 0 ) under the high pressure of gravity in high density neutron stars. The second derivative of the saturation curve gives a measure of the stiffness of the nuclear matter. The incompressibility K SNM of nuclear matter is defined by the second derivative ofĒ with respect to ρ at the saturation density ρ = ρ 0 , as such as a three-body force of a repulsive nature, although the evaluated values of K SNM ≈ 200 MeV ∼ 300 MeV are far from conclusive for determining the precise stiffness of symmetric nuclear matter and discussing the detailed properties of such an additional repulsive effect. Besides the K SNM value, what is more important is the shape of the saturation curve itself, particularly at higher densities beyond the saturation density ρ ρ 0 , since K SNM gives only the curvature of the saturation curve at or near the saturation density ρ ∼ = ρ 0 by definition. It only provides us with information about the stiffness of symmetric nuclear matter near the saturation density. Thus, it is quite important to obtain information about the stiffness of nuclear matter at densities much larger than the saturation value, say, for example, ρ > 2ρ 0 ∼ 5ρ 0 . This is just the region where one expects high density nuclear matter to be realized inside a neutron star, though it is neutron matter rather than symmetric nuclear matter. It is one of the important goals of nuclear physics to predict reliable saturation curves both of symmetric nuclear matter and neutron matter over a wide range of nucleon density, including the high density region of neutron star matter.
In the following, we will see a recent attempt to determine the strength of the three-body and four-body force effects of a repulsive nature through precise analyses of nucleus-nucleus elastic scattering in laboratories. The nucleus-nucleus interactions are derived microscopically from modern G-matrix effective NN interactions with three-body and four-body force effects included. This makes it possible to directly discuss the stiffness of high density nuclear matter, as well as that of neutron stars in the context of the effects of three-body and four-body forces that can be probed and tested by nucleus-nucleus elastic scattering in laboratories on Earth.

Microscopic models for interaction between finite nuclei
It is one of the important subjects in nuclear physics to understand microscopically the interactions between complex nuclei consisting of nucleons, starting from the bare NN interaction in free space. In constructing the nucleon-nucleus and nucleus-nucleus interactions starting from the bare NN interaction having a strong short range repulsive core, one has to correctly handle the short range correlation (SRC) originating from the short range repulsive core of the NN interaction. This is equivalent, in principle, to exactly solving the many-body dynamics in the full model space up to extremely high excitation energies, which is almost impossible to perform in practice except for few-body systems consisting of nucleons of at most A 8.
One alternative way of handling the SRC between interacting nucleons in nuclear many-body systems is to introduce an effective NN interaction that contains the SRC effectively in the multiple scattering framework, such as the T -matrix, R-matrix, and G-matrix interactions. The short range singularities in the free-space NN interaction are smoothed out in those effective interactions.

G-matrix folding model with local density approximation
Brueckner's G-matrix is one of the simplest but most reliable effective NN interaction models that includes all the ladder diagrams of the NN interaction with the medium effect from surrounding nucleons in nuclear matter being taken into account. This implies that the SRC induces virtual excitations of the interacting NN pair into high energy states outside the Fermi sphere and its effect is included, within the nuclear matter approximation, in the resultant G-matrix effective NN interaction. The G-matrix is then interpreted as an effective NN interaction in the "cold" or "frozen" nuclear matter of a uniform nucleon density. This is the basis of the idea of the folding models discussed below for constructing nucleon-nucleus or nucleus-nucleus interaction potentials in terms of the 7/25 PTEP 2016, 06A106 Y. Sakuragi G-matrix. In the folding model, the G-matrix is folded over the nucleon density distributions of the interacting finite nuclear systems in their ground states.
Since the G-matrix v NN is evaluated in nuclear matter of a uniform constant density, the evaluated G-matrix depends on the nucleon density ρ of the nuclear matter. Because of the finite size of the nucleus, the environment, or the local density, for an interacting nucleon pair embedded in finite nuclear systems changes rapidly with the spatial positions where the nucleon pair exist. Thus, the folding model adopts the local density approximation (LDA), in which the uniform density of the nuclear matter for evaluating the G-matrix is assumed to be the same as the local density at the spatial positions of the interacting nucleons in finite nuclei. In other words, we use a different strength of the G-matrix depending on the local density at the positions of interacting nucleons in finite nuclei in calculating the nuclear potential in the folding model. The G-matrix also depends on the starting energy ω, which corresponds to the incident energy per nucleon, E/A, in the case of nucleus-nucleus scattering.

Optical potential and folding model with complex G-matrix
The G-matrix is real for bound states because both interacting nucleons are bound inside the Fermi sphere in the initial and final states of the interaction, although they can be excited virtually to continuum states outside the Fermi sphere in the intermediate stages of the interaction. This kind of real G-matrix is widely used in various nuclear structure calculations. On the other hand, the G-matrix derived under the scattering boundary condition becomes complex because both of the interacting nucleons can be excited to continuum (scattering) states outside the Fermi sphere in the final states of the interaction. In nucleus-nucleus scattering, for example, the interacting nucleon pair, one in the projectile nucleus and the other in the target one, are in their relative scattering states. Thus, the G-matrix becomes complex and the complex G-matrix should be used in constructing the nucleusnucleus potential microscopically by the folding model. In fact, the nucleon-nucleus or nucleusnucleus potential evaluated phenomenologically from the elastic scattering experiment is known to be a complex potential called the optical potential. The optical potential is a complex potential having real and imaginary parts, each of which is assumed to have a simple geometrical functional form, such as the three-parameter Woods-Saxson form, , and the parameter values of the real and imaginary parts are to be optimized so that the angular distribution of the elastic scattering cross section calculated with the optical potential reproduces that of the experimentally observed elastic scattering cross section.
The imaginary part of the optical potential represents the loss of flux from the elastic scattering channel. In a nucleon-nucleus or nucleus-nucleus collision, many kinds of nuclear reactions take place and some fraction of the incident flux, which initially comes into the entrance channel (i.e. to the elastic scattering channel), escapes from the entrance elastic scattering channel to non-elastic reaction channels, such as inelastic scattering and various kinds of rearrangement, break-up and fusion reactions, etc. The flux loss can be represented approximately by the imaginary part of the interaction potential.
The complex folding model potential obtained from the complex G-matrix is a promising candidate for the microscopic interaction model to explain the complex optical potentials derived phenomenologically from nucleon-nucleus or nucleus-nucleus elastic scattering experiments, as will be shown below in some detail.
In the following sections, we will concentrate our discussion on the double folding model (DFM) for constructing interaction potentials for nucleus-nucleus systems. This is because nucleus-nucleus scattering is an ideal tool for probing the saturation property of nuclear matter at high densities beyond the saturation density, say up to twice the normal density ρ ∼ 2ρ 0 . Combining the theoretical analyses of nucleus-nucleus scattering by the DFM and BHF study of the saturation curves of nuclear matter and neutron star matter using a common G-matrix, it will be possible to determine the unknown strength of the three-body force (or many-body forces, in general) very accurately and to directly discuss the so-called "maximum mass problem on neutron stars" together with the stiffness of nuclear matter and neutron stars.

Double folding model with G-matrix interactions
In the case of a nucleus-nucleus system, the folding model potential can be written as a Hartree-Fock type potential, where v D and v EX are the direct and exchange parts of the G-matrix v NN , while ij|v D |ij and ij|v EX |ij represent the matrix elements of the direct and knock-on exchange parts, respectively, of the G-matrix interaction. Here, |ij denotes the direct product of the ground state wave functions P and T of the interacting nuclei; i.e., |ij = | P ⊗ T . The microscopic nucleus-nucleus potential of the Hartree-Fock type potential with the G-matrix effective NN interaction defined in Eq. (5) is equivalent to the DFM potential with a density dependent effective interaction.
The direct part of the Hartree-Fock type potential is rewritten in the double-folding type integration where the direct part of the G-matrix v D is doubly folded over the nucleon densities of the interacting nuclei, P (projectile) and T (target), as Here, R is the relative coordinate between the centers of mass of the interacting nuclei P and T, and s ≡ R + r 1 − r 2 is the relative coordinate between a nucleon at the spatial position r 1 with respect to the center of mass (c.m.) of nucleus P and another nucleon at the spatial position r 2 with respect to the c.m. of nucleus T. They interact with each other through the direct part of the G-matrix, v D (see Fig. 2). The nucleon densities of the nuclei P and T are defined as where A P and A T are the mass numbers of the nuclei P and T, respectively. In Eq. (6), the direct part of the G-matrix v D depends on the local density ρ. In the present DFM calculation, the local density is defined by the sum of the densities of P and T, as ρ = ρ P + ρ T evaluated at the individual positions or the mid-point of the interacting nucleons, as will be discussed in more detail in the following subsection. This prescription for evaluating the local density is called the frozen density Definitions of the coordinates appearing in Eq. (6) for calculating the direct part of the DFM potential between the projectile nucleus (P) and the target nucleus (T).
approximation (FDA) and is widely used in most DFM calculations [20][21][22][23]. The FDA is also used in evaluating the local density in calculating the exchange part below. The exchange part of the DFM potential is a non-local potential in principle but is known to be approximated with good accuracy by the following local potential form of the so-called "Brieva-Rook type" [14][15][16]: Here, k(R) is the local momentum for the nucleus-nucleus relative motion defined by where M = A P A T /(A P + A T ), E c.m. is the center-of-mass energy, E/A is the incident energy per nucleon, m is the nucleon mass, and V coul is the Coulomb potential. The exchange part is calculated self-consistently on the basis of the local energy approximation through Eq. (10). Here, the Coulomb potential V coul is also obtained by folding the NN Coulomb potential with the proton density distributions of the projectile and target nuclei. The density matrix ρ(r, r ) is approximated in the same manner as in [24]: where k eff F is the effective Fermi momentum [25] defined by where C s = 1/4 following Ref. [26]. The detailed methods for calculating U D (the direct part) and U EX (the exchange part) are the same as those given in Refs. [27] and [28], respectively.

Frozen density approximation
One has to define the local density ρ at which the strength of the G-matrix is evaluated in the DFM calculations. In a nucleus-nucleus system, there are several prescriptions for defining the local 10/25 PTEP 2016, 06A106 Y. Sakuragi density ρ that defines the nuclear environment for an interacting nucleon pair embedded in the colliding nucleus-nucleus system. One is the average prescription, in which the local density is defined by either the geometric average or the arithmetic average of the local densities corresponding to the individual positions of the interacting nucleon pair: or, alternatively, the average of the local densities corresponding to the mid-point of the positions of the interacting nucleon pair. Another prescription, widely adopted in most DFM calculations, is the sum of the local densities corresponding to the individual positions of the interacting nucleon pair, or, alternatively, to the mid-point of the nucleons, both of which are called the frozen density approximation (FDA). In the FDA, the interacting nucleon pair is assumed to feel the local density defined as the sum of the densities of the colliding nuclei. The difference between the two prescriptions, the "average" and the "sum" (FDA), becomes significant when two nuclei come close to each other with a substantial overlap of the two density distributions. In the case of the FDA, the local density can reach the sum of the central densities of the interacting nuclei, which is close to twice the saturation density, say ρ ≈ ρ P (0) + ρ T (0) ≈ 2ρ 0 . On the other hand, the local density does not exceed the saturation density ρ ≤ ρ 0 in the average prescription.
One might postulate the average prescription to be more reasonable than the FDA because of the saturation of the nucleon density in finite nuclei, and might require that the local density should not exceed the saturation density ρ ≤ ρ 0 even when two nuclei come close to each other at very short distances. However, the strict condition that the local density ρ should not exceed the saturation density ρ 0 of "stable" nuclei implies that two stable nuclei in their ground states cannot heavily overlap each other in coordinate space. Otherwise either or both of the interacting nuclei are forced to be excited to highly excited states to fulfill the Pauli principle. In other words, the strict condition, ρ ≤ ρ 0 , inevitably leads to a strong rearrangement of nucleon configurations among the two interacting nuclei when they come close at short distances. This implies that most of the incident flux is removed from the elastic scattering channel to non-elastic reaction channels. In practice, these kinds of dynamic reaction processes may occur with high probabilities in high energy HI collisions, and almost all the incident flux will be lost from the elastic scattering channel.
Such kinds of hard collision, where the incident flux is completely removed from the elastic channel to the reaction channels, corresponds to the so-called shadow scattering, which is quantummechanical elastic scattering by a black-body target. In the case of shadow scattering, all the information about the internal region of the interaction potential between colliding nuclei is completely lost and no refractive scattering due to the attractive real potential should be observed in elastic scattering. In fact, in many cases, the elastic scattering between complex nuclei is dominated by strong absorption due to a huge number of non-elastic reactions occurring, when two nuclei come close to each other, and the elastic scattering cross section often shows typical characteristics of 11/25 PTEP 2016, 06A106 Y. Sakuragi strong absorption. In such cases, there is no hope of getting any information about the interaction potentials at short distances, and the prescriptions for defining the local density, either average or sum, are not relevant to the elastic scattering observables.
However, it is often the case that experimental data for elastic scattering of light HIs, such as the 16 O + 16 O and 12 C + 12 C systems at intermediate energies, show characteristic angular distributions typical of the refractive scattering generated by the attractive real potential at short distances that strongly participates in the elastic scattering cross section at large scattering angles [29]. This indicates that HI elastic scattering can probe important information about the strength of the interaction potential at very short distances [29,30] where two nuclei, both kept in their ground states, come close to each other at very short distances. This implies that there exists a non-negligible probability of two interacting nuclei being kept in the ground state at short relative distances [31], where the local density should be a simple sum of the undisturbed densities of the individual nuclei in their ground states that can reach up to twice the saturation density ρ ∼ 2ρ 0 .
Of course, it is natural to expect that the major part of the incident flux is removed from the elastic scattering channel to various non-elastic reaction channels. However, the refractive scattering widely observed in experiments clearly tells us that the incident flux is not completely removed from the elastic scattering channel and there remains, though only a tiny fraction, the elastic scattering component of the wave function for relative motion between the interacting nuclei at very short distances. The elastic scattering experiments can probe such a tiny but non-negligible component of the high local density region that carries quite valuable information about high density nuclear matter.
Nowadays, the superiority of the FDA [defined by Eqs. (14) and (15)] over the average prescription [defined by Eqs. (13)] has been made clear through experimental and theoretical analyses of a number of such refractive scatterings of light HI systems [21,26,28,32,33]. The average prescription largely overestimates the strength of the real potential at short distances and one needs to introduce an artificial renormalization factor to strongly reduce the real potential strength so evaluated in order to reproduce the elastic scattering experiments. On the other hand, no renormalization of the real potential evaluated with the FDA is necessary. Thus, it is widely recognized that the FDA is the best prescription for defining the local density to be used in the DFM calculation of HI potentials, as long as one uses reliable, realistic density dependent effective interactions [22,23,26,34,35] such as the G-matrices [23]. Therefore, we adopt the FDA in the calculation of DFM potential with the G-matrix throughout the following discussions.

A brief history of DFM studies of HI scattering
The first successful DFM for constructing nucleus-nucleus potentials was the M3Y folding model proposed by Satchler and Love in the 1970s, and its achievements in those days were summarized in a review article [33]. The M3Y (Michigan 3-range Yukawa) interaction [36] is a kind of Gmatrix represented by a linear combination of the three-range Yukawa functions. It has no density dependence and has only the real part. Thus, the M3Y folding model only generates the real part of nucleus-nucleus potentials. Despite the simplicity of the M3Y interaction, the DFM with the M3Y interaction was quite successful in analyzing low energy HI scattering and marked a milestone in the study of microscopic theory for nucleus-nucleus interaction in those days.
However, it was soon recognized that low energy HI scattering was only sensitive to the surface part of the potential in most cases. The "success" of the M3Y folding model thus meant that it only gave a reasonable strength of the real part of HI potentials around the spatial region where two nuclei 12/25 PTEP 2016, 06A106 Y. Sakuragi just touch their surfaces, keeping a large relative distance between their centers of mass. This implies that no meaningful information about the strength and shape of the potential at short distances was obtained by the limited success of the M3Y folding model.
Since then, it took 30 years to mark the next really important milestone in the study of microscopic theory for predicting reliable complex nucleus-nucleus interaction, which was marked by Furumoto, Yamamoto, and the present author in 2008-2009 [13,23,37]. The new interaction model has been quite successful in reproducing almost all existing data for HI scattering at incident energies per nucleon of E/A ≈ 30-200 MeV [23,38] and can be applied to any HI scattering system over a wide range of incident energies up to E/A ≈ 400 MeV and any combinations of projectile and target nuclei in the nuclear chart as long as the nucleon density is available. For the sake of further global use of the successful interaction model, they have proposed a global optical potential based on the microscopic folding model that is applicable to any combinations of nucleus-nucleus systems at any incident energy within the range of E/A = 30-400 MeV [38].
What is more important is that they unveiled the decisive roles of the many-body forces, particularly of the three-body force of repulsive nature. The many-body forces play important roles both in predicting the strength and the shape of the nucleus-nucleus potential at short relative distances and in determining the stiffness of nuclear matter, namely the shape of the saturation curve, particularly in the high density side beyond the saturation density (see Fig. 1).
Before presenting the latest achievements, we briefly review the history of the development of the folding model studies of HI scattering in some detail.

M3Y and DDM3Y interactions 6.1.1. M3Y interaction
The M3Y interaction is a real interaction having no imaginary part and with a very weak energy dependence. Hence, the DFM potential calculated with the M3Y interaction is a real potential and it has to be supplemented by a phenomenological imaginary potential when applied to the analyses of HI scattering experiments [33].
The M3Y interaction is composed of a sum of three-range Yukawa potentials. The strengths of the individual ranges of the Yukawa potential were fixed so that the matrix elements of the M3Y interaction reproduced the real-part strength of the original G-matrix evaluated at a certain low value of nuclear matter density, say about one third of the saturation density ρ ≈ 1 3 ρ 0 . The M3Y is the real potential because it was originally designed for application to nuclear structure studies and it was later applied to the construction of the real part of nucleus-nucleus potentials to be used in the study of HI scattering. Moreover, the exchange part of the simplest version of the M3Y interaction is approximated by a zero-range potential of δ-function form that makes the DFM calculation of the exchange part extremely simple, as in the case of the direct part.
These simplifications make the originally complicated numerical calculation of DFM potentials [Eqs. (6)- (12)] quite light and easily accessible to anyone who is not an expert in this research area. Despite the drastic simplification, the DFM with the M3Y interaction has been quite successful in applications to HI scattering at relatively lower energies, where, say, the incident energy per nucleon E/A is less than about 25-30 MeV [33]. Figure 3 shows an example of the DFM potential calculated with the M3Y interaction and its application to the analysis of elastic scattering of the 16  of the original strength, chosen to fine tune the fit to the data. The result is shown by the dashed curve in the right panel [21].
It can be seen that the M3Y interaction produces a very deep real potential, its central depth reaching to 500 MeV. The real potential obtained by the M3Y folding supplemented by a phenomenological imaginary potential (though not displayed in the figure) produces an elastic scattering cross section that reproduces the experimental data quite well except for extremely large scattering angles θ c.m. ≥ 50 • . The characteristic shape of the observed angular distribution of the experimental data, a rapid oscillation at small angles followed by a smooth fall off up to ∼ 50 • , is reproduced well by the calculation (the dashed curve). The same quality of fit for other experimental data has also been achieved in low energy HI scattering based on the DFM potential derived from the M3Y interaction [33]. The modifications of the DFM potential strength required to attain the best fits to the data are known to be of the order of ±10%, which may indicate an acceptable level of "success" of the model [33].

DDM3Y interaction
It is known that low energy HI scattering, in most cases, is dominated by the strong absorption of the incident flux to a huge number of non-elastic reaction channels, as mentioned before, and the elastic scattering cross sections at most scattering angles are not sensitive to the depth or shape of the real potential at short distances less than the grazing distance (or so-called strong absorption radius). Therefore, it is rather difficult to probe the short range part of the HI potential by low energy HI scattering.
For example, the calculated cross sections shown by the solid curve in the right panel of Fig. 3 are obtained with the potential shown by the solid curve in the left panel (labeled "DDM3Y-ZR"), which is much shallower than the potential given by the M3Y interaction (shown by the dashed curve), particularly at short distances.
The shallower potential is obtained by the so-called DDM3Y interaction [20,39], which is a density dependent (DD) version of the M3Y interaction obtained from the original M3Y interaction g(s,Ē) 14 Here,Ē ≡ E/A. The DDM3Y interaction itself is not directly derived from any realistic G-matrix, but the density dependent factor is taken to be a suitable functional form, such as The parameters C(Ē), α(Ē), β(Ē) of the density dependent factor were determined at each energy so that the volume integral of v DDM3Y (s, ρ,Ē) matches the volume integral of the real part of a suitable G-matrix for various densities ρ of the nuclear matter [20]. The local density ρ is evaluated in the prescription of FDA as ρ = ρ P + ρ T . The effect of density dependence is quite significant at small distances between the interacting nuclei, as seen by the large reduction from the dashed curve (M3Y) to the solid one (DDM3Y) in the left panel of Fig. 3.
It is now surprising that the two very different potentials give almost the same cross sections over the entire range of scattering angles where experimental data exists. A difference only appears at the largest scattering angles over θ c.m. ≈ 50 • where the solid curve keeps following the experimental data points while the dashed curve does not. It should be noted that the potential derived from the DDM3Y interaction becomes close to that derived from the M3Y around the nuclear surface region R ≥ 6 fm. This implies that the cross sections at small and medium angles up to θ c.m. ≈ 50 • are only sensitive to the tail part (R ≥ 6 fm in the present case) of the real potential. In other words, the potential depth inside the nuclear surface, say R ≤ 5 fm, is not relevant at all to generating the characteristic angular distribution at angles less than θ c.m. ≈ 50 • . In other words, only cross sections at extremely backward angles (θ c.m. ≥ 50 • in the present case) can probe the spatial region where the two potentials have the significant difference shown in the left panel of Fig. 3.
The scattering of light heavy ions at intermediate energies above E/A ≥ 50-100 MeV becomes rather transparent and displays a refractive nature generated by the strongly attractive real potential at short distances. The HI scattering at those energies becomes very sensitive to the shape and strength of the real potential at a small separation where two nuclei strongly overlap each other. In such situations, the density dependence of the effective interaction becomes important and the scattering experiments may be able to probe meaningful information about high density nuclear matter.

DFM with a complex G-matrix: CEG07 and MP interactions
The DDM3Y interaction was first applied to the nuclear rainbow phenomenon in α-particle scattering [39], being the typical refractive scattering by the effects of attractive real potentials, and later it has been widely applied to various kinds of HI scattering (see [40] and references therein) with considerable success. However, the DDM3Y interaction has only the real part, and the generated DFM potential also becomes real. Hence, the real DFM potential has to be supplemented by a phenomenological imaginary potential in practical analyses of HI scattering experiments, as in the case of M3Y. Moreover, another serious problem is that the density dependence of the DDM3Y interaction is just introduced by hand as a density dependent factor f (ρ,Ē), which is simply multiplied by the original M3Y interaction as in Eqs. (16) and (17) and has no fundamental theoretical background, particularly in the high density domain above the saturation density. For a correct understanding of the nucleus-nucleus interactions to be used in nuclear reaction studies from the microscopic viewpoint, it is strongly expected to construct a reliable complex G-matrix having a density dependence 15/25 PTEP 2016, 06A106 Y. Sakuragi that properly represents the medium effects in nuclear matter for nucleon densities up to twice the saturation density.
Some attempts to construct a complex G-matrix to be used in nuclear reaction studies were made in the late 1970s to early 1980s by several groups. A typical example is the work by the Oxford group led by Brieva and Rook [14][15][16], and another is the work by the Miyazaki group led by Nagata [17,18]. Individual groups constructed their original complex and density dependent G-matrix effective NN interactions to be used in calculating the complex optical potential for nucleon-nucleus elastic scattering by the folding model. They were successfully applied to proton-nucleus elastic scattering at low to intermediate energies. Those effective interactions, however, were only applicable to nucleon-nucleus systems because they were evaluated at nuclear matter density only up to the saturation value ρ 0 and, hence, could not be applicable to nucleus-nucleus systems. Another problem was that those effective interactions were generated by very early prototype NN interaction models in free space, such as the Hamada-Johnston type potential [41]. This resulted in too-small binding energies in nuclear matter of less than one half of the empirical value ofĒ ∼ = −16 MeV in the BHF calculation with those G-matrices.
To overcome those problems in the pioneering works, an ideal complex G-matrix has recently been proposed by Furumoto, Yamamoto, and the present author. They have marked another important milestone and have established one of the most reliable microscopic interaction models to date for nucleus-nucleus reaction systems on a firm theoretical base starting from the free space NN interaction. It is a complex, density dependent G-matrix evaluated by the latest elaborated NN interaction model in free space. The G-matrix is applicable to a local density up to twice the saturation density and to HI reactions at energies per nucleon E/A in the range 30-400 MeV.
The new complex G-matrix is named CEG07 [42] (complex effective potential with Gaussian form factors 2007), and was derived from a modern NN interaction model in free space called ESC04 (extended soft-core model 2004) [43,44]. This is an improved version of the original CEG83 and CEG86 interactions proposed 30 years ago by the Miyazaki group [17,18]. ESC04 is designed to give a consistent description of the interactions not only for the NN system but also for nucleonhyperon (NY) and hyperon-hyperon (YY) systems in the SU3 framework. Therefore, the theoretical framework for CEG07 can easily be extended to calculate the effective interactions for NY and YY systems in nucleon-hyperon mixed matter, such as those expected in neutron stars. It can also be applied to the study of hyperon-nucleus interactions both in the bound state (hypernucleus) and in the scattering states [45].

Roles of the three-body force in nucleus-nucleus scattering
As discussed above, the empirical saturation point cannot be reproduced by the lowest order Brueckner theory (the G-matrix approximation), as long as one only uses the two-body NN interaction, whatever the two-body NN interaction is. This is also true in the case of the CEG07 G-matrix generated by the two-body force alone (which is called as CEG07a type). This deficiency can clearly be corrected by introducing the three-nucleon force (3NF) composed of the three-body attraction (TBA) and the three-body repulsion (TBR) [46][47][48]. The TBA is typically due to the two-pion exchange with excitation of an intermediate resonance, which is the famous Fujita-Miyazawa diagram [49]. An effective two-body interaction was derived from the TBA according to the formalism given in Ref. [50] and was added to the CEG07 G-matrix. The attractive 3NF effect is known to give an important contribution at low densities. However, the contribution of the TBA alone does not resolve the problem of the saturation point. 16 The role of the TBR is far more important than that of the TBA [13]. The TBR contribution becomes more and more remarkable as the density becomes higher, which plays a decisive role for the shape of the saturation curve. It is well known that such a TBR effect is indispensable to obtain the stiff EOS of neutron star matter ensuring the observed maximum mass of neutron stars. However, the origin of the TBR is not necessarily established. In the ESC approach, the TBR-like effects are represented rather phenomenologically as the density dependent two-body interactions induced by changing the vector meson masses M V in nuclear matter according to M V (ρ) = M V exp(−α V ρ) with the parameter α V . As mentioned in Ref. [44], this TBR-like contribution introduced in ESC is found to be very similar to that of the phenomenological TBR proposed in Ref. [46] given as the density dependent two-body repulsion. The CEG07 interaction supplemented with the attractive (TBA) and repulsive (TBR) three-body forces are called CEG07b and CEG07c interactions, respectively. (The detailed differences between CEG07b and CEG07c are described in Refs. [13,23].) The CEG07 models were first successfully applied to the analysis of proton-nucleus elastic scattering over a wide range of incident energies and target nuclei [13], and later to analyses of elastic scattering of nucleus-nucleus systems with great success [23,37,38,51]. They were also applied to analyses of nucleus-nucleus inelastic scattering [52][53][54][55] with the microscopic coupled-channel formalism in which not only the diagonal potentials but also the coupling potentials that induce inelastic excitations of the colliding nuclei were obtained with the CEG07 G-matrix with the 3NF effects included. respectively. The dotted curves in Fig. 4 show the DFM potential calculated with the CEG07a Gmatrix that contains no 3NF effect, while the dashed and solid curves are the results with CEG07b and CEG07c, both inluding the 3NF effects. It can clearly be seen that the 3NF strongly reduces the real potential strength in the spatial range of middle and short separations between the interacting nuclei where the local density defined by FDA (i.e. ρ = ρ P + ρ T ) exceeds the saturation density ρ 0 . The 3NF effect on the imaginary part of the DFM potential is not so significant. Figure 5 shows the elastic scattering cross sections calculated with the use of the DFM potentials shown in Fig. 4. The solid (CEG07b) and dashed (CEG07c) curves are the results with the 3NF effects that precisely reproduce the experimental data over the whole angular range. The large deviation of the dotted curves (CEG07a), which include no 3NF effect, clearly indicates the important role of the 3NF effect.
It should be noted that, for each scattering system, the imaginary part of the DFM potential is multiplied by a common renormalization factor presented in the figures (N W = 0.7 for the 16 O + 16 O system and N W = 0.5 for 12 C + 12 C) in calculating the cross section, where the N W value is chosen to optimize the fit of the solid curve to the experimental data. The same N W is used for all three kinds of calculations in each scattering system. One might suspect that the large difference between the dotted curve and the other two curves could be compensated by using a larger N W value. In fact, the use of larger N W values drastically reduces the cross section at backward angles. However, it gives completely different angular slopes from those of the other two calculations with 3NF effects, as shown in Fig. 6 [23]. Thus, the large difference between the dotted curve and the other two curves really reflects the large difference in the real part of the DFM potential that originates from the 3NF effects.
These results clearly indicate that the elastic scattering of these light HI systems at intermediate energy is quite transparent enough to probe the difference in the depth and shape of the real potential at short distances, where the local density evaluated by the FDA can exceed the saturation value, 18 The figure is taken from Ref. [23]. and the medium effect, particularly the 3NF effect, plays an important role. In other words, it will be possible to fix some unknown parameter values relevant to the medium effects, such as the strength of 3NF, at high densities through theoretical analyses of the experimental data of HI elastic scattering based on the present microscopic interaction model.

Multi-pomeron repulsion and neutron star mass
The main subject of this last section is to introduce the scenario recently proposed by Yamamoto et al. [56], who have connected nucleus-nucleus scattering experiments in laboratories on Earth with the maximum mass of neutron stars in the Universe obtained from astronomical observations via the MPP interaction model. We will briefly review their recent studies investigating whether or not the maximum mass of 2M could be obtained from the EOS for hyperon-mixed neutron star matter, when the effects of three-body and four-body repulsion among the nucleons and hyperons are taken into account. It has been made clear that elastic scattering of light HI systems such as the 16 O + 16 O and 12 C + 12 C systems at intermediate energies around E/A ≈ 100 MeV is very transparent and quite sensitive to the strength and shape of the real potential at short distances. The DFM potential with GEG07b and the 3NF effect reproduces well the experimental data of elastic scattering cross sections. The effect of 3NF, particularly its repulsive component, plays a crucial role in the proper evaluation of the interaction potential between finite nuclei at short distances where the local density exceeds the saturation value ρ 0 or reaches up to twice the saturation value ρ ∼ = 2ρ 0 . This implies that the elastic scattering can probe the medium effect of the strong interaction in dense nuclear matter. This may open up a new window for studying the properties of infinite nuclear matter as well as neutron star matter through theoretical analyses of nucleus-nucleus elastic scattering in laboratories on the basis of Brueckner's G-matrix theory and the BHF framework. 19 As mentioned in Sect. 1, the recent astronomical observations of the maximum mass for a neutron star of twice the solar mass, (1.97 ± 0.04)M [3] and (2.01 ± 0.04)M [4], impose a severe constraint on the stiffness of the EOS of neutron star matter. This suggests the existence of strong repulsive 3NF to sustain such a large neutron star mass, since it is known that the maximum mass evaluated by any BHF calculations with only the NN two-nucleon force (2NF) cannot exceed ∼ 1.5M . On the other hand, hyperon (Y) mixing in neutron star matter, which is expected to be switched on around a neutron density of ρ ≈ (2 ∼ 3)ρ 0 , is known to bring about a remarkable softening of the EOS, which cancels the repulsive 3NF effect for the maximum mass [57][58][59][60]. One of the ideas to avoid this serious problem, called the "hyperon puzzle in neutron stars," is to consider that the 3NF-like repulsions work universally for YNN, YYN, and YYY, as well as for NNN as proposed by Nishizaki, Yamamoto, and Takatsuka [59].
The repulsive component of the 3NF adopted in the CEG07b and CEG07c type G-matrices was introduced rather phenomenologically in the form of density dependent effective 2NF, and its parameter was determined so that the saturation curve reproduced the empirical saturation point. Such a phenomenological model for the repulsive 3NF is not necessarily suitable for being extended to be used as the universal three-baryon force (3BF) among nucleons and hyperons. As a possible candidate for such a universal repulsion among multi-baryon systems, the multi-pomeron exchange potential (MPP) [61] was introduced as a model of universal repulsions between three and four baryons on the basis of the extended soft core (ESC) baryon-baryon interaction model, that is, the bare NN interaction model in free space used in the construction of CEG07-type G-matrices.

Pomeron exchange model for many-body force
The pomeron is a virtual particle that simulates the multi-gluon exchange among baryons and is thought of as being related to an even number of gluon exchanges. It is a kind of substantial particle closely related to the strong repulsive core in baryon-baryon interactions.
Recently, the multi-pomeron exchange potential (MPP) model has been proposed by Rijken, Yamamoto, and their collaborators [56,61,62], and represents the repulsive component of threebody and four-body potentials among three and four baryons including nucleons and hyperons. In the MPP model, the three-and four-body local potentials are derived from the triple-and quarticpomeron vertexes. In order to be implemented in the BHF calculation for the EOS of baryonic matter, the many-body forces are represented by the density (ρ) dependent two-body potential in a baryonic medium. It is obtained by integrating over coordinates of the third (and fourth) particles in the three-body (and four-body) potentials as follows [62]: where V (x 1 , x 2 , . . . , x N ) is the N -body local potential by pomeron exchange. In the case of N = 3 (three-body) and N = 4 (four-body), the effective ρ dependent potentials have the following forms: V (4) eff (r) = g (4)  Here, m P and g N P are the pomeron mass and the two-body pomeron strength, which are taken to be the same as those in the ESC model for 2NF, and the scale mass M is taken as a pomeron mass.
The MPP strengths g P and g (4) P are unknown parameters that will be determined so that the experimental data are precisely reproduced with the DFM potential derived from the G-matrix including MPP. In other words, the G-matrix folding potentials including MPP contributions are used to analyze 16  P and g (4) P ) are adjusted so as to reproduce the experimental data. Figure 7 shows the results. In these calculations, the bare NN interaction ESC04 in free space has been replaced by the latest improved version called ESC08. The dot-dashed curve is the result with the 2NF only and largely deviates from the experimental data, while the other three curves include the many-body forces based on the MPP model. The dotted curve (MPb) includes up to the 3NF, while the solid (MPa + ) and dashed (MPa) curves include both the 3NF and 4NF contributions. In the latter three curves, the MPP strengths g (3) P and g (4) P are determined so as to reproduce the data [56]. It should be noted that, in contrast with the DFM analyses with the CEG07 interactions, no renormalization is required to the imaginary part of the DFM potentials (i.e., N W is fixed at 1.0) to attain satisfactory fits to the data shown in the figure.

Stiffness of neutron star matter fixed by nucleus-nucleus scattering
The MPP interactions, the strength of which were determined from analysis of nucleus-nucleus scattering, are then included in constructing the EOS of neutron star matter as well as symmetric nuclear matter, to see whether they result in an EOS stiff enough to give the observed large neutron star mass [3,4]. It should be noted that the MPP is defined so as to work universally, not only in NNN states but also YNN, YYN, and YYY states, with the common parameters fixed by analysis of nucleus-nucleus scattering. Therefore, there is no additional free parameter for the interaction model. Corresponding to the determined MPP, the attractive component of 3NF (TNA) is added to 21  reproduce the nuclear saturation property precisely. Since TNA is important mainly in the lower density region, the addition of TNA does not lead to any additional ambiguity in the following discussion about the stiffness of hadronic matter in high density regions. Figure 8 shows the results obtained for the saturation curves for symmetric nuclear matter as well as for pure neutron matter [56]. The upper curves are for neutron matter and the lower curves are for symmetric nuclear matter. The solid, dashed, and dotted curves are the results obtained with the MPa, MPa + , and MPb interactions, respectively. The box shows the empirical value of the saturation point. The inset shows a zoom of the region around the saturation point for symmetric nuclear matter. It can be seen that the inclusion of the quartic-pomeron exchange contributions (in MPa and MPa + ) in addition to the triple-pomeron exchange stiffens the saturation curve in the high density region in both symmetric nuclear matter and neutron matter. This is because the strengths of the effective two-body interaction derived from quartic-pomeron exchanges are proportional to ρ 2 [see Eq. (21)] and the contribution becomes sizable in the high density region.
The final step is to study properties of neutron stars with and without hyperon mixing on the basis of these multi-baryon interaction models.
Assuming mixed matter of n, p, e − , and μ − in chemical equilibrium, the TOV equation is solved for the hydrostatic structure of a spherical non-rotating star [56]. The obtained mass-radius relations of neutron stars without hyperon mixing taken into account are demonstrated in the left panel of Fig. 9. The solid, dashed, and dotted curves are for MPa, MPa + , and MPb, respectively. It can be seen that the EOSs in these cases are found to be stiff enough to give 2M . The above differences appear significantly in the inner regions of 16 O + 16 O DFM potentials in the inner region, though they cannot be seen in the cross sections in Fig. 7. It should be noted that the calculation with 2NF alone leads to EOS too soft to give 2M , though it is not shown here.
However, the hyperon mixing will be set on when the neutron density becomes as high as (2 ∼ 3) × M , which will be expected to soften the EOS substantially. The EOS of β-stable neutron star matter composed of neutrons (n), protons (p + ), electrons (e − ), muons (μ − ), and hyperons ( and − ) is derived from the G-matrix calculation with the present interaction model. Using the EOS of hyperonic neutron star matter, the Tolmann-Oppenheimer-Volkoff (TOV) equation for the 22  hydrostatic structure is solved and the mass-radius relations of neutron stars are obtained [56]. The results are shown in the right panel of Fig. 9.
The maximum mass values calculated for MPa + , MPa, and MPb are 2.07M , 1.94M , and 1.83M , respectively, being smaller by 0.61M , 0.51M , and 0.35M than the values without hyperon mixing. Although the universal repulsion works to raise the maximum mass, the hyperon mixing is also enhanced by it so that the maximum mass is reduced.
The maximum mass for MPb is considerably smaller than the observed value of ∼ 2M . On the other hand, those for MPa and MPa+ reach this value owing to the four-body repulsive contributions.

Concluding remarks
Nuclear physics provide a unique playground to survey the properties of the dense hadronic matter in the Universe, including that in neutron stars. Precise measurements of refractive nucleus-nucleus scattering in laboratories on Earth give us valuable information about high density nuclear matter and provide a useful test ground of the effective interactions in high density baryonic matter that act not only between two nucleons (NN) but also among triple and quartic baryons (BBB). The key bridge to connect the nuclear experiments in laboratories and the baryonic matter in the Universe would be a modern theory of the effective interaction in many-fermion systems. The new Brueckner G-matrix constructed with the latest models for bare NN interaction in free space together with the recently developed many-body forces effective in high density baryonic matter is a promising candidate for this bridge.
I should like to dedicate the present review article to the memory of the late Professor Yoichiro Nambu.