Probing neutronproton dynamics by pions
Abstract
In order to investigate the nuclear symmetry energy at high density, we study the pion production in central collisions of neutronrich nuclei at 300 MeV/nucleon using a new approach by combining the antisymmetrized molecular dynamics (AMD) and a hadronic cascade model (JAM). The dynamics of neutrons and protons is solved by AMD, and then pions and resonances in the reaction process are handled by JAM. We see the mechanism how the resonance and pions are produced reflecting the dynamics of neutrons and protons. We also investigate the impacts of cluster correlations as well as of the highdensity symmetry energy on the nucleon dynamics and consequently on the pion ratio. We find that the production ratio agrees very well with the neutronproton squared ratio in the highdensity and highmomentum region. We show quantitatively that production ratio, and therefore , are directly reflected in the ratio, with modification in the final stage of the reaction.
pacs:
21.65.Ef, 25.70.z, 25.80.LsI Introduction
It is one of the important subjects in nuclear physics and astrophysics to determine the nuclear symmetry energy at various densities. Constraints on the symmetry energy at densities lower than the normal density have been obtained to some degree from nuclear physics experiments Tsang:2012 ; Horowitz:2014 . On the other hand, intermediateenergy heavyion collisions are believed to be useful to investigate the symmetry energy at high densities. In fact, in central collisions at several hundred MeV/nucleon, transport calculations show that the maximum density around is reached when the system is compressed in an early stage of the reaction. The ratio of the neutron and proton densities in the compressed part is naturally sensitive to the symmetry energy at high densities BaoAnLi:2008 ; Li . However, it is necessary to understand the link between the effects in the compression stage and the final observables, in order to extract the highdensity symmetry energy from experimental data available at present and in future.
The ratio has been proposed to be a good probe to constrain the highdensity behavior of the symmetry energy Li . In heavyion collisions, the pions are produced through the resonance formation in the nucleonnucleon collisions that typically occur at early times in the compressed part of the system. If all the pions produced by the process are directly emitted, the expected ratio should be , where and are the numbers of neutrons and protons that are relevant for the production. In another case of chemical equilibrium for and reactions at a temperature , the ratio Bertsch:1980 ; Bonasera:1987 , where and are the neutron and proton chemical potentials, may also be related to the neutron and proton densities or the symmetry energy Li ; Ferini:2005 . We also notice the relation when the neutron and proton phasespace densities are compared at the same singleparticle energy satisfying . Thus, both in these extreme cases, the ratio is related to some kind of neutrontoproton squared ratio which is then supposed to be sensitive to the symmetry energy at high densities.
Some theoretical studies have been performed by different transport models to investigate the sensitivity of pion observables Li ; Gaitanos:2004 ; Reisdorf:2007 ; IBUU04 ; ImIQMD ; pBUU ; BUU . At present, however, some of these results are contradicting to each other even qualitatively. The symmetryenergy dependence of the ratio in Ref. ImIQMD is opposite to that predicted by Refs. Li ; IBUU04 , while Ref. pBUU predicts that the ratio for the total pion multiplicities does not depend on the density dependence of the symmetry energy. Moreover, the ratio in central Au + Au collisions measured by the FOPI Collaboration Reisdorf:2007 ; FOPI is actually larger than the squared neutrontoproton ratio of the total system, , at low energies such as at 400 MeV/nucleon. Since the compressed part naturally becomes less neutron rich than the total system, we inevitably have , and therefore cannot agree with , contradicting to the arguments in the previous paragraph, as long as is identified with the density ratio in the compressed part as is typically done in the literature. Most of the transport calculations underestimate the measured pion ratio, but some recent calculations ImIQMD ; pBUU predict the ratio as high as the experimental data. Nevertheless, investigations cannot be found in the literature on the mechanism to increase the pion ratio. Hence, our understanding of the pion ratio in relation to the symmetry energy and the nucleon dynamics is not complete particularly when we need to have a precise description to explain the absolute magnitude of .
In this paper, we study the pion production in central collisions of neutronrich nuclei at 300 MeV/nucleon, which is one of the systems to be measured at RIBF/RIKEN RIBF . We develop and employ a new approach by combining the antisymmetrized molecular dynamics (AMD) AMD and a hadronic cascade model (JAM) JAM . The dynamics of neutrons and protons is solved by AMD, and then pions and resonances in the reaction process are handled by JAM. AMD calculations were performed for several cases with and without cluster correlations Ono:NN2012 , and with two effective interactions corresponding to different density dependences of the symmetry energy. These different cases of AMD calculation yield different dynamics of neutrons and protons. The main aim here is to utilize these different cases to identify how the resonances and the emitted pions carry the information of the nucleon dynamics. We will also see the impacts of cluster correlations as well as of the highdensity symmetry energy on the nucleon dynamics and consequently on the pion ratio.
Ii Formulation
In order to understand the enhancement mechanism and to extract highdensity symmetry energy, we need a reliable transport model of nucleons, clusters, resonances and pions. AMD has been demonstrated to be a reliable transport model of nucleon and clusters AMD ; Ono:NN2012 . It takes account of nucleon mean field effects, twonucleon elastic scatterings, eventbyevent fluctuations and antisymmetrization of nucleon wave functions. All of these ingredients could affect the ratio during the collision, and thus have influences on the ratio. One can explain fragment mass distribution in heavyion collisions precisely, especially with cluster correlations. In the present AMD code, however, resonances and pions have not been incorporated. On the other hand, JAM is a reliable hadron transport model JAM . It has been successfully applied to and collisions in the energy range from 1 GeV/nucleon to 158 GeV/nucleon. We can describe hadron production spectra in JAM, and collective flows are also explained with the mean field effects switched on Isse:2005nk . By comparison, antisymmetrization of nucleon wave functions is not included, and cluster correlation during the time evolution are not implemented. At present, there is no transport model which includes all of the above mentioned important ingredients; mean field effects, twobody collisions including elastic and inelastic processes, eventbyevent fluctuations, antisymmetrization of nucleon wave functions, and dynamical cluster correlations. Thus one of the best ways would be to combine AMD and JAM — nucleon transport in AMD and particle production in JAM.
ii.1 Perturbative treatment of pions and
We consider here a transport model which describes the reaction dynamics by the time evolution of onebody density matrices or corresponding phasespace distributions , where the index stands for the particle species such as nucleons (), resonances () and pions (). It should be implicitly understood that the isospin (and spin) components of these particles are also distinguished by . The coupled equations for () can be written in general as
(1) 
with the meanfield term and the collision term , both of which depend on the phasespace distributions of all the particles at the time in general. The collision term includes at least the following processes
(2)  
(3)  
(4)  
(5)  
(6) 
Now we limit our consideration to the subthreshold or nearthreshold cases where only a small number of production processes occur because the incident energy is not so high. Then we expect that the production can be treated perturbatively. If the perturbation parameter is multiplied to the probability of the process, resonances and pions will appear in the first order of , and . The zeroth order equation for the nucleon distribution is
(7) 
where the collision term includes only the elastic scatterings. The distributions of resonances and pions can be solved by
(8)  
(9) 
which are correct up to the first order of .
The equation (7) for the zeroth order of the nucleon distribution can be solved assuming a system composed of only nucleons without considering productions of other particles. In our present approach, we solve the nucleon dynamics by AMD AMD ; Ono:NN2012 . Then, for the calculated , the equations (8) and (9) for resonances and pions are solved by another transport model JAM JAM which can handle particle productions. The information of nucleons in the JAM calculation is always replaced by calculated by AMD. Namely the particle production is calculated by JAM based on the nucleon dynamics calculated by AMD.
The above treatment violates some conservation laws in the higher orders . The result may be improved by introducing corrections for the conservation laws of baryon number, change and energy, modifying the nucleon information in Eqs. (8) and (9). When a resonance exists in the JAM calculation, it should have been produced by a collision of two nucleons. Such a pair of two nucleons is chosen in the AMD calculation by taking into account the distance from the resonance, the phase space to produce a resonance, and the charge conservation condition. Then one of the two nucleons is annihilated (assuming that it is replaced by the resonance) and the charge and the momentum of the other nucleon is modified for the charge and energy conservations. In the pion case, we consider only the charge conservation by modifying the charge of a nucleon for each pion if necessary. We have checked the validity of this prescription by comparing two different calculations performed by the JAM code. The first calculation is done by the JAM+JAM calculation which is the same as the AMD+JAM calculation described above but Eq. (7) is solved also by JAM by turning off all the inelastic collisions. The result is compared with the standard JAM calculation which solves all the particles as usual. With the present prescription for the conservation laws, the pion multiplicities and the pion ratios in the two calculations agree well within the errors of about 10% and 2%, respectively, at the incident energy 300 MeV/nucleon. These agreements get worse slightly at the incident energy 400 MeV/nucleon.
The incompleteness of this prescription is the dominant origin of the violation of the energy conservation in the AMD+JAM calculation, which can be estimated by the JAM+JAM calculation. It turned out that the total energy per baryon is higher than the initial value by about 2 MeV on average at fm/ in the collisions at 300 MeV/nucleon. This average increase of the total energy seems roughly consistent with the abovementioned 10%overestimation of the pion multiplicity. In fact, the pion multiplicity increases by 13% in the standard JAM calculation when the incident energy is raised from 300 to 310 MeV/nucleon, corresponding to the 2.3MeV increase of the total energy per baryon.
ii.2 Amd
AMD AMD describes the dynamics of a manynucleon system by the time evolution of a Slater determinant of Gaussian wave packets
(10) 
where the wave packet centroid is denoted by which is a complex vector, and the spinisospin state takes , , and . The width parameter is chosen to be as usual. The corresponding phasespace distribution is
(11) 
with , and .
ii.2.1 Mean field term
The mean field term is given by the equation of motion for the wave packet centroids derived from the timedependent variational principle.
In the present calculations, we employ the Skyrme SLy4 force SLy4 as the effective interaction with the spinorbit term omitted. The corresponding nuclearmatter incompressibility is MeV at the saturation density . The nuclearmatter symmetry energy at is MeV with the slope parameter MeV (called ‘asysoft’ or soft symmetry energy in Sec. III). In order to study the effect of the density dependence of the symmetry energy, we also perform calculations with a force obtained by changing the density dependent term in the SLy4 force
(12) 
to
(13) 
By choosing , we have a force corresponding to MeV (callded ‘asystiff’ or stiff symmetry energy) with the same equation of state of symmetric nuclear matter and with the same as the original SLy4 force.
The Skyrmetype prarametrization of the effective interaction is advantageous for the efficient AMD computation. However, the Skyrme forces have a quadratic momentum dependence of the mean field which is not valid at high energy collisions at several hundred MeV/nucleon. Therefore, the momentum dependence is corrected for the present calculations in a similar way to Ref. GBDG in BUUtype calculations. The detailed formulation in the case of AMD is given in the Appendix A.
ii.2.2 Collision term with and without clusters
The twonucleon collision process corresponds to the collision term . In AMD, a twonucleon collision is treated as a stochastic transition from an AMD state to another AMD state specified by the relative momentum between the scattered two nucleons . The transition rate is expressed as
(14) 
In general, medium modification is introduced for the scattering matrix elements. However, in the present calculations, we employ the matrix elements in the free space. It should be noted that some medium effect still exists in the dependence of the final state energy in a similar way to a BUU calculation of Ref. BUU . The Pauli blocking for the scattered nucleons is taken into account.
In the usual treatment of twonucleon collisions, only the states of the scattered two nucleons are changed in the final state (see Ref. AMD for the precise description of the method which employs ‘physical coordinates’). On the other hand, an extension has been introduced Ono:NN2012 to allow direct formation of light clusters with and 4 in the final state . Namely, in the calculation with cluster correlations, when two nucleons and collide, we consider the process
(15) 
in which each of the scattered nucleons () may form a cluster with a spectator particle . This process includes the collisions without cluster formation as the special case of with empty . The transition rate of the clusterforming process is given by Eq. (14) with the suitable choice of the final state . When a cluster is formed, the corresponding wave packets are placed at the same phasespace point, i.e., the cluster internal state is represented by the harmonicoscillator configuration. Denoting the initial and final states of the system by and , respectively, we have the transition rate
(16) 
where are the states after the momentum transfer to the nucleons (), and is the relative momentum between and in these states. The matrix element is the same as for the usual twonucleon collisions. We use an average value of evaluated at and that evaluated at the initial relative momentum.
The actual situation of a twonucleon collision requires more considerations because there are many possible ways of forming a cluster for each of the scattered nucleons and . For a scattered nucleon , we first consider the possibility that may form a cluster with one of the nucleons which have the same spinisospin state. This spinisospin state that is studied first is randomly decided. The clusterformed state is denoted by which is obtained, by first changing the state to by the momentum transfer to , and then moving the two wave packets of and to the same phasespace point without changing their center of mass. Since the different final states are not orthogonal , the probability that forms a cluster with one of should be calculated as
(17)  
(18) 
This probability is calculated with an approximation that the manybody state is a direct product of wave packets centered at the ‘physical coordinates’ AMD . With the calculated probability , a cluster will be formed with one of . It is somewhat arbitrary which one of should be chosen with what probability. In the present calculation, we choose with the relative weight with the parameter . With the rest of the probability (), does not form a cluster with a nucleon of this spinisospin state. The procedure is repeated for other spinisospin states for . The particle should be regarded as a cluster, instead of a scattered nucleon, if a (sub)cluster has been already formed in previous steps of the repetition. Thus the formation of light clusters is considered up to an particle. This procedure determines the probability for the combination of final clusters as a function of the momentum transfer or . It satisfies the normalization . The factor in Eq. (16) should be replaced by .
Even when the cluster formation is introduced, the manybody state is always represented by an AMD wave function which is a Slater determinant of nucleon wave packets. The time evolution of the manybody state is solved just as usual without depending on whether some of the wave packets form clusters due to collisions in the past (except for the clustercluster binding process in the next paragraph). This is in contrast to BUU by Danielewicz et al. Danielewicz:1991 where clusters are treated as new particle species. In our approach, a nucleon in a formed cluster may collide with some other nucleon so that the cluster is broken. It may be the case that the scattered nucleon forms the same cluster as before, so that an elastic scattering of the cluster is possible. All of these kinds of processes are based on the nucleonnucleon scattering matrix elements, without introducing parameters to control individual channels of cluster formation. In the present calculations with cluster correlations, however, the overall cluster production probability is suppressed, when the momentum transfer is extremely small, by a factor .
It has been found that we should take account of the correlations to form heavier fragments via coalescence of light clusters on top of the usual time evolution of AMD Ono:NN2012 . This option of improvement has been turned on in the present calculations even though it will not strongly influence the following discussions in high energy collisions. The details are described in Appendix B for completeness.
It is experimentally clearly known that the clusters are important in many situations of heavyion collisions even though the incident energy is relatively high. For example, the FOPI data FOPI show that only 21% of the total protons in the Au + Au system are emitted as free protons in central collisions at 250 MeV/nucleon, and all the other protons are bound in light clusters and heavier fragments. AMD without cluster correlations overestimates the proton multiplicity as many other transport models do. On the other hand, the AMD calculation with cluster correlations can reproduce this feature very well as shown in Ref. Ono:NN2012 for systems including the Au + Au central collisions at several hundred MeV/nucleon. Therefore, the calculation with clusters is believed to be much closer to the realistic case than without clusters. In the following section, however, we are going to show the results of both calculations with and without clusters, which are useful for the purpose to study the dependence of the pion production on the nucleon dynamics.
ii.2.3 Test particles
The information of nucleon dynamics calculated by AMD is sent to the JAM calculation in the form of a set of test particles , ,…,. One test particle per nucleon is generated following the distribution function defined by Eq. (11). It should be noted that all kinds of quantum effects from the antisymmetrization of the manybody state are contained in this distribution function (or the Wigner function). In particular, it is not positive definite, and therefore in the phasespace region of the probability has to be replaced by zero, which can potentially introduce some inaccuracy of the testparticle representation. To check the accuracy, we compared the density profile for the ground state of the Au nucleus, to find no visible difference between the distribution of the generated test particles and the exact density profile. Therefore, this method of test particles should be sufficiently accurate in highly excited situations during heavyion collisions. The method to generate the test particles is described in Appendix C.
The set of test particles is sent to the JAM calculation at every 2 fm/. We have checked that the result does not change when it is sent at every 1 fm/.
The AMD calculation is much more time consuming than other transport models such as JAM. We typically generated 1000 AMD events for each case of the present calculations. However, we improve the statistics for the pion production by generating 500 JAM events from the same AMD event. As we will see in the next section, a sufficient statistical accuracy is obtained for the pion production with this limited number of AMD events.
ii.3 Jam
JAM is a transport model which is developed by Nara et al. JAM . This model has been successfully applied to highenergy collisions up to more than one hundred GeV/nucleon. In this model, in the energy domain relevant for the present work, the hadronhadron reactions are treated by the cross sections based on experimental data and the detailed balance. In particular, in the present work, isospin symmetry is assumed. The cross section for the process at the c.m. energy is written as Teis:1996kx ; Weil:2012ji
(19) 
where and are the initial and final momenta in the c.m. frame, and the ClebshGordan factor is or . The matrix element is assumed as
(20) 
with GeV and GeV. This is a similar parametrization to UrQMD in Ref. UrQMD . The constant mb GeV is determined to fit the data of the cross sections. The distribution of the mass of the resonance is determined by the integrand of Eq. (19). is the decay width for parametrized as in Refs. JAM ; UrQMD . In order to fit the data precisely near the pion threshold, nonresonant contributions and/or some components that violates isospin symmetry are necessary for the processes, which are ignored, however, in the present work so that the pions are produced only through the formation of resonances. In our JAM calculation, the mean field for nucleons is not included.
To know how the JAM calculation reproduces the pion multiplicity of the experimental data, we have calculated Au + Au collisions at various incident energies for the impact parameters fm. Figure 1 shows the energy dependence of the pion multiplicity. The solid lines indicate the pion multiplicities calculated by JAM. The points show the experimental data taken from the FOPI measurement of Au + Au collisions at 0.4, 0.8 and 1.5 GeV/nucleon FOPI . We find that the JAM calculation almost reproduces the experimental data of pion multiplicities reasonably well. However, the calculation overestimates the pion multiplicities in particular at lower energies. This is probably because the JAM calculation here does not include the mean field potential, and therefore more energy is available for particle production without the cost of energy to compress the system Danielewicz:1991 ; Kruse:1985hy .
We have also checked the the time evolution of the numbers of pions and resonances in the JAM calculations as shown in Fig. 2 for the central Au + Au collisions at 1 GeV/nucleon. The result shows a similar behavior to a relativistic meanfield transport calculation reported in Fig. 1 of Ref. Ferini .
Thus JAM calculation provides a reasonable description for the pions and resonances produced in collisions at around 1 GeV/nucleon. It should be noted that in the AMD+JAM calculation the mean field for nucleons is taken into account in AMD and the ignored mean field in JAM does not actually influence the results.
The production and absorption reactions for and pions occur in the JAM calculation as in the free space without medium modification for the thresholds, while nucleons feel potential in the AMD calculation. This corresponds to assuming that the potentials and for and pions are related to the isospindependent nucleon potential as
(21) 
where is the isospin component. This is equivalent to the choice in the pBUU calculation of Ref. pBUU if the momentum dependence is ignored. It should be noted that other transport calculations take a different choice Li ; Gaitanos:2004 . The details of the inmedium effects for and pions may influence the pion yields as investigated in equilibrium calculations JunXu:2010 ; JunXu:2013 . On the other hand, recently, the effects of the completely unknown symmetry (isovector) potential of the resonance on the pions in heavyion collisions have been studied in Ref Li:2015hfa . It has been reported that these effects are negligible except for at deeply subthreshold energies and thus the pion ratio is still a good probe to investigate the highdensity symmetry energy.
Iii Results
We calculated collisions of at 300 MeV/nucleon for the impact parameters fm. At the initial time , the centers of the two nuclei are separated by 15 fm. In order to investigate the relation between the highdensity symmetry energy, the nucleon dynamics in a compressed neutronrich system and the emitted pions, we are going to compare the results from five different cases as follows,

