Introduction

The isospin asymmetric nuclear equation of state is crucial for understanding isospin-asymmetric objects such as the structure of neutron-rich nuclei, the mechanism of neutron-rich heavy ion collisions (HICs), and the properties of neutron stars, including neutron star mergers and core-collapse supernovae [1-10]. The symmetric part of the isospin-asymmetric equation of state has been well constrained using flow and kaon condensation [11]. However, the symmetry energy away from the normal density still has a large uncertainty, making the constraint of symmetry energy becomes one of the important goals in nuclear physics [12, 13].

For probing the symmetry energy at suprasaturation density using HICs, the isospin-sensitive observables, such as the ratio of elliptic flow of neutrons to charged particles, hydrogen isotopes, or protons ($v_2^{\text{n}}/v_2^{\text{ch}}$, $v_2^{\text{n}}/v_2^{\text{H}}$, or $v_2^{\text{n}}/v_2^{\text{p}}$) [14-18] and the multiplicity ratio of charged pions (i.e., M(π^-)/M(π⁺) or denoted as π^-/π⁺) [19-30], were mainly used. By comparing the calculations with transverse momentum-dependent or integrated elliptic flow data of nucleons, hydrogen isotopes, and charged particles from the FOPI/LAND and ASY-EOS experimental collaborations, a moderately soft to linear symmetry energy was obtained with the UrQMD [31, 14, 16] and Tübingen quantum molecular dynamics (TüQMD) model [15, 18]. The lower limit of L obtained from the flow ratio data is L>60 MeV [32], which overlaps with the upper limits of the constraints from the nuclear structure and isospin diffusion, i.e., $L\approx 60\pm 20$ MeV [33-35]. It shows strong model dependence for the symmetry energy constraints from π^-/π⁺ [21-23, 25-28, 36] and the extracted value of L ranges from 5 MeV to 144 MeV. This may be due to the different treatments of nucleonic potential, Δ potential, threshold effects, pion potential, Pauli blocking, in-medium cross-sections, and by using the different numerical techniques for solving transport equations.

To understand the model dependence of the symmetry energy constraints of HICs and improve it in the future, two aspects should be considered. One is to find and fix the deficiencies of transport models, which can be achieved through a transport model evaluation project. The other is to test the model by simultaneously describing the multi-observable data and then providing the constraints of symmetry energy at their probed densities.

The transport model evaluation project has made important progress in benchmarking the treatment of particle-particle collision [37, 38] and nucleonic mean field potential [39] in both the Boltzmann-Uehling-Uhlenbeck (BUU) type and quantum molecular dynamics (QMD) type models. In the collision part, a time-step-free method is suggested [37, 38] for simulating the collisions or decay of resonance particles because it automatically determines whether the resonance will collide or decay according to the sequence of the collision time and decay time. In the UrQMD model, the time-step-free method is adopted in the collision part [37, 38], and the nucleonic potential is also involved in extending its application to low-intermediate energy HICs [16, 28].

Despite the successful applications of the UrQMD model in studying heavy-ion collisions (HICs) across a range of energies from low-intermediate to high energy [40, 41, 16, 28, 32], previous calculations have revealed that the data of pion multiplicity and the nucleonic flow observables for Au+Au at 0.4A GeV were not simultaneously described. This ‘inconsistency’ may be attributed to the momentum-dependent interaction (MDI) form and near-threshold ${\sigma }_{NN\to N\text{Δ}}$ cross-sections used in the UrQMD model [16, 28].

The MDI form used in the previous analyses on the elliptic flow of neutrons, protons, hydrogen isotopes, charged particles [16] or the pion multiplicity [28] is $t_4{{\displaystyle \ln }}^2\left(1+t_5{\left({\text{p}}_1-{\text{p}}_2\right)}^2\right)\delta \left({\text{r}}_1-{\text{r}}_2\right)$, in which the MDI parameters were extracted by fitting the Arnold’s optical potential data [42]. By using this MDI form in the UrQMD, the transverse momentum-dependent elliptic flow for neutrons and hydrogen isotopes is underestimated by 40% in the high p_T/A region [31]. In the 1990s, the real part of the global Dirac optical potential (Schrodinger equivalent potential) was published by Hama et al. [43], in which the angular distribution and polarization quantities in proton-nucleus elastic scattering in the 10 MeV to 1 GeV range were analyzed. Based on Hama’s data, a Lorentzian-type momentum-dependent interaction [44] was generated and used in the IQMD model in Ref. [44], and in many versions of transport model, such as the Boltzmann-Nordheim-Vlasov (BNV) [45], IBUU [46-48], the jet AA microscopic transportation model + relativistic version of the QMD model (JAM+RQMD) [49], previous version of the UrQMD [50], the Giessen-BUU model (GiBUU) [51], the antisymmetrized molecular dynamics approach (AMD) +JAM [52], RQMD [53], TüQMD [54] model for studying intermediate-high energy HICs. The MDI [44] generated from Hama’s data provides a stronger momentum-dependent potential than that from Arnold’s data in the high-momentum region. Thus, checking whether Hama’s MDI form can refine the nucleonic flow description in the UrQMD model is important.

For the cross sections of the $NN\to N\text{Δ}$ channel used in the UrQMD model, they are obtained by fitting CERN8401 data [55]. This fitting formula underestimates the data for ${\sigma }_{NN\to N\text{Δ}}$ near the threshold energy, which will be shown in Fig. 2 in Sect. 2, and thus leads to an underestimation of the pion productions in the UrQMD calculation. The other transport models, such as TüQMD, pBUU, and RVUU, use the ${\sigma }_{NN\to N\text{Δ}}$ cross-section obtained from the one-boson exchange model by fitting CERN8301 data. However, one should note that the data from CERN8301 and CERN8401 were different, particularly near the threshold energy. Thus, investigations of the different formulas of ${\sigma }_{NN\to N\text{Δ}}$ near the threshold energy are necessary for describing the pion observables.

