Pion Productions with Isospin-Dependent In-Medium Cross Sections

: Total pion yields and π − / π + ratios in two Sn+Sn collision systems, 132 Sn+ 124 Sn (neutron rich) and 108 Sn+ 112 Sn (neutron poor) at E = 270 A MeV, are signiﬁcant observables in SAMURAI Pion-Reconstruction and Ion-Tracker (S π RIT) experiments. Based on a recently developed transport model, DaeJeon Boltzmann–Uehling–Uhlenbeck (DJBUU), we investigate the isospin-dependent in-medium effects by including density- and isospin-dependent cross sections for ∆ baryon productions. We compare our results with the S π RIT data by considering these effects. We ﬁnd that the yields as well as the ratios strongly depend on the in-medium effect, especially isospin-dependent in-medium cross sections for ∆ resonances.


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 (ρ > ρ 0 ) 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, δ = (ρ n − ρ p )/ρ [10][11][12][13][14][15]. Symmetry energy for a symmetric nuclear matter (ρ n = ρ p ) 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 ρ 0 , 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 (ρ n /ρ p ) 2 .
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 E = 270A MeV. The results are compared with SπRIT experiments. Finally, we summarize our results and present discussions in Section 5.

Relativistic Mean-Field Model
The relativistic Lagrangian density without an isovector-scalar field is given by [28] L where U(σ) = 1 3 aσ 3 + 1 4 bσ 4 , Ω µν = ∂ µ ω ν − ∂ ν ω µ , and R µν = ∂ µ ρ ν − ∂ ν ρ µ . The nucleon, isoscalar-scalar, isoscalar-vector, and isovector-vector fields are represented by ψ, σ, ω µ , and ρ µ , respectively. Masses of these particles are m N , m σ , m ω , and m ρ , and the coupling constants of mesons to nucleons are given by g σ , g ω , and g ρ . 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 ρ 0 = 0.16 fm −3 , the binding energy E/A = −16 MeV, the incompressibility K = 240 MeV, the nucleon effective mass m * = 0.75 m N , and the symmetry energy S(ρ 0 ) = 30.5 MeV, which is summarized in Table 1 [28]. For particles with an electric charge, we also include the electromagnetic field (A µ ). In our mean-field treatment, we assume that the spatial components of the vector meson fields, ω µ and ρ µ 3 , 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 m * B = m B + g σ σ and p 0 a = m * 2 B + p 2 . The factor g a 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