AMD + JAM with clusters (asysoft)

AMD + JAM with clusters (asystiff)

AMD + JAM without clusters (asysoft)

AMD + JAM without clusters (asystiff)

JAM (no mean field).
The first two cases are calculated with AMD with cluster correlations, and the next two cases are calculated without cluster correlations. For each of them, we have calculated with two different effective interactions (‘asysoft’ and ‘asystiff’) for different density dependence of symmetry energy. We have also performed a standard JAM calculation without combining with AMD in order to clarify the effect of the mean field and the symmetry energy by the comparison with the AMD+JAM calculations.
iii.1 Neutronproton dynamics
To see the dynamics of neutrons and protons in the five different calculations, Fig. 3 shows some information on the time evolution of the densities of protons and neutrons. The upper panels show the neutron and proton densities in the radius of 2 fm from the centerofmass of the system. The maximum density is reached at fm/. However, we find that the maximum density is higher with cluster correlation than without it. The time of the maximum density also depends on the cluster correlation. One of the possible reasons for this is that the cluster formation has an effect to gather nucleons spatially so that the compression of the central part of the system continues longer. In the case of JAM calculation, the higher maximum density is reached at an earlier time than in the AMD without clusters, which is reasonable because the meanfield potential for nucleons is not included in JAM.
The lower panels of Fig. 3 show the time evolution of the neutrontoproton ratio calculated for the central region within a radius from the center of mass. The radius is determined at each time in each event by the condition , where the spherically averaged density is evaluated by using the set of test particles. The results with the soft and stiff symmetry energies are shown by the solid and dashed lines, respectively, for each case of with and without cluster correlation. In all the cases, the ratio of the compressed part (at fm/) becomes smaller than the that of the total system , which is consistent with the symmetry energy effect that does not favor highdensity neutronrich matter. We can see clearly that this effect to reduce of the compressed part is stronger with the stiff symmetry energy. This symmetry energy effect is consistent with the results of other transport models Li at least qualitatively. However, there may be model dependence in the quantitative values of . In fact, in our calculations here, the symmetry energy effect is stronger without cluster correlation than with cluster correlation.
iii.2 and pions
In Fig. 4, we show the reaction rates of production () and absorption () in the upper and middle panels, respectively. The time evolution of these rates for , , and are shown by the four lines. We can see that about 70% of the produced resonances are absorbed and turned back to nucleons. We also show the numbers of existing resonances and pions in the lower panel of Fig. 4.
In Fig. 5, we show the time evolution of the ratio for the five different cases. The particles are defined, including the resonances depending on the branching ratio to decay into the pion, as
(22)  
(23)  
(24) 
Our calculation predicts that the evolution of the pionlike ratio has a dependence on the symmetry energy. The pionlike ratio calculated with the soft symmetry energy is larger than that with the stiff symmetry energy in both of calculations with and without cluster correlations. This result seems to be similar to the predictions reported in Refs. Li ; IBUU04 qualitatively. We also find that the pionlike ratio depends on the cluster correlations and that the symmetryenergy effect appears in the pionlike ratio stronger in the case without cluster correlations. This is consistent with the result of the neutronproton dynamics as shown in the lower panels of Fig. 3, suggesting a possibility that the pion production is really related to the neutronproton dynamics.
The pionlike ratio reaches the final value at around fm/. In all the five cases, the predicted final ratios become larger than of the total system, which is consistent with the experimental observation for Au + Au system at 400 MeV/nucleon Reisdorf:2007 and is suggesting that the relation does not hold if is taken from Fig. 3. On the other hand, the behaviors of the pionlike ratio before fm/ are complicated. The origin of these behaviors will be better understood through the analysis in the next subsections.
iii.3 Relation of nucleon dynamics and production
As mentioned above, the calculated pion ratio becomes larger than the ratio of the compressed part of the system. To find the origin of this effect, we investigate what kind of information of neutrons and protons is carried by resonances. Since and are produced only by the and reactions, respectively, we expect that the ratio of the production rates of these resonances should be most directly linked to some kind of ratio of nucleons. Figure 6 shows the ratio of the production reaction rates of and as functions of time for the five different calculations. The production is peaked around 20 fm/ as shown in Fig. 4, and the ratio does not have a significant meaning at very early and late times. This production ratio is compared with Fig. 7(a) which shows the same information as in the lower panels of Fig. 3 but shows the squared ratio . We find clearly that the production ratio is much larger than this ratio in the highdensity region of the system, except for the JAM case in which both ratios are close to the ratio of the total system. The relative ordering of the ratios for the five cases also disagrees between Fig. 6 and Fig. 7(a). This result shows that the production, and therefore the pion production, are not simply linked to the ratio of the highdensity part of the system.
Each panel of Fig. 7 shows the time evolution of the squared neutrontoproton ratio calculated for the nucleons that satisfy the following condition.