Another method for reducing model uncertainties is to simultaneously describe the multi-observables data (or doing so-called combination analysis), including the isospin-independent and isospin-dependent nucleonic collective flow and pion observables. For the combination analysis on the isospin sensitive nucleonic flow and pion observables, there were few works to simultaneously investigate them, except for the TüQMD model [27] and IBUU model [17]. In the TüQMD model, the medium correction on the cross-sections, energy conservation, and momentum-dependent symmetry potential have been considered, and four observables, such as $M_{{\pi }^{+}}$, $M_{{\pi }^{-}}$, π^-/π⁺, and integral $v_2^{\text{n}}/v_2^{\text{p}}$, were analyzed. In the IBUU model, the nucleon-nucleon short-range correlations and isospin-dependent in-medium inelastic baryon-baryon scattering cross-sections were considered, and six integrated flow and pion observables were analyzed.

In the last decades, ASY-EOS, FOPI-LAND, and FOPI have published 17 datasets on nucleonic collective flows and pion observables, as listed in Table 3, which provides a significant opportunity to benchmark the model and understand the contributions of the different physical phenomena. In this study, we attempted to use 17 observables to limit the physical uncertainties and improve the ability of the UrQMD model.

The paper is organized as follows: In Sect. 2, we briefly introduce the MDI, symmetry energy, and ${\sigma }_{NN\to N\text{Δ}}$ that will be refined in the UrQMD model. In Sect. 3, the impacts of MDI, symmetry energy, and the refined ${\sigma }_{NN\to N\text{Δ}}$ on the 17 observables, such as nucleonic flow and pion observables, are presented and discussed. In Sec. 4, the symmetry energy constraints at the flow and pion characteristic densities are obtained and the model dependence of them are discussed. Sect. 5 concludes this study.

UrQMD model and its refinements

The UrQMD model version we used is the same as that used in Ref. [28], in which the cross-sections of the $N\text{Δ}\to NN$ channel are replaced with a more delicate form by considering the Δ-mass dependence of the M-matrix in the calculation of the $N\text{Δ}\to NN$ cross-section [56]. This differs from the version used to describe only the flow data and constraints in Refs. [16, 32]. To distinguish them, we named the current version as UrQMD-CIAE and previous version as UrQMD-HZU. The main differences are the momentum dependence potential and $NN\leftrightarrow N\text{Δ}$ cross-section.

We focused on the impacts of different forms of the MDI, symmetry energy, and ${\sigma }_{NN\to N\text{Δ}}$ and have briefly introduced them in the subsequent sections. The nucleonic potential energy U was calculated from the potential energy density u, i.e., $U={\displaystyle \int }u{\text{d}}^{3}r$. The u reads as $\begin{array}{c}u=\frac{\alpha }{2}\frac{{\rho }^2}{{\rho }_0}+\frac{\beta }{\eta +1}\frac{{\rho }^{\eta +1}}{{\rho }_0^{\eta }}\\ +\frac{g_{\text{sur}}}{2{\rho }_0}{\left(\nabla \rho \right)}^2+\frac{g_{\text{sur,iso}}}{{\rho }_0}{\left[\nabla \left({\rho }_{\text{n}}-{\rho }_{\text{p}}\right)\right]}^2\\ +u_{\text{md}}+u_{\text{sym}}\mathrm{.}\end{array}$ (1) The parameters α, β, and η are related to the two and nonlinear density-dependent interaction term. The third and fourth terms are the isospin-independent and isospin-dependent surface term, respectively. The u_md is from the MDI term, and two forms were adopted in this work. u_sym denotes the symmetry potential energy term.

The energy density associated with the MDI, i.e., u_md, is calculated according to the following relationship: $u_{\text{md}}={\displaystyle \sum _{ij}{\displaystyle \int }{d}^3}{p}_1{d}^3{p}_2{f}_i\left(\overrightarrow{r}\mathrm{,}{\overrightarrow{p}}_1\right)f_j\left(\overrightarrow{r}\mathrm{,}{\overrightarrow{p}}_2\right)v_{\text{md}}\left(\text{Δ}{p}_{12}\right)\mathrm{.}$ (2) $f_i\left(\overrightarrow{r}\mathrm{,}{\overrightarrow{p}}_1\right)$ denotes the phase space density of nucleon i. The MDI form, that is, $v_{\text{md}}\left(\text{Δ}{p}_{12}\right)$, is assumed to be $v_{\text{md}}\left(\text{Δ}{p}_{12}\right)=t_4{{\displaystyle \ln }}^2\left(1+t_5\text{Δ}{p}_{12}^2\right)+c\mathrm{,}$ (3) where $\text{Δ}{p}_{12}=\left|{\text{p}}_1-{\text{p}}_2\right|$, and the parameters t₄, t₅, and c are obtained by fitting the data of the real part of the optical potential. In detail, we fit the data of the real part of the nucleon-nucleus optical potential $V_{\text{md}}\left(p_1\right)$ according to the following ansatz: $V_{\text{md}}\left(p_1\right)={\displaystyle {\int }_{p_2<p_F}}{v}_{\text{md}}\left(p_1-p_2\right){\text{d}}^3{p}_2/{\displaystyle {\int }_{p_2<p_{\text{F}}}}{\text{d}}^3{p}_2\mathrm{.}$ (4) This method is the same as that used in Ref. [44].