Pion Production in Daejeon BUU Transport Model
The transport model describes the time evolution of the nucleon phase-space distribution function f (x, p; t) under the relativistic mean-field potentials as follows: where the left hand side represents the propagation of distribution function f (x, p; t) 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 particleŝ 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 m = 2 and n = 3. 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 σ elastic BB→BB = 40 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 √ s ≥ 2.015 GeV [23]. Because the beam energy we study is just above the threshold for the production of ∆ resonances, we only consider the inelastic reaction NN → N∆. Details of the isospin-dependent cross sections are summarized in Table 2. Table 2. Isospin-dependent channels for ∆ resonances and pion productions with the corresponding cross sections (σ), spin-isospin factors (g), and the decay widths (Γ).
The mass distribution of ∆ resonance is assumed to be in the Breit-Wigner form with density-dependent effective masses [32,33] whereÃ = 0.95 to satisfy dm 2π A = 1 and m * 0 is the effective mass that is the peak mass (m 0 = 1.232 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 m * The cross section of the inverse reaction, which satisfies the detailed balance condition, is given as where p 2 f = s * /4 − m * 2 N , and g is the spin-isospin factor. In the pion production from ∆ resonance decays, pions are treated as free particles with the bare mass (m π = 139 MeV) and their momenta are sampled isotropically in the rest frame of the progenitor particles (∆ baryons). For the process N + π → ∆, the isospin-averaged cross section from the detailed balance condition is where p cm is the nucleon or pion momentum in their cm frame and √ s * 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 NN 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 σ NN→N∆ (0) is ∆ cross sections in the free space as given in Equation (11), ρ B 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 x ±, 0 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 (ρ n /ρ p ) in Equation (19) instead of (N/Z) sys . However, we see in our simulations that the (N/Z) sys and the average of (ρ n /ρ p ) local are consistent within a few percent in the regions where ∆ resonances are produced. As such, we will use (N/Z) sys for simplicity.
In the previous studies [36,37] with no x ±, 0 , 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 x + , x 0 , and x − correspond to the pp, np, and nn channels, respectively. For the asymmetric collisions with a high n/p 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 NN 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 dt = t coll. − t now . Otherwise, we use the original time interval, dt = t next − t now . The probability of decay within a time step is then given by where Γ is the mass-dependent decay width of ∆ resonance and dt 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

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 N test = 100, the time interval dt = 0.2 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 ( 132 Sn+ 124 Sn) and neutron poor ( 108 Sn+ 112 Sn). We take the impact parameter b = 3 fm and the incident beam energy E = 270A 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 (C, x + , x 0 , x − ) 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ρ 0 near t = 15 fm/c. During the compressing stage (t = 10 ∼ 20 fm/c) where the density increases, ∆ baryons are actively produced from NN scatterings. After the peak of ∆ production around t = 20 fm/c, ∆ production decreases rapidly and the ∆'s mostly decay to nucleons and pions. This is clearly seen at t = 20 ∼ 40 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 (case1, case2, and case3) 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 x ±, 0 critically affect the π − production of neutron-rich collision systems as we expected. However, for neutron-poor systems (N/Z = 1.2), pion yields rarely depend on the different cases. In Figure 2b, the pion single ratios are summarized. For the neutron-rich system with N/Z = 1.56 ( 132 Sn+ 124 Sn), differences in the single ratios are much larger compared to the neutron-poor system with N/Z = 1.2 ( 108 Sn+ 112 Sn). All pion ratios with N/Z = 1.2 are close to the experiment data regardless of the parameter sets. For N/Z = 1.56, case4 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 132 Sn+ 124 Sn collisions. On the other hand, all produced π ± are slightly smaller than the experimental results for 108 Sn+ 112 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. . We consider four difference cases as given in Table 3. (b) Charged pion single ratios as a function of N/Z for four cases. All the empty boxes in (a,b) are experimental data with uncertainties. (c) Double ratios for 132 Sn+ 124 Sn and 108 Sn+ 112 Sn for four different cases. Experimental data are marked as a gray box. Note that the statistical errors are quite large because the number of produced pions per run is very small. Table 4. Pion multiplicities Y(π ± ), single multiplicity ratio SR(π − /π + ), and double multiplicity ratio DR(π − /π + ) of Sn+Sn collisions for the four cases. The experimental data of Sn+Sn collisions are from the SπRIT Experiment [38].
132 Sn+ 124 Sn @ 270A MeV 108 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 ∆, dN/dt(NN → N∆) − dN/dt(N∆ → NN) is summarized in the bottom panel. The negative values of net ∆ indicate that the inverse reaction (N∆ → NN) is predominant in the time interval t = 20 ∼ 40 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 (N∆ → NN). In the early stage (up to about t = 20 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 t = 20 fm/c, while NN → N∆ peaks near t = 17 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 N/Z = 1.2 ( 108 Sn+ 112 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 (N/Z = 1.56) and (b) neutron-poor (N/Z = 1.2) systems. In general, single ratios are high at low p T and decrease as p T increases (even though there exist significant deviations at high p T 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 p T . For both collision systems, most of the theoretical predictions are higher than the experiment data at p T > 100 MeV/c. In the high p T 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 case2 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 (NN → N∆ −,0 → N pπ − ), while one proton is consumed in the π + production channel (NN → N∆ +,++ → Nnπ + ). 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 p T .
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].   (c) Transverse momentum spectra of the double ratio of two Sn collision systems, Equation (22). Circle markers with uncertainties are experiment data taken from Ref. [40] and the curves are the ratios corresponding to four different parameter sets.

Summary and Discussion
In this work, using the recently developed transport model DJBUU, we have studied the isospin-dependent in-medium effects by introducing density-and isospin-dependent cross sections for ∆ baryon production. Because pions are produced as the results of ∆ resonance decay, the number of pions is closely related to the number of ∆ baryons. We compare our theoretical pion yields and single and double ratios with the measurements of charged pions in Sn+Sn systems using the SπRIT Time Projection Chamber, revealing values of <4% in multiplicity and ≈ 2% in ratio [38]. We found that the effect of dense nuclear matter as well as isospin asymmetries are very important to explain the pion productions. We found that the density-dependence parameter C has to be about 2.5, which is slightly larger than previous results, and the isospin dependency parameter x ±, 0 is crucial in ∆ and pion productions.