Nucleons in the sphere centered at the centerofmass of the system.

Nucleons with high momenta in the sphere centered at the centerofmass of the system. We take MeV/. The collective radial momentum is subtracted from the nucleon momentum .
For the condition (a), we choose the nucleons in the highdensity central region within a radius from the center of mass. The radius is determined by the condition as described in the subsection A. For the condition (b), we choose only the nucleons with momenta larger than MeV/ in the centerofmass system, in addition to the condition (a). The collective radial momentum is subtracted for this condition, where is the radial momentum component averaged for the nucleons on the sphere of the radius . It is natural to consider this kind of momentum condition since a sufficient energy is required to excite a resonance in a twonucleon collision. Our choice of corresponds to MeV while MeV. We have also checked that the changes of by MeV/ result in differences in (for the average values shown in Fig. 8).
The panels of Fig. 7 correspond to the results of the time evolution of the ratio calculated for the nucleons satisfying the conditions (a) and (b), respectively. We find that when the high momentum condition is imposed, the ratio changes drastically compared to that without the momentum condition. The ratio in condition (b) becomes larger than that in condition (a).
By comparing the production ratio in Fig. 6 and the ratio in Fig. 7(a) for each condition, we have already seen that and do not agree if the nucleons are selected only by the highdensity condition (a). On the other hand, in the result of the condition (b), is quite similar to , i.e., the relation holds as a function of time. Thus, we can conclude that resonances, and hopefully pions, carry direct information on nucleons in high density and high momentum region of the onebody phase space. We also mention that the agreement is not as perfect as in the case (b) at 30 fm/ if the collective radial momentum is not subtracted to define the condition.
iii.4 From nucleons to pion ratios
To discuss the relation between the dynamics of nucleons and the final pion ratio, we show in Fig. 8 a summary of different ratios in five calculations.
For the timedependent nucleon ratio shown in Fig. 7, we here define a representative ratio as
(25) 
where and indicate the numbers of neutrons and protons as functions of time which satisfy the conditions described in the previous subsection. This ratio carries the information in the compression stage because and take large values only if the system has a large highdensity region . In the first two columns of Fig. 8, the and ratios show the calculated representative values when the nucleons are selected by the conditions (a) and (b), respectively. As we have already seen in the previous subsection, the increases by choosing the high momentum part of the phase space. The effect is stronger in the calculations with cluster correlations. This is understandable because, with cluster correlations, many clusters are formed which contain the same number of neutrons and protons and have relatively low momentum per nucleon, and therefore the remaining part of highmomentum nucleons becomes neutron rich. Another important point found here is that the symmetry energy effect, namely the difference between stiff and soft symmetry energies, is smaller when the cluster correlation is turned on. This is also reasonable because the cluster correlation forces some neutrons and protons to move together, and therefore the different forces acting on neutrons and protons are averaged out to some degree.
As a representative value of the production ratio, we show, in the third column of Fig. 8, the ratio of the total production numbers
(26) 
where and indicate the reaction rates of the production at each time shown in the upper panel of Fig. 4. We can see that the ratio is different from the ratio, while the ratio is almost equal to the ratio. These results are consistent with the comparison of Fig. 6 and Fig. 7.
It might not be straightforward, in principle, how the production ratio is related to the pion ratio because many of the produced resonances are absorbed by reactions as we have seen seen in Fig. 4. However, we find that the behavior of the production ratio before fm/ in Fig. 6 is similar to the pionlike ratio in Fig. 5 calculated from the resonances and pions that exist at each time. The ratio of numbers of resonances at fm/ is in the lower panel of Fig. 4 for the case of soft symmetry energy with cluster correlations. The to ratio is significantly larger than the production ratio because the neutronrichness of the system influences the isospin dependence of the absorption rates^{1}^{1}1 If the system were in chemical equilibrium with a temperature and neutron and proton chemical potentials and , one would expect the multiplicity ratio might follow for . In our simulation, however, is not as large as . This implies that dynamical effects in the neutron and proton distributions are important and/or the full chemical equilibrium for resonances is not achieved. . The pionlike ratio calculated from only these resonances is 2.86 which happens to be similar to the production ratio. The fourth column of Fig. 8 shows the pionlike ratio at fm/ (see Fig. 5). An interesting observation is that the dependence of the pionlike ratio on the symmetry energy and the cluster correlations is quite similar to that of the production ratio and therefore to that of . This suggests that the information on the highdensity nucleon dynamics remains in the pionlike ratio at fm/, without being much influenced by the absorption.
The pionlike particles have to go through the exterior region of the expanding system. The symmetryenergy effect on the nucleon ratio in the exterior region should be opposite to that in the inner region because the total numbers of neutrons and protons are (almost) conserved. As shown in the rightmost column of Fig. 8, in the final stage of the reaction, the pion ratio is modified to some degree from the pionlike ratio at fm/ to the final ratio. In each case of with and without clusters, the symmetry energy effect at fm/ is reduced in the final ratio to about 70% of the value at fm/. We also find the effect of clusters tends to raise the pion ratio in the final stage, which is probably because the interior/exterior part becomes less/more neutron rich when more clusters, with the same number of neutrons and protons contained, are formed with relatively low velocities.
iii.5 Pion spectra
Finally, we investigate the pion spectra. Fig. 9 shows the pion spectra with or without Coulomb force for charged pions. Coulomb force evidently changes the pion spectrum because is accelerated and is decelerated. The spectral ratio shown in Fig. 10 can become very large at low pion momenta due to the Coulomb effect.
The symmetryenergy effect, however, appears in our present results as a simple normalization factor for the pion spectra and the spectral ratio, and thus we cannot obtain information on the symmetry energy effect more than found in the ratio of the total pion multiplicities. This point does not agree with what has been found in the pBUU calculation pBUU ; Tsangprivate in which the symmetryenergy effect is strong in the high momentum part of the spectra.
Iv Summary
The mechanism of pion production was studied with a new approach by combining two transport models AMD and JAM. For central collisions at 300 MeV/nucleon, the production of resonances and pions are treated as perturbation. Two different AMD calculations with and without cluster correlations were performed, not only to investigate the effect of clusters but also to study the correlations between the nucleon dynamics and and pion production. We found that the production ratio agrees very well with the neutronproton squared ratio in the highdensity and highmomentum region of the onebody phase space. We also found that the production ratio, and therefore , are directly reflected in the ratio. The effect of the highdensity symmetry energy in the ratio is modified in the final stage of the reaction, with a large part of the effect still remaining, which is qualitatively similar to the case of BUU in the literature Li .
If the AMD calculations with and without clusters are regarded as two different models, the present results show the value of the pion ratio is modeldependent as well as the nucleon dynamics is. Evidently models should be constrained by other observables such as the multiplicities and the spectra of clusters in order to extract the symmetry energy from the pion ratio. It may also possible to find a combination of observables that has less model dependence, such as by taking the ratio of an observable from different reaction systems BaoAnLi:2008 . It is, nevertheless, preferable to clarify the origin of the different predictions of different models by investigating the dynamics in detail as has been done here.
It is, of course, interesting and necessary to extend the present study to other reaction systems with different neutronproton asymmetries and for different energies. Such calculations are in progress, and the results will be reported elsewhere.
Acknowledgments
This work was supported by JSPS KAKENHI Grant Numbers 24105008, 15K05079 and 15H06413.
Appendix A Momentum dependent term
The momentumdependent term of the interaction energy density for a Skyrme effective interaction can be written by using the phasespace distribution function () as
(27) 
The coefficients are related to the Skyrme parameters by
(28) 
Employing the fact that the momentum dependence is quadratic, it is possible to write
(29) 
where several kinds of densities are defined by
(30)  
(31)  
(32) 
Although can be arbitrarily chosen for Eq. (29) to hold, we will choose the local average momentum
(33) 
so that the term in Eq. (29) can be neglected keeping the Galilei invariance.
For the distribution function of Eq. (11) in the case of AMD, we can analytically perform the momentum integration to have
(34) 
For the application to high energy collisions, we now modify the quadratic momentum dependence to the same momentum dependence as Ref. GBDG by modifying to
(35) 
and then use the modified momentumdependent part of the interaction energy density
(36) 
In the present work, we choose the parameter MeV/ for the momentum scale.
Appendix B Correlations to bind clusters
Many of light nuclei (Li, Be etc.) have only one or a few bound states which may be regarded as bound states of internal clusters. The quantummechanical probability of forming such a nucleus is not consistent with the semiclassical phase space with which it can be formed in the standard treatment of AMD. Therefore, for a better description, intercluster correlation is introduced as a stochastic process of binding clusters.
The basic idea is to replace the radial component of the relative momentum between clusters by zero if moderately separated clusters ( fm) are moving away from each other with a small relative kinetic energy [ and MeV]. In addition to these conditions, linking is allowed only if each of the two clusters is one of the three closest clusters of the other when the distance is measured by , so that linking usually occurs in dilute environment. Nonclustered nucleons are treated here in the same way as clusters but two nucleons are not allowed to be linked. Two clusters also should not be linked if they can form an or lighter cluster due to the combination of their spins and isospins. It is possible that more than two clusters are linked by this condition. However, only in the case that the mass number of the linked system is , the binding is performed for the linked system by eliminating the radial velocities of clusters in the centerofmass frame of the linked system.
The energy conservation should be achieved by scaling the relative radial momentum between the centerofmass of the linked system and a third cluster. A reasonable way to choose a third cluster may be to find a cluster which has participated in a collision that formed one of the clusters in the linked system. However, since we do not keep the full history of collisions in our computation, we choose a cluster that has the minimal value of
(37) 
as the third cluster for energy conservation, where and are the distance and the radial component of the kinetic energy for the relative motion between the linked system and the third cluster. The factor with the angle between the relative coordinate () and velocity () is introduced so as to favor the case of .
Appendix C Test particles for the AMD phasespace distribution
Here we describe a method to generate test particles following the onebody phasespace distribution function given by Eq. (11). We generate test particles, where is the number of nucleons (for each spinisospin state) in the system.
Let us first consider a distribution given by a sum of Gaussian distributions
(38) 
with
(39) 
where the ‘physical coordinates’ AMD are used as the centroids . The case of corresponds to the usual wavepacket molecular dynamics (MD) without antisymmetrization. A natural idea in MD is to sample a test particle from each Gaussian distribution . It should be noted that this MD sampling introduces manybody correlations among test particles, which is different from sampling test particles independently following the total distribution .
Our aim here is to extend the MD sampling, with reasonable correlations, for the total onebody distribution . We may decompose into terms as
(40) 
with . The average number of test particles to be generated for each term should be
(41) 
which is evaluated by sampling many points () from the Gaussian distribution . With these average numbers , the actual numbers of test particles , which should be integers, are randomly determined in such a way that and or .
For each term of Eq. (40), test particles should be sampled with the relative weight function , which is a straightforward numerical procedure. However, we may introduce additional correlations among test particles by modifying as
(42) 
when a test particle is generated. Test particles are generated sequentially in a random order, and the modification of by the th test particle influences only the test particles generated after .
In the present work, we have chosen the parameters , and .
References
 (1) M. B. Tsang et al., Phys. Rev. C 86, 015803 (2012).
 (2) C. J. Horowitz, E. F. Brown, Y. Kim, W. G. Lynch, R. Michaels, A. Ono, J. Piekarewicz, M. B. Tsang and H. H. Wolter, J. Phys. G: Nucl. Part. Phys. 41 093001 (2014).
 (3) B. A. Li, L. W. Chen and C. M. Ko, Phys. Rep. 464, 113 (2008).
 (4) B. A. Li, Phys. Rev. Lett. 88, 192701 (2002); B. A. Li, Nucl. Phys. A 708, 365 (2002); B. A. Li, Phys. Rev. C 67, 017601 (2003).
 (5) G. F. Bertsch, Nature 283, 280 (1980).
 (6) A. Bonasera and G. F. Bertsch, Phys. Lett. B 195, 521 (1987).
 (7) G. Ferini, M. Colonna, T. Gaitanos and M. Di Toro, Nucl. Phys. A 762, 147 (2005).
 (8) T. Gaitanos, M. Di Toro, S. Typel, V. Baran, C. Fuchs, V. Greco and H. H. Wolter, Nucl. Phys. A 732, 24 (2004).
 (9) W. Reisdorf et al. (FOPI Collaboration), Nucl. Phys. A 781, 459 (2007).
 (10) Z. Xiao, B. A. Li, L. W. Chen, G.C. Yong, and M. Zhang, Phys. Rev. Lett. 102 (2009) 062502.
 (11) Z. Q. Feng and G. M. Jin, Phys. Lett. B 683 (2010) 140.
 (12) J. Hong and P. Danielewicz, Phys. Rev. C 90 (2014) 024605.
 (13) WenMei Guo, GaoChan Yong and Wei Zuo, Phys. Rev. C 90 (2014) 044605.
 (14) W. Reisdorf et al. (FOPI Collaboration), Nucl. Phys. A 848, 366 (2010).
 (15) Symmetry Energy Project, https://groups.nscl.msu.edu/hira/sepweb/pages/home.html
 (16) A. Ono, H. Horiuchi, T. Maruyama, and A. Ohnishi, Prog. Theor. Phys. 87, 1185 (1992).
 (17) Y. Nara, N. Otuka, A. Ohnishi, K. Niita, and S. Chiba, Phys. Rev. C 61, 024901 (1999).
 (18) A. Ono, J. Phys.: Conf. Ser. 420, 012103 (2013).
 (19) M. Isse, A. Ohnishi, N. Otuka, P. K. Sahu and Y. Nara, Phys. Rev. C 72, 064908 (2005).
 (20) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
 (21) C. Gale, G. Bertsch, and S. Das Gupta, Phys. Rev. C 35 (1987) 1666.
 (22) P. Danielewicz and G.F. Beartsch, Nucl. Phys. A 533, 712 (1991).
 (23) S. Teis, W. Cassing, M. Effenberger, A. Hombach, U. Mosel and G. Wolf, Z. Phys. A 356, 421 (1997).
 (24) J. Weil, H. van Hees and U. Mosel, Eur. Phys. J. A 48, 111 (2012) [Eur. Phys. J. A 48, 150 (2012)].
 (25) S. A. Bass et al., Prog. Part. Nucl. Phys. 41, 255 (1998).
 (26) H. Kruse, B. V. Jacak and H. Stoecker, Phys. Rev. Lett. 54, 289 (1985).
 (27) G. Ferini, T. Gaitanos, M. Colonna, M. Di Toro and H. H. Wolter, Phys. Rev. Lett. 97, 202301 (2006).
 (28) J. Xu, C. M. Ko, and Y. Oh, Phys. Rev. C 81, 024910 (2010).
 (29) J. Xu, L. W. Chen, C. M. Ko, B. A. Li and Y. G. Ma, Phys. Rev. C 87, 067601 (2013).
 (30) B. A. Li, Phys. Rev. C 92, 034603 (2015).
 (31) M. B. Tsang, private communication.