Two sets of the real part of optical potential data were used in this work. One is from Arnold et al. [42], which was used in the previous version of the UrQMD model[16, 28]. Another is from Hama et al. [43], which is widely used in many transport models, such as BNV[45], IBUU[46-48], JAM+RQMD[49], the previous version of the UrQMD[50], GiBUU[51], AMD +JAM[52], RQMD[53], and the TüQMD[54] model. The momentum dependence of $v_{\text{md}}^{\text{Hama}}\left(\text{Δ}{p}_{12}\right)$ is stronger than that of $v_{\text{md}}^{\text{Arnold}}\left(\text{Δ}{p}_{12}\right)$, and the value of $v_{\text{md}}^{\text{Hama}}\left(\text{Δ}{p}_{12}\right)$ is larger than that of $v_{\text{md}}^{\text{Arnold}}\left(\text{Δ}{p}_{12}\right)$ at high momentum region. The corresponding single-particle potentials are presented in Fig. 1 (a). Because the MDI can influence the EOS, parameters α, β, and η should be readjusted to keep the desired shape. The parameters α, β, and η were readjusted to maintain the incompressibility of the symmetric nuclear matter K₀=231 MeV for the two different MDIs, and the values of the parameters and the corresponding effective mass m^*/m are listed in Table 1.

Fig. 1Full-screen View

(Color online) (a) The parametrization of the bare interaction v_md (lines) as compared to the data points of the real part of optical potential (symbols) from Arnold et al. [42] (green) and Hama et al. [43] (red). (b) Density dependence of the symmetry energy with different S₀ and L values

