1. Introduction
The nature of dense nuclear matter can be studied in various ways: theoretical approaches, heavy ion experiments, and astrophysical observations. Theoretical approaches include estimation of the equation of state (EOS) of dense nuclear matter and heavy-ion collision simulations with transport models. The information on nuclear EOS can be extracted from heavy-ion collisions (HICs), because they can provide opportunities to explore densities above nuclear saturation density (
) in laboratory experiments. Many accelerator facilities for heavy-ion experiments such as the Rare isotope Accelerator complex for ON-line experiment (RAON) in Korea, the Facility for Rare Isotope Beams (FRIB) in the United States, the GSI Facility for Antiproton and Ion Research (FAIR) in Germany, the Radioactive Isotope Beam Factory (RIBF) at RIKEN in Japan, the Heavy Ion Research Facility in Lanzhou (HIRFL) in China, and the Second Generation System Online Production of Radioactive Ions (SPIRAL2) at GANIL in France [
1,
2,
3,
4,
5,
6,
7] are in operation or under construction. In the astrophysical side, in addition to the estimation of the masses and radii of neutron stars by electromagnetic wave observations, the tidal deformability of a neutron star has been estimated by gravitational wave detections, GW170817 [
8,
9]. These astrophysical multi-messenger observations can give strong constraints on the dense nuclear matter EOS.
Symmetry energy is a significant quantity in nuclear matter EOS, which describes how EOS depends on the isospin asymmetry,
[
10,
11,
12,
13,
14,
15]. Symmetry energy for a symmetric nuclear matter (
) is relatively well constrained [
16,
17,
18]. However, symmetry energy for asymmetric nuclear matter is still quite uncertain and it can be challenging to better understand it. Pion production in HICs can be a significant and sensitive process which can provide valuable constraints on the symmetry energy in an asymmetric nuclear matter [
19]. To produce dense nuclear matter whose densities are above
, the beam energy should be high enough to exceed the pion-production threshold. In this case,
baryons can be produced by inelastic collisions of nucleons. They will then decay to nucleons by emitting pions. Hence, understanding the process of
production and decay is essential in estimating not only the multiplicities of charged pions but also the
ratio. For instance, in Refs [
15,
19], it was claimed that the pion ratio is proportional to
.
The transport model has an advantage in that it can describe dynamical processes in a system created by two colliding nuclei [
20,
21,
22,
23]. In this work, we employ a recently developed BUU-type transport model, DaeJeon Boltzmann–Uehling–Uhlenbeck (DJBUU), to study the pion production in intermediate energy HICs [
24,
25]. The first version of DJBUU with only elastic collisions was successful in describing Au+Au collisions, box calculation using collisions, and mean-field dynamics [
25,
26,
27]. In this work, DJBUU has been extended to include pion productions via inelastic collisions.
In
Section 2, we introduce the formalism of a relativistic mean-field model. In
Section 3, we describe the pion production mechanism in the DJBUU model with numerical schemes including basic transport equations under the relativistic mean-field potentials. Inelastic collision for
baryons and their decay process are also described in this section. In
Section 4, we summarize results of Sn+Sn collisions with different combinations of isotopes (
n-rich and
n-poor) at
MeV. The results are compared with S
RIT experiments. Finally, we summarize our results and present discussions in
Section 5.
2. Relativistic Mean-Field Model
The relativistic Lagrangian density without an isovector-scalar field is given by [
28]
where
,
, and
. The nucleon, isoscalar-scalar, isoscalar-vector, and isovector-vector fields are represented by
, and
, respectively. Masses of these particles are
, and
, and the coupling constants of mesons to nucleons are given by
, and
. The nonlinear potential of
-meson self-interaction is specified by the strengths
a and
b. The parameter set we use satisfies the properties of the nuclear matter such as the saturation density of
fm
, the binding energy
MeV, the incompressibility
MeV, the nucleon effective mass
, and the symmetry energy
MeV, which is summarized in
Table 1 [
28]. For particles with an electric charge, we also include the electromagnetic field (
). In our mean-field treatment, we assume that the spatial components of the vector meson fields,
and
, are zero and we only take the time component. In principle, spatial components of meson fields have to be included. However, in this study the spatial components of all vector mesons are neglected as we only deal with uniform nuclear matter.
The baryon scalar density, baryon current, and the third component of isospin current are
with
and
. The factor
represents the degeneracy of the baryon species
a. We assume that not only nucleons but also
baryons in the calculations have the same amount of shift from their vacuum masses due to the
-meson interactions. For the full treatment of
baryons in medium, Equation (1) has to be extended by including interaction terms of
baryons. However, the in-medium mass shift of the
baryon is still an open question and the quark components of
baryons are the same as nucleons. As such, we assume that the
baryons couples to the
field in the same way the nucleons do. In our mean-field approximation, the meson field equations become
where
3. Pion Production in Daejeon BUU Transport Model
The transport model describes the time evolution of the nucleon phase-space distribution function
under the relativistic mean-field potentials as follows:
where the left hand side represents the propagation of distribution function
in the relativistic mean-field potentials and the right hand side describes all possible collision processes, elastic and inelastic scatterings, and decay of resonances. We use the test particle ansatz to solve the integro-differential non-linear BUU equation numerically [
29]. The phase-space density and cross sections are scaled with the number of test particles
DJBUU has a distinct feature in calculating the phase-space density distribution. In most transport models, Gaussian shape is used as the particle profile [
26,
27,
30], but in DJBUU, a polynomial shape is used. The polynomial shape is designed to have finite endpoints in the profile and it is integrable. The form of the polynomial function used in DJBUU is
with
and
. We verified that the polynomial shape with these parameters gives nearly the same distribution as a Gaussian shape in the density region where the particle production is dominant.
For nucleons and
resonances, we use an isotropic constant cross section
mb for elastic scatterings, which is in good agreement with the decay width of the isoscalar giant quadrupole resonance [
31]. The energy-dependent isospin-averaged isotropic cross section for
resonance production by nucleon–nucleon collisions is
where the cm energy is above the threshold
GeV [
23]. Because the beam energy we study is just above the threshold for the production of
resonances, we only consider the inelastic reaction
. Details of the isospin-dependent cross sections are summarized in
Table 2.
The mass distribution of
resonance is assumed to be in the Breit–Wigner form with density-dependent effective masses [
32,
33]
where
to satisfy
and
is the effective mass that is the peak mass (
GeV) shifted by the
-meson mean field. The mass-dependent width in Equation (12) is taken from phenomenological approaches [
34,
35];
where the momentum
q of the pion in the
rest frame is given as
The effective mass of
resonance is sampled according to
where
All asterisk marks represent the effective masses and canonical energy-momentum by meanfield potentials. The lower and upper limits of
mass are
. The cross section of the inverse reaction, which satisfies the detailed balance condition, is given as
where
, and
g is the spin-isospin factor.
In the pion production from
resonance decays, pions are treated as free particles with the bare mass (
MeV) and their momenta are sampled isotropically in the rest frame of the progenitor particles (
baryons). For the process
, the isospin-averaged cross section from the detailed balance condition is
where
is the nucleon or pion momentum in their cm frame and
is the cm energy of the nucleon and pion. Details of isospin-dependent cross sections, spin-isospin factors, and decay widths are summarized in
Table 2.
In the relatively low-beam-energy region, the mean fields are more likely attractive and result in the increase in collision rates compared to the case without mean fields (Cascade mode) [
25,
26]. In this sense, the number of
collisions is simply increased due to the presence of mean fields and the probability of
resonance production is increased accordingly. In Ref. [
36], it was already estimated that the total
yields with meanfield potentials are almost three times larger than those in the Cascade calculations (without mean-field potentials). One way to reduce both
and
numbers (to be consistent with experimental results) is to introduce density-dependent cross sections that can include the nuclear medium effect indirectly [
36,
37]. In addition to the density-dependent effect, we also introduce the isospin-dependent effect for
cross sections because
ratios in HICs are strongly correlated with isospin asymmetry of the reaction system [
19]. Our cross sections for
resonances are
where
is
cross sections in the free space as given in Equation (11),
is the baryon number density,
C is a fitting parameter for the density dependence,
Z and
N are the total number of protons and neutrons of a collision system, and
are the fitting parameters for the collision system dependence. In principle, given an underlying Lagrangian, one should be able to calculate cross sections that depend on the medium density, isospin asymmetry, and the mean fields. Unfortunately, such calculation is currently not very feasible. What we would like to do with the parametrization in Equation (19) is to see if this parametrization is flexible enough to capture the essential physics. It is possible to use the time-dependent local density ratio (
) in Equation (19) instead of
. However, we see in our simulations that the
and the average of (
are consistent within a few percent in the regions where
resonances are produced. As such, we will use
for simplicity.
In the previous studies [
36,
37] with no
,
C was taken to lie between 1.65 and 1.9 in order to compare pion-related quantities with the FOPI experiment data. The isospindependent parameters
,
, and
correspond to the
,
, and
channels, respectively. For the asymmetric collisions with a high
ratio, these factors will give strong constraints to the pion productions. The cross sections for the inverse reactions also have density and isospin dependences as in Equation (17).
Once
baryons are produced from
scatterings, they should eventually decay into one nucleon and one pion with corresponding electric charges. With a given decay width Equation (13), we determine whether a
baryon decays or propagates at a given time step according to the decay probability within that time step. To avoid the spurious decay of resonance, two different time intervals are considered in our calculation. If the resonance is scheduled to collide, we use the modified time interval
. Otherwise, we use the original time interval,
. The probability of decay within a time step is then given by
where
is the mass-dependent decay width of
resonance and
depends on whether the resonance is scheduled to undergo a collision or not.
The dominant collisions from nucleon–nucleon elastic scatterings are governed by the Pauli blocking effect. This can be interpreted as an in-medium effect of the elastic scattering. We neglect the Pauli blocking effect for baryons, but we fully consider the blocking effect for nucleons. For example, the typical Pauli blocking factor for a nucleon–nucleon scattering is , while the probability for not being Pauli blocked for or is .
4. Numerical Results
In our transport model, baryons propagate under the mean-field potentials according to the classical equations of motion
Under the mean-field potentials, each baryon species is affected by different vector potentials, because nucleons and
baryons have different isospin values and electric charges. In the absence of the mean-field potentials, baryons propagate in the phase space in only the coordinate space unless a collision occurs. However, under the mean-field potentials, the time evolution of the phase-space functions of baryons are governed by Equation (21). In this case, both the spatial and momentum coordinates change. For pions, because they are treated as free particles, and the only way to change their momentum is a collision regardless of the existence of mean fields. We take the number of test particles
, the time interval
fm/
c in all calculations to reduce spurious collisions or decays, and 40 independent runs for a better statistics. The density distributions from relativistic Thomas–Fermi calculations are used for initializing the projectile and the target nuclei.
Two combinations of Sn+Sn collision systems are considered in this work: neutron rich (
Sn+
Sn) and neutron poor (
Sn+
Sn). We take the impact parameter
fm and the incident beam energy
MeV which corresponds to the beam energy and impact parameter of the S
RIT experiment [
38]. For the cross section of
baryons, we examine four sets of (
) combinations as in
Table 3.
Figure 1 shows the number of
resonances, pions, and the central density divided by the saturation density as a function of time. Regardless of the parameter sets we are using, the maximum densities of Sn+Sn collisions reach about 1.6
near
fm/
c. During the compressing stage (
fm/
c) where the density increases,
baryons are actively produced from
scatterings. After the peak of
production around
fm/
c,
production decreases rapidly and the
’s mostly decay to nucleons and pions. This is clearly seen at
fm/
c. At the later stage of the time evolution, most pions are produced from
decay process and the number of pions almost saturates. Although the other results (
, and
) are not shown in
Figure 1, the trends of all quantities with respect to time are not so sensitive to the choice of parameter sets.
Numerical results for pion yields for the four studied cases are shown in
Figure 2a. Empty boxes in this figure are experimental results from the S
RIT experiment. In
Figure 2a, the variations in the
prediction (lower lines) are smaller than those of the
prediction (upper lines). Note that the isospin-dependent parameters
critically affect the
production of neutron-rich collision systems as we expected. However, for neutron-poor systems (
), pion yields rarely depend on the different cases. In
Figure 2b, the pion single ratios are summarized. For the neutron-rich system with
(
Sn+
Sn), differences in the single ratios are much larger compared to the neutron-poor system with
(
Sn+
Sn). All pion ratios with
are close to the experiment data regardless of the parameter sets. For
,
is most consistent with the experiment. The errors in
Figure 2b are estimated by considering the error propagation of each pion yields in
Figure 2a. In the experimental aspect of pion measurement, the pion double ratio (DR),
can reduce sensitivities to specific details on both collisions and mean-field potentials, as well as the treatment of the Coulomb interaction. As we use Sn isotopes for both the projectile and the target, the uncertainties due to the Coulomb effects can be significantly reduced in the double ratio. Our estimation of the double ratios for the four different cases are shown in
Figure 2c. The gray horizontal bar represents the experimental value. In Ref. [
38], the authors considered two different cases for symmetry energy (soft and stiff) and the height of each colored box represents the difference between the soft and stiff cases. Comparing with these results, our results are in better agreement with the experiment data. Our predictions of the pion double ratio are between 2.2 and 2.6. All results in the figures are directly comparable to the predictions from seven widely used transport models.
Table 4 summarizes the pion multiplicities, single ratios, and double ratios for the four different density- and isospin-dependent parameter sets and the S
RIT experiment data. All produced
and all
are slightly larger than the experiment results for
Sn+
Sn collisions. On the other hand, all produced
are slightly smaller than the experimental results for
Sn+
Sn collisions. Nevertheless, by introducing the isospin-dependent term in the
cross sections, we have improved agreement with the experiment data compared to previous studies.
In the transport simulation, one can keep track of the processes in HICs as a function of time. In
Figure 3, we show the collision rates of the total (not divide by charge states)
baryons (left panels) and pions (right panels). The collision rate of the produced particles (by collisions or decays) are summarized in the top panel, and the rate of the absorbed particles (by the inelastic collisions) are summarized in the middle panel. The net
,
) −
) is summarized in the bottom panel. The negative values of net
indicate that the inverse reaction (
) is predominant in the time interval
fm/
c. From the results of the
production in the left panels, we can estimate that more than 70% of the produced
baryons are returned to nucleons by participating in the inverse reactions (
). In the early stage (up to about
fm/
c) where projectile and target are compressed, production of
baryons wins over decays. The situation is reversed afterwards. For pion cases, this trend is somewhat different from that of the
baryon. Both pion production and the absorption peak around
fm/
c, while
peaks near
fm/
c. It is obvious that
production should be dominant over other processes such as inverse or decay processes, and other processes should occur after
production as a chain reaction. There are no negative values in the net pion production (right bottom panel) at any time. The results for the neutron-poor system with
(
Sn+
Sn) show a similar tendency, even though the absolute values are different.
In
Figure 4, we analyze the single-pion spectral ratios with the four different parameter sets for both (a) neutron-rich (
) and (b) neutron-poor (
) systems. In general, single ratios are high at low
and decrease as
increases (even though there exist significant deviations at high
above 300 MeV/c). This trend can be interpreted as an effect of Coulomb interaction [
39]. The Coulomb interaction changes momenta of charged pions by acting as an attractive force for
and repulsive one for
due to the presence of many protons around them. the results are rather flat with respect to
. For both collision systems, most of the theoretical predictions are higher than the experiment data at
MeV/
c. In the high
region above 300 MeV/c, the discrepancies between the results with different parameter sets become large for both collision systems. Much larger single ratios above 300 MeV/c in
for the neutron-poor system can be caused either by the large number of
or by the small number of
or by both. One proton is produced in the
production channel (
), while one proton is consumed in the
production channel (
). Because Coulomb effects are very important for pion production, the change in the density of protons can cause changes in the momentum distribution of pions, especially at high
.
As in
Figure 4c, the uncertainties in the double ratios are much larger than those in the single ratios, giving less predictive power on the effect of medium and isospin dependencies. More accurate treatment of Coulomb interactions and proper pion potentials is required to resolve the theoretical overpredictions and discrepancies [
40]. Subtracted and isoscaling ratios are alternative quantities to understand the mechanism of pion production [
41].