For the potential energy density of the symmetry energy part, i.e., u_sym, only the local interaction has contribution since the nonlocal term is isospin-independent momentum dependent interaction. We take two forms of $u_{\text{sym}}^{\text{pot}}$ in our calculations: the Skyrme-type polynomial form ((a) in Eq. (5)) and the density power law form ((b) in Eq. (5)), which read as $\begin{array}{c}{u}_{\text{sym}}=S_{\text{sym}}^{\text{pot}}\left(\rho \right)\rho {\delta }^2\\ =\{\begin{array}{ll}\left(A\left(\frac{\rho }{{\rho }_0}\right)+B{\left(\frac{\rho }{{\rho }_0}\right)}^{{\gamma }_{\text{s}}}+C{\left(\frac{\rho }{{\rho }_0}\right)}^{5/3}\right)\rho {\delta }^2\mathrm{,}\hfill & (a)\hfill \\ \frac{C_{\text{s}}}{2}{\left(\frac{\rho }{{\rho }_0}\right)}^{{\gamma }_i}\rho {\delta }^2\mathrm{.}\hfill & (b)\hfill \end{array}\end{array}$ (5) Correspondingly, the density dependence of the symmetry energy is $S(\rho )=\frac{{\hslash }^2}{6m}{\left(\frac{3{\pi }^2\rho }{2}\right)}^{2/3}+S_{\text{sym}}^{\text{pot}}\left(\rho \right)\mathrm{.}$ (6) In Eq.(5), A, B, C, C_s, and γi are the parameters of the symmetry potential directly used in the UrQMD model. They are determined by the symmetry energy values at the saturation density S₀, the slope of the symmetry energy L, and the parameters in Table 1 according to the relationship described in Refs. [35, 57], where S₀=S(ρ₀) and $L=3{\rho }_0\partial S(\rho )/\partial \rho {|}_{{\rho }_0}$. The ranges of S₀ and L are listed in Table 2. In this work, we varied S₀ and L to investigate the influence of the symmetry energy on HIC observables.

For the L<35 MeV case, we used the Skyrme-type polynomial form of $S_{\text{sym}}^{\text{pot}}\left(\rho \right)$ because the simple power law form of the symmetry energy cannot provide reasonable values at a subnormal density. Furthermore, the L < 5 MeV sets are not adopted because the corresponding symmetry energy becomes negative at densities above 2.7ρ₀ and the EOS is not favored by the neutron star properties. Thus, the lower limit of L in our calculations was 5 MeV. For L>35 MeV, we used a simple power law form of $S_{\text{sym}}^{\text{pot}}\left(\rho \right)$. As an example, we present the density dependence of the symmetry energy in Fig. 1 (b) for L = 20, 144 MeV at S₀=30 and 34 MeV.

The symmetry potential of Δ resonance was calculated from the symmetry potential of the nucleon as same as that in Refs. [20, 22, 26-28, 58, 59]. The effects of different strengths of Δ potential on pion production were also investigated in Refs. [28, 60], and the total and differential π^-/π⁺ ratios in heavy ion collisions above the threshold energy were weakly influenced by the completely unknown symmetry (isovector) potential of the Δ(1232) resonance, owing to the very short lifetimes of the Δ resonances.

In the collision term, medium-modified nucleon-nucleon elastic cross-sections were used, as same as that in our previous works [32]. For the $NN\to N\text{Δ}$ cross-sections used in UrQMD [40], a formula used to fit the CERN8401 data was adopted and denoted as ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$. After zooming out the figure of the $NN\to N\text{Δ}$ cross-sections near the threshold in Ref. [40], we found that ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ underestimated the data [55] by approximately 3 mb at E=0.4A GeV. This discrepancy is illustrated in Figs. 2 (a), where the blue line is the default fitting formula in Ref. [40] and the solid symbols represent the CERN8401 data obtained from Ref. [55]. One can expect that the default formula ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ would underestimate the pion multiplicities. Thus, we used an accurate form of ${\sigma }_{NN\to N\text{Δ}}$ near the threshold energy to describe pion production at 0.4A GeV. A Hubbert function form was used to refit the $NN\to N\text{Δ}$ cross-sections data at $\sqrt{s}<2.21$ GeV. That is, $\begin{array}{c}{\sigma }_{NN\to N\text{Δ}}\left(\sqrt{s}\right)=A_1+\frac{4A_2\times e^{-\left(\sqrt{s}-A_3\right)/A_4}}{{\left(1+e^{-\left(\sqrt{s}-A_3\right)/A_4}\right)}^2}\mathrm{,}\\ \sqrt{s}<2.21\text{ GeV}\mathrm{.}\end{array}$ (7) In which, A₁=-1.11 mb, A₂=26.30 mb, A₃=2.24 GeV, and A₄=0.05 GeV. We denote this as ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ to distinguish it from the default form in Ref. [40]. The fitting result is represented by red line in Fig. 2 (a). Above 2.21 GeV, the default fitting function was used.

Fig. 2Full-screen View

(Color online) (a) The $NN\to N\text{Δ}$ cross-section with default parameterization in the UrQMD model ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ (blue line), Hubbert parameterization ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ (red line), and that obtained based on the OBE model ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ [61] (orange line). (b) The ratio of ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ over ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ (red line) and ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ over ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ (orange line) as a function of energy $\sqrt{s}$

As shown in Fig. 2 (a), the ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ is closer to the experimental data than the ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$. The right panel shows that the ratio of $R={\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}/{\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$, and one can see that the cross-section ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ increases by a factor of 8.56 at a beam energy of 0.4A GeV. Consequently, one can expect a higher pion multiplicity with ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ than the one with ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$.

Additionally, we present another form of ${\sigma }_{NN\to N\text{Δ}}$ (indicated by the orange line in Fig. 2 (a)), which has been widely used in the other transport models [23, 27, 54, 30, 58]. This was obtained by fitting the results from the one-boson exchange model [61], in which the model parameters were obtained by fitting the CERN8301 experimental data [62], referred to as ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ in this work. At 0.4A GeV, the data from CERN8301 was smaller than the one from CERN8401, and the difference between the data from CERN8301 and CERN8401 was approximately 3 mb. Correspondingly, the ratio $R={\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}/{\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ decreased to 1.88, and one can expect that the pion multiplicity will be underestimated in transport models with this form. To enhance the pion multiplicity in the transport model calculations, the threshold effects may needed when the ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ is used in the model.

For the $N\text{Δ}\to NN$ cross-sections, they were obtained based on the detailed balance, in which the Δ mass-dependent $N\text{Δ}\to NN$ cross-sections were also considered as in Refs. [56, 28].

The descriptions of the collective flow and pion observables

The collective flow reflects the directional features of the transverse collective motion, which can be quantified in terms of the moments of the azimuthal angle relative to the reaction plane, i.e., $v_n=\langle \cos (n\varphi )\rangle$, $n=\mathrm{1,}\mathrm{2,}\mathrm{3,}\cdots$. Among the vn, the elliptic flow v₂ has been used to determine the MDI [63], and the elliptic flow ratios, such as $v_2^{\text{n}}/v_2^{\text{p}}$, $v_2^{\text{n}}/v_2^{\text{H}}$, and $v_2^{\text{n}}/v_2^{\text{ch}}$ are proposed to determine the symmetry energy at suprasaturation density [16, 31, 18]. Pions are known to be mainly produced through Δ resonance decay in the suprasaturation density region at an early stage, and the multiplicity ratio of charged pions, i.e., π^-/π⁺, is also considered as a probe for constraining the symmetry energy at the suprasaturation density and has been widely studied [19-23, 27, 28].

In this work, we perform the calculations of Au+Au collision at 0.4A GeV with UrQMD model, and 200,000 events are simulated for each impact parameter. The final flow observables as functions of p_T/A and $y_0^{\text{lab}}$ are obtained by integrating over b with a Gaussian weight. For the one with p_T/A, the integration range is b from 0 to 10 fm [31, 64-66]; for the one with $y_0^{\text{lab}}$, it is from 5 to 7 fm. The pion observable was also obtained by integrating over b from 0 to 2 fm with a certain weight, which is the same as that in Ref. [67]. The 17 observables listed in Table 3 are investigated in the following analysis.

3.1

Directed flow and elliptic flow

Figure 3 (a) and (b) show the directed flow as a function of the transverse momentum per particle p_T/A for neutrons $v_1^{\text{n}}\left(p_{\text{t}}/A\right)$ and for charged particles $v_1^{\text{ch}}\left(p_{\text{t}}/A\right)$ at given rapidity regions and angle cuts. The symbols represent the ASY-EOS data from Ref. [31]. The lines correspond to UrQMD calculations using $v_{\text{md}}^{\text{Arnold}}$ (green) and $v_{\text{md}}^{\text{Hama}}$ (blue and red) for the L = 144 MeV (solid lines) and L = 20 MeV (dashed lines) cases at S₀ = 32.5 MeV. In which, green and blue lines correspond to ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$, and red lines correspond to ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$. By comparing the green and blue lines, the effects of MDI can be understood, while the comparison between blue and red lines can be used to study the effects of ${\sigma }_{NN\to N\text{Δ}}$. The calculations show that $v_1^n\left(p_{\text{t}}/A\right)$ and $v_1^{\text{ch}}\left(p_{\text{t}}/A\right)$ increase from negative to positive as pt/A increases, and the sign of v₁ changes around $p_{\text{t}}/A\approx 0.5$ GeV/c. Furthermore, the calculations show that there are no sensitivities of v₁ to L, MDI, and ${\sigma }_{NN\to N\text{Δ}}$ in the selected rapidity region owing to the spectator matter-blocking effect. In addition, calculations with different combinations of L, MDI, and ${\sigma }_{NN\to N\text{Δ}}$ fall within this data region.

Fig. 3Full-screen View

(Color online) (a) $v_1\left(p_{\text{t}}/A\right)$ for neutrons; (b) $v_1\left(p_{\text{t}}/A\right)$ for charged particles; (c) $v_2\left(p_{\text{t}}/A\right)$ for neutrons, and (d) $v_2\left(p_{\text{t}}/A\right)$ for charged particles for ¹⁹⁷Au+¹⁹⁷Au collisions. The dash and solid lines correspond to the results with L=20 MeV and L=144 MeV at S₀=32.5 MeV, respectively. The gray lines are the results in Ref. [31]. The black symbols are ASY-EOS data[31]

Figure 3 (c) and (d) show the elliptic flow of neutrons $v_2^{\text{n}}\left(p_{\text{t}}/A\right)$ and charged particles $v_2^{\text{ch}}\left(p_{\text{t}}/A\right)$, with different L, MDI, and ${\sigma }_{NN\to N\text{Δ}}$. The symbols and lines have the same meaning as in panels (a) and (b), respectively. The gray lines represent the results in Ref. [31]. Both $v_2^{\text{n}}$ and $v_2^{\text{ch}}$ have negative values and decrease as p_T/A increases, indicating a preference for particle emission out of the reaction plane towards 90 and 270°. Notably, both $v_2^{\text{n}}$ and $v_2^{\text{ch}}$ at high p_t regions are highly sensitive to the strength of MDI and L but are hardly influenced by the forms of ${\sigma }_{NN\to N\text{Δ}}$. This is because only 6% of NN collisions belong to $NN\to N\text{Δ}$ collisions in the presently studied beam energy [28]. The values of v₂ obtained with $v_{\text{md}}^{\text{Hama}}$ are always lower than those obtained with $v_{\text{md}}^{\text{Arnold}}$ because the momentum dependence of $v_{\text{md}}^{\text{Hama}}$ is stronger than that of $v_{\text{md}}^{\text{Arnold}}$. The calculations of $v_2^{\text{n}}$ and $v_2^{\text{ch}}$ with $v_{\text{md}}^{\text{Hama}}$ are closer to the ASY-EOS experimental data than those obtained using the previous version of UrQMD [31] (gray lines), in which $v_{\text{md}}^{\text{Arnold}}$ is used.

In addition to MDI, the $v_2^{\text{n}}$ and $v_2^{\text{ch}}$ both exhibit some sensitivity to the stiffness of the symmetry energy. As shown in Fig. 3 (c), the values of $v_2^{\text{n}}$ obtained with the L=144 MeV (stiff) case are lower than those obtained with the L=20 MeV (soft) case. This is because the stiff symmetry energy provides a stronger repulsive force on the neutrons at suprasaturation density than that with the soft symmetry energy case. For charged particles, as shown in panel (d), the $v_2^{\text{ch}}$ obtained in the stiff symmetry energy case is higher than that obtained in the soft symmetry energy case. This is because the emitted charged particles are mainly composed of free protons, which feel a stronger attractive interaction for the stiff symmetry energy case than that for the soft symmetry energy case at suprasaturation density.

Figure 4 shows the calculated results on the elliptic flow of neutrons, protons, and H isotopes. Panels (a) and (b) show $v_2^{\text{n}}\left(p_{\text{t}}/A\right)$ and $v_2^{\text{H}}\left(p_{\text{t}}/A\right)$, and panels (c) and (d) show $v_2^{\text{n}}\left(y_0^{\text{lab}}\right)$ and $v_2^p\left(y_0^{\text{lab}}\right)$. The red lines have the same meaning as those shown in Fig. 3. The gray lines represent the calculations using $v_{\text{md}}^{\text{Arnold}}$ in the previous UrQMD model [14]. The black symbols represent elliptic flow data from the FOPI-LAND experiment [14]. The calculations with $v_{\text{md}}^{\text{Hama}}$ can nearly reproduce the FOPI-LAND data. However, the strength of $v_2^{\text{n}}$ and $v_2^{\text{H}}$ is slightly overestimated at p_T/A=0.85 GeV/c and an underestimation of the strength of $v_2^{\text{p}}$ at $y_0^{\text{lab}}>0.6$, which may be caused by isospin splitting of the proton and neutron effective masses [68, 69].

Fig. 4Full-screen View

(a) $v_2\left(p_{\text{t}}/A\right)$ for neutrons; (b) $v_2\left(p_{\text{t}}/A\right)$ for H isotopes; (c) $v_2\left(y_0^{\text{lab}}\right)$ for neutrons; and (d) $v_2\left(y_0^{\text{lab}}\right)$ for protons for ¹⁹⁷Au+¹⁹⁷Au collisions. The red lines are for $V_{\text{md}}^{\text{Hama}}$ and ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ at S₀=32.5 MeV. The gray lines are the results in Ref. [14], and the black symbols are the FOPI-LAND data[14]

To single out the contributions of the isovector potential and cancel those of the isoscalar potentials, $v_2^{\text{n}}/v_2^{\text{ch}}$, $v_2^{\text{n}}/v_2^{\text{H}}$, and $v_2^{\text{n}}/v_2^{\text{p}}$ ratios were proposed to probe the symmetry energy. Figure 5 (a) shows the calculations for $v_2^{\text{n}}/v_2^{\text{ch}}$ as a function of p_t/A obtained with $v_{\text{md}}^{\text{Hama}}$ and ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$. The lines represent the UrQMD calculations with L=20 MeV (dash) and L=144 MeV case (solid) at S₀ = 30 MeV (violet) and S₀ = 34 MeV (red). The lower calculation limit is L=5 MeV represented by the rectangular box. The calculations show that $v_2^{\text{n}}/v_2^{\text{ch}}$ is sensitive to L at the low p_t region where the mean field plays a more important role. The $v_2^{\text{n}}/v_2^{\text{ch}}$ values obtained with the stiff symmetry energy cases are larger than that with the soft symmetry energy case, which is consistent with the results by using the UrQMD model [14, 16, 31], IBUU model [17], or TüQMD model [18, 54]. This behavior can be understood from Fig. 3 (c) and (d). By comparing the calculations of $v_2^{\text{n}}/v_2^{\text{ch}}$ with the ASY-EOS experimental data [31] represented by the symbols and doing a χ² analysis, one can find the parameter sets favored by data. Within the framework of UrQMD and setting S₀= 30–34 MeV, the parameter sets with L = 5–70 MeV can describe the data.

Fig. 5Full-screen View

(Color online) (a) $v_2^{\text{n}}/v_2^{\text{ch}}$ as a function of p_t/A for ¹⁹⁷Au+¹⁹⁷Au collisions with L = 20 MeV and 144 MeV at S₀ = 30 MeV (violet line) and S₀ = 34 MeV (red line). The black symbols represent the ASY-EOS experimental data. (c) $v_2^{\text{n}}/v_2^{\text{H}}$ as a function of p_t/A, (e) $v_2^{\text{n}}/v_2^{\text{p}}$ as a function of $y_0^{\text{lab}}$ for ¹⁹⁷Au+¹⁹⁷Au collisions with L = 20 MeV and 144 MeV at S₀ = 32.5 MeV (red lines). The black symbols represent FOPI-LAND experimental data. (b) χ² of $v_2^{\text{n}}/v_2^{\text{ch}}$ as a function of L at S₀ = 30 MeV (violet line) and S₀ = 34 MeV (red line). (d), (f) χ² of $v_2^{\text{n}}/v_2^{\text{H}}$, $v_2^{\text{n}}/v_2^{\text{p}}$ as a function of L at S₀ = 32.5 MeV (red line)

Figure 5 (c) and (e) depict the calculated $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$ and $v_2^{\text{n}}/v_2^{\text{p}}\left(y_0^{\text{lab}}\right)$ results, respectively. The L dependence of the χ² value at S₀ = 32.5 MeV for $v_2^{\text{n}}/v_2^{\text{H}}$ and $v_2^{\text{n}}/v_2^{\text{p}}$ are presented in Fig. 5 (d) and (f). By comparing with the FOPI-LAND data [14] and conducting a χ² analysis, we obtain L constraints from $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$ in a broader range of L = 5-95 MeV, owing to the larger uncertainty in the $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$ data. The L constraints from $v_2^{\text{n}}/v_2^{\text{p}}\left(y_0^{\text{lab}}\right)$ calculations are L = 5–60 MeV, which is narrower than those from $v_2^{\text{n}}/v_2^{\text{ch}}\left(p_{\text{t}}/A\right)$ and $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$ since the more pronounced sensitivity of the symmetry potential effects on free protons. As our calculations can reproduce the differential data of elliptic flow, it naturally expects that the ratios of integral value of the elliptical flow, i.e., $v_2^{\text{n}}/v_2^{\text{H}}$, $v_2^{\text{n}}/v_2^{\text{p}}$, and the difference $v_2^{\text{n}}-v_2^{\text{H}}{v}_2^{\text{n}}-v_2^{\text{p}}$, can describe the experimental data.

The constraint by flow ratio and flow difference in this work is lower than those with the previous UrQMD model [31] or TüQMD model [27]. The discrepancy between our results and those of the previous UrQMD [31] is caused by using the different forms of v_md and K₀. In Ref. [31], $v_{\text{md}}^{\text{Arnold}}$ and K₀=200 MeV were used, and they provided a weaker repulsive force than that in our study. Consequently, the magnitude of the elliptical flow was underestimated. Consequently, their study required a more repulsive symmetry potential at high density, which has large L values. The difference between our results and TüQMD results [27] may be caused by using different K₀ and impact parameter ranges in calculations. In Ref. [27], they used K₀ = 214 MeV and the impact parameter b<7.5 fm. In our case, K₀=231 MeV, the impact parameter is distributed from 0 to 10 fm, and the weight of the impact parameter has a Gaussian form, inferred from the experimental event selection [31]. Both led to the L value constraints being greater than our study results. In addition, the different treatments for the medium effect of elastic cross-sections and isovector neutron-proton effective mass splitting may also have some effects; however, understanding the difference requires further study.

3.2

Pion productions and charged pion yield ratios

Figure 6 (a) shows the calculated pion multiplicity per participant Mπ/A_part as a function of L with different ${\sigma }_{NN\to N\text{Δ}}$ and symmetry energy forms. A_part is the number of nucleons in the participant, which constitutes 90% of the total mass of the system. The red and violet lines represent the calculations with ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ at S₀=30 and 34 MeV, respectively. The orange and blue lines represent the results obtained using ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ and ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ at S₀ = 32.5 MeV. The results obtained with ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ are lower than the data because the ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ is 50% lower than the ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ at E=0.4A GeV. In this case, one may expect that an obvious threshold effect is required to enhance the production of pions to describe the data. The calculation obtained using ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ underestimated Mπ/A_part by approximately 30%, relative to the data. This discrepancy can be understood from the underestimation of the $NN\to N\text{Δ}$ cross-sectional data using the default formula ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$, as shown in Fig. 2 (a). Surprisingly, the results obtained with the ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ fall into the data region since the ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ enhances the cross-sections by a factor of 8.56 at 0.4A GeV relative to the ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$. The above conclusion is not modified by using the different values of S₀ and L, i.e., in the range of S₀ = 30,34 MeV and L=5-144 MeV. Thus, Mπ/A_part alone cannot be used to distinguish between the different forms of symmetry energy.

Fig. 6Full-screen View

(Color online) (a) Mπ/A_part, and (b) π^-/π⁺ as a function of L for ¹⁹⁷Au+¹⁹⁷Au collisions with ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ at S₀=30 MeV (violet lines) and S₀=34 MeV (red lines), with ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ (blue line) and ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ (orange line) at S₀=32.5 MeV. The shaded region is the FOPI data [67]

In Fig. 6 (b), we present the calculated ratios π^-/π⁺ as a function of L with different forms of ${\sigma }_{NN\to N\text{Δ}}$ and S₀. These calculations indicate that π^-/π⁺ is sensitive to L and ${\sigma }_{NN\to N\text{Δ}}$. Using ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ or ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ leads to fewer pions, but increase the values of π^-/π⁺ which fall into the data region for L=5-144 MeV. Even calculations with ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ or ${\sigma }_{NN\to N\text{Δ}}^{\text{OBEM}}$ can reproduce the π^-/π⁺ data (the blue and orange lines); however, one cannot make this conclusion because the pion multiplicity is underestimated relative to the data. For the calculations with ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$, the data of both Mπ/A_part and π^-/π⁺ can be reproduced with the L=5-70 MeV and S₀=30-34 MeV parameter sets.

However, one should keep in mind that the integral observable, i.e., the pion multiplicity, is less influenced by pion potential due to the cancellation effects from pion potential and threshold effects [70]. To deeply understand the effect from pion-nucleon potential, a differential observable, such as the energy spectral of pion yields and charged pion rations or pionic flow, is suggested [71, 30] and it should be further studied in both the theoretical and experimental sides.

The symmetry energy constraints and its model dependence

4.2

The symmetry energy at characteristic densities and its model dependence

In Fig. 7 (a) and (b), we present the symmetry energy values at their characteristic densities, i.e., S(1.2ρ₀) and S(1.5ρ₀), obtained in this study (red symbols with errors). The uncertainties are the differences between the lower and upper boundaries of the favored symmetry energy parameter sets. The upper boundary of the symmetry energy was obtained with the symmetry energy with (S₀, L)=(34,70) MeV, and the lower boundary was obtained with the symmetry energy with (S₀, L)=(30,5) MeV. For the historical constraints on the symmetry energy [14-16, 21-23, 25-31, 70, 79], we calculate the density dependence of symmetry energy according to the constraints given in previous studies. Subsequently, the values of S(1.2ρ₀) and S(1.5ρ₀) and their uncertainties were similarly obtained. Th e results are represented by blue symbols with errors.

Fig. 7Full-screen View

(a) and (b) are the constraints of S(1.2ρ₀) and S(1.5ρ₀) in this study (red symbols) and from the previous study (blue symbols)

The S(1.2ρ₀) values obtained in this study were between 30 and 38 MeV. This is slightly lower than the constraints from the analyses of elliptic flow ratios or elliptic flow differences using the previous version of the UrQMD [16, 31] and TüQMD models [15, 18], which are in the 34-48 MeV range. The S(1.5ρ₀) value obtained in this study ranged from 28-44 MeV, which can overlap with the recent constraints by comparing SπRIT data with dcQMD [30] (S(1.5ρ₀)=38-72 MeV) and IBUU models [29] (S(1.5ρ₀)=35-47 MeV) within their uncertainties. Our results can also overlap with the previous constraints from the FOPI data by using the previous UrQMD version [28] by Liu et al, TüQMD [27] by Cozma et al, RVUU [70, 26] by Zhang et al, χBUU [79] by Zhang et al, pBUU [25] by Hong et al, and IBUU [21] by Xiao et al; however, they can overcome the constraints by using the isospin-dependent Boltzmann-Langevian (IBL) [23] model by Xie et al and the Lanzhou quantum molecular dynamics (LQMD) model [22] by Feng et al.

To quantitatively describe the theoretical uncertainties caused by the model dependence, a quantity, ${\delta }_{\text{model}}\left({\rho }^{*}\right)=\frac{S_{\text{max}}\left({\rho }^{*}\right)}{S_{\text{min}}\left({\rho }^{*}\right)}\mathrm{,}$ (10) is adopted. $S_{\text{max}}\left({\rho }^{*}\right)$ and $S_{\text{min}}\left({\rho }^{*}\right)$ are the largest and smallest values of the symmetry energy constraints at ρ^* among the different models. The larger the model dependence, the larger the δ_model. If there is no model dependence, δ_model will be one. For the symmetry energy constraints at 1.2ρ₀ using the FOPI-LAND and ASY-EOS flow data, that is, S(1.2ρ₀), the δ_model is 1.45. For the symmetry energy constraints at 1.5ρ₀ using the FOPI and SπRIT pion data, that is, S(1.5ρ₀), the δ_model is 2.75. These values are smaller than the model dependence described by the extrapolated symmetry energy at 3ρ₀, that is, S(3ρ₀), which is δ_model(3ρ₀) = 170. This clearly indicates that simply extrapolating the symmetry energy constraints from the characteristic density to other densities may lead to a misunderstanding of the symmetry energy constraints via HICs.

4.3

Remarks on the symmetry energy constraints at 0.1-3.0 ρ₀

Notably, presenting the symmetry energy only at 1.2ρ₀ and 1.5ρ₀ is incomplete because the probed density region using flow and pion observables is in a wide density region, that is, in 1.2±0.6ρ₀ for flow observables and 1.5± 0.5ρ₀ for pion observables. In Fig. 8 (a), we present the constrained symmetry energy in the flow characteristic density region (0.6-1.8ρ₀) as a pink shaded region, and the constraints in the pion characteristic density region (1.0–2.0ρ₀) with a violet shaded region. This completely overlaps with the constraints from the theoretical calculation using the chiral effective field theory (χEFT)[80] (green region); however, the uncertainty is larger than that from χEFT. Compared with the analyses of the SπRIT data obtained using dcQMD [30], the symmetry energy constraint in the high-density region is relatively small. However, it can overlap with the uncertainty.

Fig. 8Full-screen View

(Color online) (a) The density dependence of symmetry energy constraints in this study at 1.2± 0.6 ρ₀ region (wink region) and at the 1.5± 0.5 ρ₀ region (violet region). Other constraints are obtained from Ref. [72] (cyan region) and Ref.[80] (green region). (b) The density dependence of the corresponding pressure of neutron star matter in this study (violet region) and the extrapolation (violet dash lines), and the pressure constraints of neutron star matter obtained by Drischler et al.[80] (green region), Legred et al. [87] (pink region), and Huth et al. [73] (cyan region)

For symmetry energies below 0.6ρ₀ and above 2ρ₀, one can only infer the symmetry energy values by extrapolation because the symmetry energy information in these density regions is beyond the flow capability and pion observables at 0.4A GeV. The extrapolated symmetry energy below 0.6ρ₀ is consistent with the results from the neutron to proton yield ratios in HICs (HIC(n/p)) [81], the isospin diffusion in HICs (HIC(isodiff)) [82], the nuclear mass calculated by the Skyrme energy-density functional (Skyrme-EDF[A12])(Mass(Skyrme)) [83] and density functional theory (DFT) (Mass(DFT))[84], isobaric analog state (IAS) [85], and electric dipole polarization α_D[86], decoded by Lynch and Tsang in Ref. [72]. However, the uncertainties of the constraints using HICs in this study were larger than those of these observables. The extrapolated symmetry energy above 2ρ₀ is weaker than that obtained from the neutron star by Drischler et al. [80], Legred et al. [87], and Huth et al. [73], as shown in Fig. 8 (b). This discrepancy may be related to the momentum-dependent symmetry potential uncertainties, which may provide the same symmetry energy density dependence but with different effects on the isospin-sensitive observables [46, 88-93]. Thus, investigating the form of the momentum-dependent symmetry potential is very important in HICs.

Summary and outlook

In summary, we investigated the influence of different momentum-dependent interactions, symmetry energy and $NN\to N\text{Δ}$ cross-sections on nucleonic and pion observables, such as $v_1^{\text{n}}\left(p_{\text{t}}/A\right)$, $v_1^{\text{ch}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{ch}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{H}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}\left(y_0^{\text{lab}}\right)$, $v_2^{\text{p}}\left(y_0^{\text{lab}}\right)$, $v_2^{\text{n}}/v_2^{\text{ch}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}/v_2^{\text{p}}\left(y_0^{\text{lab}}\right)$, $v_2^{\text{n}}/v_2^{\text{H}}$, $v_2^{\text{n}}/v_2^{\text{p}}$, $v_2^{\text{n}}-v_2^{\text{H}}$, $v_2^{\text{n}}-v_2^{\text{p}}$, Mπ, and π^-/π⁺, using the UrQMD model for Au+Au collision at a beam energy of 0.4A GeV. Our results confirm that the elliptic flows of neutrons and charged particles, i.e., $v_2^{\text{n}}$ and $v_2^{\text{ch}}$, are sensitive to momentum-dependent interactions. The ASY-EOS and FOPI-LAND flow data favor calculations with strong momentum-dependent interactions, that is, $v_{\text{md}}^{\text{Hama}}$. However, calculations with $v_{\text{md}}^{\text{Hama}}$ underestimate the pion multiplicity by approximately 30% relative to FOPI data if the ${\sigma }_{NN\to N\text{Δ}}^{\text{UrQMD}}$ is adopted. Our calculations illustrate that the underestimation can be fixed by considering the accurate $NN\to N\text{Δ}$ cross-sections ${\sigma }_{NN\to N\text{Δ}}^{\text{Hub}}$ in the UrQMD model.

Furthermore, the symmetry energy constraints at the flow and pion characteristic densities were investigated using the updated UrQMD model. The characteristic density probed by the flow is approximately 1.2ρ₀, which is smaller than the pion characteristic density of 1.5ρ₀ [28]. By simultaneously describing the data of $v_2^{\text{n}}/v_2^{\text{ch}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}/v_2^{\text{H}}\left(p_{\text{t}}/A\right)$, $v_2^{\text{n}}/v_2^{\text{p}}\left(y_0^{\text{lab}}\right)$, $v_2^{\text{n}}/v_2^{\text{p}}$, $v_2^{\text{n}}/v_2^{\text{H}}$, $v_2^{\text{n}}-v_2^{\text{p}}$, $v_2^{\text{n}}-v_2^{\text{H}}$, and π^-/π⁺ with UrQMD calculations, the favored effective interaction parameter sets are obtained and we got the S(1.2ρ₀)=34± 4 MeV and S(1.5ρ₀)=36± 8 MeV. The extrapolated values of L in this work are in 5-70 MeV within 2σ uncertainty for S₀=30-34 MeV, which is below the analysis of the PREX-II results with a specific class of relativistic energy density functional [94], but is consistent with the constraint from the charged radius of ⁵⁴Ni [95], resulting from combining the astrophysical data with PREX-II and χEFT [96], and from the SπRIT pion data for the Sn+Sn collision at 0.27A GeV [30].

For the model dependence of the symmetry energy constraints, our calculations show that the strengths of the model dependence among the different transport models are 1.45 and 2.75 for the symmetry energy at the flow and pion characteristic density, respectively. These values are obviously smaller than the strength of the model dependence described by the symmetry energy at three times the normal density, which is 170.

Finally, simultaneously describing the ASY-EOS and FOPI data provides a rigorous limit on the UrQMD model and a solid foundation to further understand the effects of unsolved physics problems, such as the threshold effect, the pion potential, and the momentum-dependent symmetry potential.

Notably, the discrepancies in $v_2^{\text{n}}$ and $v_2^{\text{p}}$ at high p_t and rapidity relative to the data demonstrate the importance of the momentum dependence of the symmetry potential, as mentioned in Refs. [68, 77, 69], which should be investigated using the momentum and rapidity distributions of the nucleonic and pionic probes in the future. Another important direction for developing the transport model and limiting its uncertainties is describing the nucleonic and pionic flow observables and their spectra at subthreshold energies and above 1 GeV/u. This will help to further understand the pion production mechanism and provide the symmetry energy constraints twice beyond the normal density with HICs.