Progress of quantum molecular dynamics model and its applications in heavy ion collisions

Abstract

In this review article, we first briefly introduce the transport theory and quantum molecular dynamics model applied in the study of the heavy ion collisions from low to intermediate energies. The developments of improved quantum molecular dynamics model (ImQMD) and ultra-relativistic quantum molecular dynamics model (UrQMD), are reviewed. The reaction mechanism and phenomena related to the fusion, multinucleon transfer, fragmentation, collective flow and particle production are reviewed and discussed within the framework of the two models. The constraints on the isospin asymmetric nuclear equation of state and in-medium nucleon–nucleon cross sections by comparing the heavy ion collision data with transport models calculations in last decades are also discussed, and the uncertainties of these constraints are analyzed as well. Finally, we discuss the future direction of the development of the transport models for improving the understanding of the reaction mechanism, the descriptions of various observables, the constraint on the nuclear equation of state, as well as for the constraint on in-medium nucleon–nucleon cross sections.

Graphical abstract

Keywords

quantum molecular dynamics model / low energy heavy ion collisions / low-intermediate energy heavy ion collisions / fusion / multinucleon transfer reaction / multifragmentation / collective flow / isospin asymmetric equation of state / in-medium nucleon–nucleon cross sections

1 Introduction

Heavy ion collisions (HICs) provide a unique way to investigate the basic nuclear physics problems in laboratory, such as where is the end of nuclear chart? how can we reach it? what is the properties of neutron-rich matter in cosmos and on the earth? and the origin of elements heavier than Fe? In the past 50 years, a large number of experimental data for HICs have been accumulated and lots of interesting results have been obtained [1-24]. It has been found that the phenomena and the physics processes in HICs at various energy domains are very rich and complicated. Thus, understanding the mechanism behind them is very helpful for us to obtain the knowledge of related nuclear phenomena and physics.

For low energy heavy ion reactions (in this review paper, we consider this energy domain to be E_{beam}≤ 10 MeV/nucleon), the synthesis of superheavy nuclei and new neutron (proton)-rich isotopes are the hot topics [25-61]. Related to them, there are also a lot of interesting problems in the reactions, such as the role of the dynamical effect, the nuclear structure effect and the interplay between them, and the dependence of the barrier on the reaction systems and incident energies, need to be clarified. Multi-nucleon transfer (MNT) reactions between very heavy nuclei at near or above barrier energies, attracted a lot of attention as well [52-54] since it could lead to production of new heavy and possibly superheavy neutron-rich isotopes, especially those nuclei at waiting point for r-process which are important to the nuclear astrophysics. In addition, the study of MNT reaction is highly required due to the new facilities of High Intensity heavy ion Accelerator Facility (HIAF) and Argonne Tandem Linac Accelerator System (ATLAS).

As the beam energy increases, the interplay of one-body mean field and two-body scattering makes the phenomena in this energy domain fascinating and complicated. The ternary and four-fragment breakup process may take place accompanying the binary process. At the low-intermediate energies region, i.e., ~20 MeV/nucleon to ~300 MeV/nucleon (the definition of energy domain of intermediate energies, usually differs among authors), multifragmentation process which was thought as a signal of liquid–gas phase transition in finite systems, appears and provides the opportunity to study the equation of state at sub-saturated densities [22, 23, 62-65] and mechanism of liquid-gas phase transition in the finite system. The phenomena relating to the multifragmentation, such as isospin distillation, neutron-rich neck emission, bimodality, …, evolved with the isospin asymmetry and the size of reaction system and the beam energy as well, are very rich and complicated, and it provides us with the hints of equation of state (EOS), mechanism of fragmentation and information of liquid-gas phase transition for two components finite system.

When the beam energies increase from around 100 MeV per nucleon up to few GeV per nucleon, the collective motion of emitted nucleons and light particles appears, and it was named as collective flow [3, 9, 11-14, 16]. The measurement of collective flow also provides a possibility for the study of the nuclear EOS at supersaturated densities as well as the information of in-medium nucleon–nucleon (NN) cross sections [24], because the changing of the momentum of emitted nucleons is related to the gradient of pressure between the participant and spectator of the reaction system, which is caused by the nucleonic mean field and nucleon–nucleon collisions. The sub-threshold and above threshold production of mesons take place in addition to the multifragmentation and the collective flow effects. The yields of mesons and their yield ratios between different charge states as well as the flow effects evolve with energy and isospin asymmetry of the system [66-84], contain a lot of information of the isospin asymmetric EOS of nuclear matter, the isospin dependent in-medium NN cross sections and the reaction dynamics. The study of the phenomena at the beam energy from low to intermediate energies, such as the production of new isotopes, multifragmentation, and flow effects closely relating to equation of state stimulate the building of next generation of rare isotope facilities, such as the Intensity Heavy ion Accelerator Facility (HIAF) at Institute of Modern Physics (IMP) in China [85-87], Facility for rare isotope beams (FRIB) [88] in USA, Rare isotope Accelerator complex for ON-line experiment (RAON) [89] in Korea, RI beam factory (RIBF/RIKEN) in Japan [90], SPIRAL2/GANIL in France [91], Facility for Antiproton and Ion Research (FAIR) at GSI in Germany [92], Selective Production of Exotic Species (SPES/LNL) in Italy [93], Nuclotron-based Ion Collider facility at the Joint Institute for Nuclear Research(NICA/JINR) in Russia [94], and the proposed Beijing Isotope-Separation-On-Line Neutron-Rich Beam Facility (BISOL) [95] in China.

In order to describe the phenomena and processes in HICs at each energy domain, kinds of successful models have been developed for explaining and investigating the related physics problems. For instance, for low beam energy reactions, both the diffusion model by using Langevin equation and the di-nuclear model can describe fusion cross sections [46, 52, 53, 53, 54, 96-117]; for the intermediate energy region, the statistical multifragmentation model (SMM) is quite successful in describing the charge distribution of multifragmentations at intermediate energy reactions [118-125], etc. However, heavy ion reaction or collision process is essentially a non-equilibrium process. It is highly demanded to develop a model or theory which can describe the time evolution of reaction, so that one can understand the reaction dynamics and the transient states in the compressing and expanding stages which are important for studying the equation of state of nuclear matter. Also, the theories (models) describing heavy ion collisions in a more unified way from low energies to higher energies are considerably demanded for insight in depth of the properties of nuclear many-body system and reaction dynamics. Thus, the microscopic dynamics models, namely, the microscopic transport models for heavy ion reactions were developed and have been used to study HICs with a great success [63-65, 73, 74, 126-174].

Generally, the transport models adopt different philosophy to solve the transport equations, which can be roughly divided into two types, the Boltzmann–Uehling–Uhlenbeck (BUU) model and the quantum molecular dynamics (QMD) model.

In the BUU approach, the goal is to describe the evolution of the one-body phase space occupation probability f(r, p; t) as a function of the time under the action of a mean field potential U[f], usually derived from a density functional, and two-body collisions specified by an in-medium cross section \text{d}{\sigma }^{\text{med}}/\text{d}\text{Ω}. The non-relativistic BUU equation reads

\begin{array}{}\text{(1)}& \left(\frac{\partial }{\partial t}+\frac{p}{m}\cdot {\nabla }_r-{\nabla }_{r}U\cdot {\nabla }_p\right)f(r\mathrm{,}p;t)=I_{coll}(r\mathrm{,}p;t)\mathrm{.}\end{array}

I_{coll} is the collision term, which accounts for the particles entering and leaving the phase space {\text{d}}^{3}r{\text{d}}^{3}p. The integro-differential non-linear BUU equation is solved numerically by using the test particle method (TP) [175], as

\begin{array}{}\text{(2)}& f(r\mathrm{,}p;t)=\frac{{(2\pi \hslash )}^3}{g{N}_{\text{TP}}}{\displaystyle \sum _{i=1}^{A{N}_{\text{TP}}}}G(r-r_i(t))\tilde{G}(p-p_i(t))\mathrm{,}\end{array}

where N_{\text{TP}} is the number of test particles (TP) per nucleon, r_i and p_i are the time-dependent coordinate and momentum of the test particle i, G and \tilde{G} are the shape functions in coordinate and momentum space, respectively, with a unit norm (e.g., δ-functions or normalized Gaussians). For nucleons, the degeneracy factor g=4 is to define f(r, p; t) as the spin–isospin averaged phase space occupation probability. It is also possible to express the distribution function for each isospin (or spin) state in a similar way. In the BUU approach the phase space distribution function is seen as a one-body quantity and a smooth function of coordinates and momenta and it can be approximated better by increasing the number of test particles in the solution. In the limit of N_{\text{TP}}\to \infty, the BUU equation is solved exactly. In this limit the solution is deterministic and does not contain fluctuations. In case of studying the cluster and fragment production, it will need to introduce the fluctuation which actually exist in low-intermediate energy heavy ion collisions. There are some effort to do it through the Boltzmann–Langevin equation (BLE) which adds a fluctuation term on the right hand side of Eq. (1),

\begin{array}{}\text{(3)}& \begin{array}{l}\left(\frac{\partial }{\partial t}+\frac{p}{m}\cdot {\nabla }_r-{\nabla }_{r}U\cdot {\nabla }_p\right)f(r\mathrm{,}p;t)\\ =I_{coll}(r\mathrm{,}p;t)+\delta I_{coll}\mathrm{.}\end{array}\end{array}

This equation is solved approximately in the stochastic mean field (SMF), Boltzmann–Langevian one body (BLOB) [126, 127, 140-142, 176-182].

In the QMD approach, the evolution of N-body phase space density distribution f_N(r_1\mathrm{,}\dots \mathrm{,}{r}_N\mathrm{,}{p}_1\mathrm{,}\dots \mathrm{,}{p}_N;t) is formulated. Thus, in the philosophy of QMD, the effect of go-beyond the mean field approach is realized by including correlations and fluctuations from the beginning, and it is at the expense of more rapidly destroying the fermionic character of the system and of reverting to a classical system. The QMD model can be seen as an event generator, where the time evolution of different events is solved independently. The fluctuations in QMD-type codes are regulated and smoothed by choosing the parameter {\sigma }_r, the width of the wave packet. The mean field part in the QMD approach can also be viewed as derived from the Time-Dependent Hartree method with a product of trial wave function of single particle states in Gaussian form that we will mention in Section 2. Also the collision term, which relocates nucleon wave packets in momentum space, introduces more fluctuations than those for the collision term in BUU. The fluctuations among events are not suppressed in QMD approach even in the limit of infinite number of events. Taking into account all fluctuations and correlations has basically two advantages: i) many-body processes, in particular the formation of complex fragments are explicitly treated, and ii) the model allows for an event-by-event analysis of heavy ion reactions similar to the methods which are used for the analysis of exclusive high acceptance data [148-150].

Stimulated by studying the dynamical effects on the heavy ion reactions near the barrier, we have made a series of improvements on nucleon propagation in the mean field part, nucleon–nucleon collision in two-body collision part, the initialization and the cluster recognization [155, 156, 159, 162] based on the original QMD model during last 20 years. It was named as improved quantum molecular dynamics model (also known as ImQMD). One of the important changes is that we adopt the potential energy density functional to determine the nucleon propagation, which was stimulated by a steady transition taken place during the past several years from the mean-field description of nuclear properties in terms of effective forces to an density functional approach [183-187]. The basic idea of density functional is that the ground-state energy of a stationary many-body system can be represented in terms of the ground state density alone, and thus, energy density functional theory calculations are comparatively simple to implement yet often very accurate and computationally feasible even for systems with large particle numbers. Following this transition, we replaced the mean filed potential part in terms of the effective nuclear interaction by the potential energy density functional, and we will mention it in the following chapter.

Up to now, we have developed three versions of ImQMD designed for different purposes. The details will be given in Section 2. A series of applications of the ImQMD model to the fusion reaction, deep inelastic scattering (multi-nucleon transfer reaction) at near Coulomb barriers and to multifragmentation, collective flow effects and other important phenomena at intermediate and high energies will be presented in this review paper. The study of liquid-gas phase transition in finite nuclear systems connecting to multifragmentation, and the constraints on symmetry energy for asymmetric nuclear matter by comparison of the model calculations with the experimental data, are also presented. Further more, the phenomena at intermediate and high energies related to particle production and collective flow are presented and discussed within the framework of UrQMD model.

This article is organized as follows. In Section 2, we will briefly review the many-body transport theory and its solution in the quantum molecular dynamics approaches. In Section 3, we will focus on the study of phenomena in heavy ion reaction (collision), and its reaction mechanism. Section 4 will present the investigation on the in-medium NN cross sections by using the closed time-path green function method and one-boson exchange model, and the efforts of extracting the in-medium NN cross section by comparing the QMD type model calculations to heavy ion collisions data. Section 5 will present the investigation on the isospin asymmetric nuclear equation of state, symmetry energy and the uncertainties of these constraints. Section 6 will give discussions and prospect on the development of transport theory and its applications in the nuclear physics in the future.

2 Transport theory and the quantum molecular dynamics model

Utility of the quantum–mechanical phase space distributions for the formulation and solution of scattering and production problems is an important approach in heavy ion collisions. In 1932, Wigner [188] discovered an interesting version of the density matrix which allows the expression of quantum dynamics in a form directly comparable with the classical analog while maintaining the quantum integrity of the conjugacy of the variables p and x. This approach can also be extended to the particle production problem by expanding the field operator φ in terms of the annilation operator a(p\mathrm{,}t) [189], and the N-particle covariant distribution functions are directly connected with the inclusive differential cross sections. The advantages of this technique are clearly exhibited in application to quantum transport theory [190-192].

The transport equation can be derived by means of the nonequilibrium Green’s function technique, i.e., the closed time-path Green’s function(CTPGF) technique which is based on the theoretical concepts for a proper many-body description in terms of a real time nonequilibrium field theory initiated by Schwinger in the early sixties. By using the CTPGF technique, the transport equation and the analytical expressions of in-medium two-body scattering cross sections applicable for heavy ion collisions are simultaneously obtained. As the whole theory is complicated but very useful, here we refer to the literatures (see Refs. [193-201]).

In this section, we will shortly mention the transport theory which is used for QMD approach.

2.1 Transport theory for N-body system

The derivation of Boltzmann’s N-body phase space distribution f_N(r_1\mathrm{,}\dots \mathrm{,}{r}_N;p_1\mathrm{,}\dots \mathrm{,}{p}_N) provide an intuitive picture of complex collisions process based on the quantum mechanics. As is known that the simultaneous probability for position r and momentum p is forbidden in quantum mechanics by Heisenberg’s uncertainty principle, i.e., \text{Δ}{r}_i\text{Δ}{p}_i⩾\hslash /2. One way to map a quantum variable to a classical one is to use coarse-grain method, such as Wigner transformation [188] in which both respects the rules of quantum mechanics and recaptures most of the desired features of Boltzmann function. If a wave function \psi (r_1\mathrm{,}\dots \mathrm{,}{r}_N) is given, one may build the density matrix with the representation as following:

\begin{array}{}\text{(4)}& \begin{array}{l}{f}_N(r_1\mathrm{,}\cdots \mathrm{,}{r}_N;p_1\mathrm{,}\cdots \mathrm{,}{p}_N)\\ ={\left(\frac{1}{2\pi \hslash }\right)}^N{\displaystyle {\int }_{-\infty }^{\infty }}\cdots {\displaystyle \int }\text{d}{y}_1\cdots \text{d}{y}_N\\ \cdot {\psi }^{*}(r_1-y_1\mathrm{,}\cdots \mathrm{,}{r}_N-y_N)\\ \cdot \psi (r_1+y_1\mathrm{,}\cdots \mathrm{,}{r}_N+y_N){\text{e}}^{-2\text{i}(p_1\cdot y_1+\cdots +p_N\cdot y_N)/\hslash }\mathrm{,}\end{array}\end{array}

to express the probability-function of the simultaneous values of r_1\mathrm{,}\dots \mathrm{,}{r}_N and p_1\mathrm{,}\dots \mathrm{,}{p}_N. Eq. (4) is real, but not everywhere positive and it means that the f_N(r_1\mathrm{,}\dots \mathrm{,}{r}_N;p_1\mathrm{,}\dots \mathrm{,}{p}_N) cannot be really interpreted as the simultaneous probability for coordinates and momenta. However, the lack of positivity will not hinder the use of it since we are mainly concerned with computation of positive definite asymptotic quantities. When Eq. (4) is integrated with respect to p, the correct probabilities {\left|\psi (r_1\mathrm{,}\dots \mathrm{,}{r}_N)\right|}^2 is given; if we integrate Eq. (4) with respect to r, the correct probabilities {\left|C(p_1\mathrm{,}\dots \mathrm{,}{p}_N)\right|}^2 can also be verified. Hence, one may get the correct expectation values of mechanical quantities O as a function of coordinate or momenta for the state ψ as

\begin{array}{}\text{(5)}& \begin{array}{lll}\langle O\rangle \hfill & =\hfill & {\displaystyle {\int }_{-\infty }^{\infty }}\cdots {\displaystyle \int }\text{d}{r}_1\cdots \text{d}{r}_N\text{d}{p}_1\cdots \text{d}{p}_{N}O(r_i\mathrm{,}{p}_i)\hfill \\ \hfill & \hfill & \cdot f_N(r_1\mathrm{,}\cdots \mathrm{,}{r}_N;p_1\mathrm{,}\cdots \mathrm{,}{p}_N)\hfill \\ \hfill & =\hfill & {\displaystyle {\int }_{-\infty }^{\infty }}\cdots {\displaystyle \int }\text{d}{r}_1\cdots \text{d}{r}_N{\psi }^{*}(r_1\mathrm{,}\cdots \mathrm{,}{r}_N)\hfill \\ \hfill & \hfill & \cdot O(r_i\mathrm{,}-\text{i}\hslash {\partial }_{r_i})\psi (r_1\mathrm{,}\cdots \mathrm{,}{r}_N)\mathrm{,}\hfill \end{array}\end{array}

where f_N(r_1\mathrm{,}\dots \mathrm{,}{r}_N;p_1\mathrm{,}\dots \mathrm{,}{p}_N) is the probability-function described above.

For the equation of N-body system phase space distribution f_N(r_1\mathrm{,}{r}_2\mathrm{,}\dots \mathrm{,}{r}_N;p_1\mathrm{,}{p}_2\mathrm{,}\dots \mathrm{,}{p}_N), it has been derived by Aichline in Ref. [150] and is written as follows:

\begin{array}{}\text{(6)}& \begin{array}{l}\left(\frac{\partial }{\partial t}+{\displaystyle \sum _i}\frac{p_i}{m}\cdot {\nabla }_{r_i}\right)f_N(r_1\mathrm{,}\cdots \mathrm{,}{r}_N\mathrm{,}{p}_1\mathrm{,}\cdots \mathrm{,}{p}_N;t)\\ ={\displaystyle \int }{\text{Π}}_i{\text{d}}^3{p}_i{\text{d}}^3{Q}_i{\text{d}}^3{q}_i{\text{e}}^{\text{i}{r}_i\cdot p_i}\\ \cdot f_0^{(n)}(Q_1\mathrm{,}\dots \mathrm{,}{Q}_N\mathrm{,}{q}_1\mathrm{,}\dots \mathrm{,}{q}_N;t)\\ \cdot [I_1(T)+I_2(T)+I_3(T)]\mathrm{,}\end{array}\end{array}

where Q_i and q_i are the momenta of final states. T={\displaystyle \sum }{t}_i+{\displaystyle {\sum }_{k\ne m}}{\displaystyle {\sum }_m}{t}_{im}{G}_0^{†}{t}_{ik}+\dots, t_{im} is the sum of all possible transition matrix combinations, and G_0^{†} are the on shell propagators. f_0^{(n)} is the time evolved free wave packet. The definition of terms I₁, I₂, and I₃ can be found in Ref. [150]. The real part of I_1+I_2 acting as an effective potential has been replaced by two-body potential in this case, and it can be easily related to the nuclear equation of state. The I₃ term can be reduced to a sum of terms which contains only absolute squares of transition matrices, and it is assumed to be proportional to products of the cross sections [150].

For the Vlasov equation of Eq. (6), i.e., without I₃ term, the time evolution of the phase space density of particles which moves on classical orbits is specified by the Hamilton equation. The equations of motion for the expectation values of \langle r_i\rangle and \langle p_i\rangle are

\begin{array}{}\text{(7)}& \begin{array}{l}\frac{\partial \langle p_i\rangle }{\partial t}=-\langle {\nabla }_{i}V(r_1\mathrm{,}\cdots \mathrm{,}{r}_N)\rangle \mathrm{,}\\ \frac{\partial \langle r_i\rangle }{\partial t}=\frac{\langle p_i\rangle }{m}\mathrm{.}\end{array}\end{array}

Supposing the potential V(r_1\mathrm{,}\dots \mathrm{,}{r}_N) is slowly varying, and we can expand \frac{\partial V}{\partial r_i} as a Taylor series about \langle r_i\rangle as following:

\begin{array}{}\text{(8)}& \begin{array}{l}\frac{\partial V(r_1\mathrm{,}\cdots \mathrm{,}{r}_N)}{\partial r_i}\\ =\frac{\partial V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)}{\partial \langle r_i\rangle }\\ +\frac{{\partial }^{2}V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)}{\partial {\langle r_i\rangle }^2}(r_i-\langle r_i\rangle )\\ +\frac{1}{2}\frac{{\partial }^{3}V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)}{\partial {\langle r_i\rangle }^3}{\left(r_i-\langle r_i\rangle \right)}^2+\cdots \mathrm{.}\end{array}\end{array}

Substitution of the above expansion into Eq. (7), and we have

\begin{array}{}\text{(9)}& \begin{array}{l}\frac{\partial \langle p_i\rangle }{\partial t}\\ =\frac{\partial \langle V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)\rangle }{\partial \langle r_i\rangle }\\ +\frac{{\sigma }^2}{2}\frac{{\partial }^3\langle V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)\rangle }{\partial {\langle r_i\rangle }^3}+\cdots \mathrm{,}\end{array}\end{array}

with \langle (r_i-\langle r_i\rangle )\rangle =0 and {\sigma }^2=\langle {(r_i-\langle r_i\rangle )}^2\rangle. If

\begin{array}{}\text{(10)}& \begin{array}{l}\left|\frac{1}{2}{\nabla }_{\langle r_i\rangle }^{3}V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N){\sigma }^2\right|\\ \ll \left|{\nabla }_{\langle r_i\rangle }V(r_1\mathrm{,}\cdots \mathrm{,}\langle r_i\rangle \mathrm{,}\cdots \mathrm{,}{r}_N)\right|\mathrm{,}\end{array}\end{array}

the second term in the right hand of Eq. (9) can be neglected. While, \langle {\nabla }_{\langle r_i\rangle }V(r_1\mathrm{,}\dots \mathrm{,}\langle r_i\rangle \mathrm{,}\dots \mathrm{,}{r}_N)\rangle \approx {\nabla }_{\langle r_i\rangle }U(\langle r_1\rangle \mathrm{,}\dots \mathrm{,}\langle r_N\rangle ) and the equations in Eq. (7) will be the same as it in the classical equation. Here, the potential energy U is a function of parameters {\langle r_i\rangle, \langle p_i\rangle} and thus U can also be thought as the potential at {\langle r_i\rangle, \langle p_i\rangle} in practical calculations. If the gradients of the potential are strong, the high order terms in Eq. (7) cannot be neglected and it causes not only strong force but also large fluctuations around the mean trajectories.

For the collision term, in the actual calculations, the particle collisions are simulated when the particles are sufficiently close [202] with the Monte-Carlo method, and the scattering angle is chosen randomly according to the differential cross section. The outgoing states of collisions also need to be checked whether the states of outgoing particles have been occupied or how much is the probability of the occupation by other particles. If the outgoing states have been fully occupied, the collision will not happen and it is named as Pauli blocking. Otherwise, the collision will happen with certain probability. More details will be given in following sections.

2.2 Quantum molecular dynamics approach

In the quantum molecular dynamics approach, each nucleon is represented by a Gaussian wave packet,

\begin{array}{}\text{(11)}& {\psi }_i(r_i)=\frac{1}{{(2\pi {\sigma }_r^2)}^{3/4}}{\text{e}}^{-\frac{{(r_i-r_{i0})}^2}{2{\sigma }_r^2}+\text{i}(r_i-r_{i0})\cdot p_{i0}/\hslash }\mathrm{,}\end{array}

here, {\sigma }_r and r_{i0} are the width and centroid of wave packet, respectively. Its Wigner density reads

\begin{array}{}\text{(12)}& \begin{array}{l}{f}_i(r\mathrm{,}p;t)\\ =\frac{1}{{(2\pi {\sigma }_r^2)}^{3/2}}{\text{e}}^{-{(r-r_{i0})}^2/(2{\sigma }_r^2)}\frac{1}{{(2\pi {\sigma }_p^2)}^{3/2}}{\text{e}}^{-{(p-p_{i0})}^2/(2{\sigma }_p^2)}\\ =\frac{1}{{(\pi \hslash )}^3}\exp \left[-\frac{{(r_i-r_{i0})}^2}{2{\sigma }_r^2}-\frac{{(p_i-p_{i0})}^2}{2{\sigma }_p^2}\right]\mathrm{,}\end{array}\end{array}

where {\sigma }_r{\sigma }_p=\hslash /2. The features of f_i from Gaussian wave function are as follows: (i) f_i does not spread with time for fixed p or r, even though the underlying wave function does. But it can reproduce the fine structure of nucleon distance in fragmentation. (ii) An equation f(r\mathrm{,}p;t)=\delta (r-r_{i0})\delta (p-p_{i0}) is out of reach, because the uncertainty relation \text{Δ}x\text{Δ}{p}_x⩾\frac{\hslash }{2}. (iii) As {\sigma }_r tends to infinity, f_i spread uniformly over space, while the momentum factor is sharp in momentum.

Since the nuclei of colliding is a N-body system, its wave function should be \psi (r_1\mathrm{,}\dots \mathrm{,}{r}_N). In the QMD approach, the system wave function is assumed as a direct product of N coherent states, which is in the Hartree approximation,

\begin{array}{}\text{(13)}& \psi (r_1\mathrm{,}\cdots \mathrm{,}{r}_N)={\varphi }_{k_1}(r_1){\varphi }_{k_2}(r_2)\cdots {\varphi }_{k_N}(r_N)\mathrm{,}\end{array}

\begin{array}{}\text{(14)}& \begin{array}{ll}{\varphi }_{k_i}(r_i)=\hfill & \frac{1}{{(2\pi {\sigma }_r^2)}^{3/4}}\hfill \\ \hfill & \cdot \exp \left[-\frac{{(r_i-r_{i0})}^2}{2{\sigma }_r^2}+\frac{\text{i}{p}_{i0}\cdot (r_i-r_{i0})}{\hslash }\right]\mathrm{,}\hfill \end{array}\end{array}

where {\varphi }_{k_i}(r_i) is the wave function of the ith particle at state k_i (p_i=\hslash k_i form). {\varphi }_{k_i} is chosen as Gaussian wave packet to avoid the negative values of phase space distribution (f_N).

For the N-body Wigner function within the QMD assumption, it reads

\begin{array}{}\text{(15)}& \begin{array}{l}{f}_N(r_1\mathrm{,}\cdots \mathrm{,}{r}_N;p_1\mathrm{,}\cdots \mathrm{,}{p}_N)={\displaystyle \prod _{i=1}^N}f(r_i\mathrm{,}{p}_i)\\ ={\displaystyle \prod _{i=1}^N}\frac{1}{{(\pi \hslash )}^3}\exp \left[-\frac{{(r_i-r_{i0})}^2}{2{\sigma }_r^2}-\frac{{(p_i-p_{i0})}^2}{2{\sigma }_p^2}\right]\mathrm{,}\end{array}\end{array}

r_{i0}=\langle r_i\rangle and p_{i0}=\langle p_i\rangle are the centroids of wave packets in coordinate and momentum space, respectively. The time evolution of the f_N, i.e., Eq. (6), is a highly non-linear integral-differential equation and is difficult to be solved exactly.

Since the width of wave packet is fixed during the time evolution in the QMD approcah, the time evolution of the phase space density can be determined from the time evolution of the centroids of the wavepacket in the coordinate and momentum spaces, which are driven by the mean field potential and nucleon–nucleon collisions. Two ingredients in Eq. (6), such as the mean field part and the collision part, are solved separately rather than self consistently in practical calculations. The initialization is also very important for simulating the heavy ion collisions, and we will mention it in the next section.

In the Vlasov model (i.e., only the mean field and without collision part), the time evolution of the centroid of wave packets in the coordinate and momentum space is derived in previous section and is written as

\begin{array}{}\text{(16)}& \frac{\partial \langle p_i\rangle }{\partial t}\approx -\frac{\partial U(r_{10}\mathrm{,}\cdots \mathrm{,}{r}_{N0})}{\partial r_{i0}}\mathrm{,}\end{array}

\begin{array}{}\text{(17)}& \frac{\partial \langle r_i\rangle }{\partial t}=\frac{p_{i0}}{m}\mathrm{.}\end{array}

The time evolution of the centroids of p_{i0} and r_{i0} has also been derived by using the Euler–Lagrange equations as in Ref. [150],

\begin{array}{}\text{(18)}& \frac{\text{d}}{\text{d}t}\frac{\partial \mathcal{L}}{\partial {\dot{p}}_{i0}}-\frac{\partial \mathcal{L}}{\partial p_{i0}}=0\to {\dot{r}}_{i0}=\frac{p_{i0}}{m}+{\nabla }_p{U}_i\mathrm{,}\end{array}

\begin{array}{}\text{(19)}& \frac{\text{d}}{\text{d}t}\frac{\partial \mathcal{L}}{\partial {\dot{r}}_{i0}}-\frac{\partial \mathcal{L}}{\partial r_{i0}}=0\to {\dot{p}}_{i0}=-{\nabla }_p{U}_i\mathrm{.}\end{array}

The potential energy U in the quantum molecular dynamics model can be directly calculated from the potential operator \widehat{V}=v_{ij}+v_{ijk}+\cdots as follows based on the above assumption,

\text{d}{\text{Γ}}_i={\text{d}}^3{r}_i{\text{d}}^3{p}_i, v_{ij}, v_{ijk} are two-body interaction, three-body interaction, respectively. As in Ref. [150], the interactions may consist of local interaction, Yukawa and Coulomb interactions.

The local interaction has the form

\begin{array}{}\text{(21)}& v_{ij}=t_1\delta (r_1-r_2)\mathrm{,}{v}_{ijk}=t_2\delta (r_1-r_2)\delta (r_1-r_3)\mathrm{,}\end{array}

and one has

\begin{array}{}\text{(22)}& U_{ij}=t_1\tilde{\rho }(r_{i0}\mathrm{,}{r}_{j0})=\frac{t_1}{{(4\pi {\sigma }_r^2)}^{3/2}}{\text{e}}^{-{(r_{i0}-r_{j0})}^2/(4{\sigma }_r^2)}\mathrm{,}\end{array}

\begin{array}{}\text{(23)}& \begin{array}{l}{U}_{ijk}=\frac{t_2}{{(2\pi {\sigma }_r^2)}^3\cdot 3^{3/2}}\\ \cdot {\text{e}}^{-[(r_{i0}-r_{j0}{)}^2+{(r_{i0}-r_{k0})}^2+{(r_{k0}-r_{j0})}^2]/(6{\sigma }_r^2)}\\ \approx \frac{t_2}{{(2\pi {\sigma }_r^2)}^3\cdot 3^{3/2}}{\text{e}}^{-[(r_{i0}-r_{j0}{)}^2+{(r_{i0}-r_{k0})}^2]/(4{\sigma }_r^2)}\mathrm{.}\end{array}\end{array}

In the spin saturated nuclear matter, the three-body interaction can be viewed as the density dependent two-body interaction due to the hard core, and one can replace the three-body interaction to its effective form

\begin{array}{}\text{(24)}& t_2\delta (r_1-r_2)\delta (r_1-r_3)=\frac{t_2}{6}\delta (r_1-r_2)\rho \left(\frac{r_1+r_2}{2}\right)\mathrm{.}\end{array}

The Yukawa interaction is as

\begin{array}{}\text{(25)}& V^{Yuk}=t_3\frac{{\text{e}}^{-|r_1-r_2|/\mu }}{|r_1-r_2|/\mu }\mathrm{,}\end{array}

with \mu =1.5 fm and t_3=-6.66 MeV which can give the best preservation of nuclear surface for certain parameter sets, and it also gives the contributions to two-body terms.

The parameter t₁ and t₂ can be determined by fitting the nuclear matter potential,

\begin{array}{}\text{(26)}& V=\alpha \left(\frac{\rho }{{\rho }_0}\right)+\beta {\left(\frac{\rho }{{\rho }_0}\right)}^2\mathrm{.}\end{array}

α term is related to the two-body interaction term, and β term is related to the three-body interaction. The two free parameters α and β can be determined by the requirement of average binding energy (E₀) and compressibility [K_0=9{\rho }_0^2{(\frac{{\partial }^{2}E/A}{\partial {\rho }^2})}_{{\rho }_0}] at the normal density. In order to study the effects of different incompatibilities, one generalizes the potential to be

\begin{array}{}\text{(27)}& V=\alpha \left(\frac{\rho }{{\rho }_0}\right)+\beta {\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }\mathrm{.}\end{array}

The additional third parameter γ can allow us to fix the compressibility independent of other quantities.

The momentum dependent interaction is also important for describing flows over a wide incident energy range with reasonable compressibility. There are the Logarithm-type [203, 204] momentum dependent interaction,

\begin{array}{}\text{(28)}& \begin{array}{l}U(\text{Δ}p)\delta (r-r^{\prime })\\ =1.57{\left[\ln \left(1+5\times {10}^{-4}\text{Δ}{p}^2\right)\right]}^2\rho /{\rho }_0\delta (r-r^{\prime })\mathrm{,}\end{array}\end{array}

with \text{Δ}p=|p_1-p_2| in units of MeV/c and U is in MeV. Another form is the Lorentzian-type momentum-dependent nucleonic mean field V_{md}, which was used in RQMD/S and JAM [205, 206]

\begin{array}{}\text{(29)}& V_{md}={\displaystyle \sum _{k=\mathrm{1,2}}}\frac{C_{ex}^{(k)}}{{\rho }_0}{\displaystyle \int }\text{d}{p}^{\prime }\frac{f(r\mathrm{,}{p}^{\prime })}{1+{\left[(p-p^{\prime })/{\mu }_k\right]}^2}\mathrm{.}\end{array}

C_{ex}^k and {\mu }_k are the parameters of momentum dependent interaction. The mean field potential [Eq. (29)] leads to the following potential energy,

\begin{array}{}\text{(30)}& U_{md}={\displaystyle \sum _{k=\mathrm{1,2}}}\frac{C_{ex}^{(k)}}{{\rho }_0}{\displaystyle \int }\text{d}r\text{d}p\text{d}{p}^{\prime }\frac{f(r\mathrm{,}p)f(r\mathrm{,}{p}^{\prime })}{1+{\left[(p-p^{\prime })/{\mu }_k\right]}^2}\mathrm{.}\end{array}

Exact calculation of integral of Eq. (30) is time consuming in the QMD-type model calculations. Thus, in the actual calculations [206], the momentum dependent potential which used in the relativistic QMD framework [174] is

\begin{array}{}\text{(31)}& {\displaystyle \sum _{k=\mathrm{1,2}}}\frac{C_{ex}^{(k)}}{2{\rho }_0}{\displaystyle \sum _{j(\ne i)}}\frac{1}{1+{\left[{\tilde{p}}_{ij}/{\mu }_k\right]}^2}{\rho }_{ij}\mathrm{,}\end{array}

where {\rho }_{ij}={\displaystyle \int }{\text{d}}^{3}r{\rho }_i(r){\rho }_j(r). In their formulas, the relative distance r_{ij}=r_i-r_j and p_{ij}=p_i-p_j (for convenience, we use r_i and p_i to represent r_{i0} and p_{i0} in the following description) in the potentials were replaced by the squared four-vector distance with a Lorentz scalar,

\begin{array}{}\text{(32)}& \begin{array}{l}{\tilde{r}}_{ij}^2=r_{ij}^2+{\gamma }_{ij}^2{\left(r_{ij}\cdot {\beta }_{ij}\right)}^2\mathrm{,}\\ {\tilde{p}}_{ij}^2=p_{ij}^2-{\left(p_i^0-p_j^0\right)}^2+{\gamma }_{ij}^2{\left(\frac{m_i^2-m_j^2}{p_i^0+p_j^0}\right)}^2\mathrm{.}\end{array}\end{array}

The parameters in Eqs. (28) and (31) are determined by reproducing the real part of the global Dirac optical potential (Schrödinger equivalent potential) of Hama et al. [207], in which angular distribution and polarization quantities in proton–nucleus elastic scatterings are analyzed in the range of 10 MeV to 1 GeV. One should note, α, β, γ should be readjusted by fitting the EOS in the uniform nuclear matter after including the momentum dependent interaction term.

In the collision part, only binary collisions (two-body level) are considered. The collisions are performed in a point-particle sense with a similar way as in cascade [208] without considering the shape of nucleons. Each pair of particles within an same event is tested for a collision at every time step. In details, there is possible collision between particles 1 and 2, if their minimum distance d_{12} in center of mass of colliding pair satisfies

\begin{array}{}\text{(33)}& d_{12}⩽\sqrt{{\sigma }_{tot}(\sqrt{s})/\pi }\mathrm{,}\end{array}

where {\sigma }_{tot}(\sqrt{s}) is the total cross section of incoming the particles 1 and 2 at the center of mass energy \sqrt{s}. {\sigma }_{tot}={\displaystyle {\sum }_i^{N_c}}{\sigma }_i, and i represents the outgoing channel and N_c is the maximum number of outgoing channels. Another condition is that the particles 1 and 2 can move long enough for colliding in the time interval -\delta t/2 to \delta t/2. The channel of outgoing particles (ic) is chosen randomly according to the relative weights of the different reaction cross sections, such as

ξ is a random number. The momenta of the outgoing particles are generated randomly according to the angular differential cross sections and in agreement with the energy–momentum conservation laws. Thus, the cross sections constitute another major part of the model. In the original QMD model, the experimental values of nucleon–nucleon (baryon–baryon) cross sections in free space are used and the medium correction on the cross section is based on it. The Pauli blocking is considered as that in VUU [148, 209-212].

In addition, fragmentation is an important mechanism for intermediate energy heavy ion collisions. Thus, one also needs to identify the fragments at the end of simulations. A reasonable method for identifying the fragments from the simulation results of QMD calculations is needed for obtaining the reaction observables. It is found that the analyzing code also contains important physical contents and influence the final results to a certain extent.

2.3 Improved quantum molecular dynamics model

The original version of QMD code we used was developed in the Frankfurt [213], and there are lots of developments accompanied with appearance of the new generation facilities along the beam energy, as well as isospin degree of freedom. Here, we briefly introduce the improvements we have made in past decades, and it was also known as improved quantum molecular dynamics model (ImQMD).

In our following description, for convenience, we use r_i and p_i to represent the r_{i0} and p_{i0}. The adjustable parameters in the model can be divided into two types. One is related to the numerical calculations, and another is related to physics parameters. The adjustable parameter related to the numerical calculations is the width of Guassian wavepacket, the time step in computation, etc. The adjustable parameters related to the physics are the mean field parameters and in-medium correction parameters, which totally have about 5–15 parameters. The exact number depends on the physics we study.

2.3.1 Nucleonic mean field

With the beam energy decreasing to the Coulomb barrier, the fermion properties of nucleons become more and more important and it naturally requires to develop the original transport model for including or mimicking these effects, such as including the effects from antisymmetric wave function and time evolution of width of wave function, which has led the antisymmetric molecular dynamics model (AMD) [143, 144], Fermionic molecular dynamics model (FMD) [214, 215], extended quantum molecular dynamics model (EQMD) [216]. However, these treatments face a extremely large cost in calculation and thus hardly apply to heavy nuclear systems. One has to find a way to balance the efficiency of computation and physics before the revolution of computation ability.

In the framework of QMD approach, it can also be refined by improving the mean-filed part in two sides. One is to adopt a reasonable energy density functional, which can well reproduce the properties of finite nuclei, in the Hamiltonian equation. For example, one can include the Pauli potential in Hamiltonian [149, 217-220], which is a phenomenological repulsive momentum dependent potential, for mimicking the fermions’ properties. Thus, the equation of motion of centroid of wave packet for particle i reads

\begin{array}{}\text{(35)}& {\dot{r}}_i=\frac{\partial H}{\partial p_i}+\frac{\partial H^{Pau}}{\partial p_i}\mathrm{,}\end{array}

\begin{array}{}\text{(36)}& {\dot{p}}_i=-\frac{\partial H}{\partial r_i}-\frac{\partial H^{Pau}}{\partial r_i}\mathrm{.}\end{array}

Another ad hoc method to mimic this effect is to adopt the phase space constraints [152-154]. This method requires to check whether the phase space occupation number of each particle ({\overline{f}}_i) violates the fermi-dirac distributions, i.e., {\overline{f}}_i>1, during the evolution. If the phase space occupation number is greater than 1, the momentum direction of nucleons will be rearranged to let the phase space occupation number less than or equal to 1. More details can be found in Refs. [152-154].

In the ImQMD model, we mainly improve the mean field part based on the concept of energy density functional. It comes from the approximation used in the QMD approach, where the potential energy as a function of centroid of wave function are used [i.e., \langle {\nabla }_{r_i}V(r_1\mathrm{,}\dots \mathrm{,}{r}_N)\rangle \approx {\nabla }_{\langle r_i\rangle }U(\langle r_1\rangle \mathrm{,}\dots \mathrm{,}\langle r_N\rangle )]. Thus, one can also directly calculate the potential energy U from its energy density functional U={\displaystyle \int }u[\rho ]{\text{d}}^{3}r, by using the \rho (r)={\displaystyle \sum }{\rho }_i=\frac{1}{{\left(2\pi {\sigma }_r^2\right)}^{3/2}}{\text{e}}^{-(r-r_{i0})/(2{\sigma }_r^2)}. In this case, the self contribution to density is considered.

Till now, there are three kinds of Skyrme-type energy density functional used in the ImQMD which depends on the energy region we used.

1) In order to study the heavy ion reaction at low beam energy, we adopt a reasonable energy density functional (EDF) derived from the Skyrme EDF and the Fermi constraints (similar concept as phase space constraints) is used. It is usually named as ImQMD-v2 [155, 162, 221-226], and mainly used at the beam energy above Coulomb barriers and less than Fermi energy. In the simulations, calculations are stopped when the dynamical processes are finished, for example, whether the composite system reach equilibrium. The corresponding stop time is later 1000–5000 fm/c than the projectile and target contacting time.

The Hamiltonian reads

\begin{array}{}\text{(37)}& H=T+U={\displaystyle \sum }\frac{p_i^2}{2m}+{\displaystyle \int }{u}_{\rho }{\text{d}}^{3}r+U_{Coul}\mathrm{.}\end{array}

The potential energy density functional u_{\rho } used in the ImQMD-v2 for low energy heavy ion collisions (near and above the Coulomb barrier energy domain) reads

\begin{array}{}\text{(38)}& \begin{array}{ll}{u}_{\rho }=\hfill & \frac{\alpha }{2}\frac{{\rho }^2}{{\rho }_0}+\frac{\beta }{\gamma +1}\frac{{\rho }^{\gamma +1}}{{\rho }_0^{\gamma }}+\frac{g_{sur}}{2{\rho }_0}{(\nabla \rho )}^2\hfill \\ \hfill & +\frac{C_s}{2{\rho }_0}[{\rho }^2-{\kappa }_s{(\nabla \rho )}^2]{\delta }^2+g_{\rho \tau }\frac{{\rho }^{\eta +1}}{{\rho }_0^{\eta }}\mathrm{,}\hfill \end{array}\end{array}

where the asymmetry \delta =({\rho }_n-{\rho }_p)/({\rho }_n+{\rho }_p), {\rho }_n and {\rho }_p are the neutron and proton densities. The g_{sur} is the coefficient related to the density gradient. C_s is the symmetry potential coefficient, and {\kappa }_s is the parameter related to isospin dependent density gradient term. g_{\rho \tau } term is obtained from the contribution of the \rho \tau term in Skyrme energy density functional, where τ is the kinetic energy density and expressed with the density ρ. The details can be found in Refs. [155, 162, 221, 224, 222]. The parameters in the u_{\rho } as well as the width of the wave-packet

\begin{array}{}\text{(39)}& {\sigma }_r={\sigma }_0+{\sigma }_1{A}^{1/3}\end{array}

in the coordinate space, here A is the number of nucleons in projectile or target, are given in Table 1 which are determined by fitting the properties (including the stability) of nuclei at ground state, the fusion excitation functions of a number of heavy ion fusion reactions at energies around the Coulomb barrier and the charge distributions in multifragmentation process at Fermi energies.

For the lower beam energies, the excitation of system is low and the momentum distribution of the reaction system is not far from the Fermi–Dirac distribution at zero temperature. Thus, we roughly use \eta =5/3 in the g_{\rho \tau } term to approximately describe the contribution from the \rho \tau term. Actually, the exact calculations of \rho \tau terms in transport models should be directly based on the relative momentum of nucleons, and it will give the obviously momentum dependent interaction and effective mass. It will be interesting to check how the momentum dependent interaction influences the low energy reactions.

2) When the beam energies are high enough to trigger the multifragmentation, which has close relation with the EOS in a wide density range, investigation on the nuclear EOS becomes possible, especially for isospin asymmetric nuclear equation of state with the building of new generation rare facility in last couple decades. It requires to develop the transport codes which can incorporate kinds of effective interactions such as the widely used Skyrme interactions (Skyrme potential energy density functionals). At this energy region, the calculations are stopped after 100–200 fm/c when the projectile and target contact, and the exact values depend on the beam energy we studied.

In the version of ImQMD05 [156-158], it is mainly applied in the energy ranging from 20 MeV/nucleon to 300 MeV/nucleon. The Coulomb interaction is as same as the previous treatments, but the nucleonic potential energy density functional, i.e., u=u_{\rho }+u_{md}, reads as

\begin{array}{}\text{(40)}& \begin{array}{ll}{u}_{\rho }=\hfill & \frac{\alpha }{2}\frac{{\rho }^2}{{\rho }_0}+\frac{\beta }{\gamma +1}\frac{{\rho }^{\gamma +1}}{{\rho }_0^{\gamma }}+\frac{g_{sur}}{2{\rho }_0}{(\nabla \rho )}^2\hfill \\ \hfill & +\frac{g_{sur\mathrm{,}iso}}{{\rho }_0}{[\nabla ({\rho }_n-{\rho }_p)]}^2+g_{\rho \tau }\frac{{\rho }^{8/3}}{{\rho }_0^{5/3}}+u_{\rho }^{sym}\mathrm{.}\hfill \end{array}\end{array}

The last term in Eq. (40) is the symmetry potential energy density functional,

\begin{array}{}\text{(41)}& u_{\rho }^{sym}\text{}=\text{}\left[\text{}{A}_{sym}\frac{\rho }{{\rho }_0}\text{}+\text{}{B}_{sym}{\left(\text{}\frac{\rho }{{\rho }_0}\right)}^{\gamma }\text{}+\text{}{C}_{sym}{\left(\text{}\frac{\rho }{{\rho }_0}\right)}^{5/3}\text{}\right]{\delta }^2\rho \mathrm{,}\end{array}

which makes it possible to investigate the different Skyrme forms of density dependence of symmetry energy. For the Skyrme interactions, symmetry potential energy terms come from two-body, three-body and momentum dependent interaction terms, and its related parameters are A_{sym}, B_{sym}, and C_{sym}, respectively. All the parameters in ImQMD, such as α, β, η, g_{sur}, g_{sur\mathrm{,}iso}, g_{\rho \tau } and A_{sym}, B_{sym}, C_{sym} can be derived from the standard Skyrme parameters, x₀, x₁, x₂, x₃, t₀, t₁, t₂, t₃, σ. The relationship between the standard Skyrme parameters and parameters in the ImQMD can be found in Ref. [160].

If one sets A_{sym}=C_{sym}=0 and B_{sym}=\frac{C_s}{2}, u_{\rho }^{sym} becomes

\begin{array}{}\text{(42)}& u_{\rho }^{sym}=\frac{C_s}{2}{\left(\frac{\rho }{{\rho }_0}\right)}^{{\gamma }_i}\rho {\delta }^2\mathrm{,}\end{array}

and one can easily investigate the power law form of density dependence of symmetry energy. In the following text, we use {\gamma }_i to denote the symmetry potential parameters in the power law form of density dependence of symmetry energy.

According to Eq. (41) and Eq. (42), one can write the density dependence of the symmetry energy of cold nuclear matter as

\begin{array}{}\text{(43)}& \begin{array}{ll}S(\rho )\hfill & =\frac{{\hslash }^2}{6m}{\left(\frac{3{\pi }^2\rho }{2}\right)}^{2/3}+A_{sym}\frac{\rho }{{\rho }_0}\hfill \\ \hfill & +B_{sym}{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }+C_{sym}{\left(\frac{\rho }{{\rho }_0}\right)}^{5/3}\mathrm{,}\hfill \end{array}\end{array}

and

\begin{array}{}\text{(44)}& S(\rho )=\frac{{\hslash }^2}{6m}{\left(\frac{3{\pi }^2\rho }{2}\right)}^{2/3}+\frac{C_s}{2}{\left(\frac{\rho }{{\rho }_0}\right)}^{{\gamma }_i}\mathrm{,}\end{array}

respectively. In Eq. (43) and Eq. (44), the first term is kinetic symmetry energy term which comes from the kinetic energy contributions. The rest terms are the symmetry potential energy.

The energy density associated with the mean-field momentum dependence is represented as

\begin{array}{}\text{(45)}& \begin{array}{ll}{u}_{md}=\hfill & \frac{1}{2{\rho }_0}{\displaystyle \sum _{N_1\mathrm{,}{N}_2}}\frac{1}{16{\pi }^6}{\displaystyle \int }{\text{d}}^3{p}_1{\text{d}}^3{p}_2{f}_{N_1}(p_1)f_{N_2}(p_2)\hfill \\ \hfill & \cdot 1.57{\left[\ln \left(1+5\times {10}^{-4}{\left(\text{Δ}p\right)}^2\right)\right]}^2\mathrm{,}\hfill \end{array}\end{array}

f_N are nucleon Wigner functions, \text{Δ}p=|p_1-p_2|. The energy is in MeV and momenta are in MeV/c. With the interaction as in Eq. (45), the Skyrme-EOS will be modified with additional repulsion from momentum dependent interaction. In order to get the reasonable Skyrme EOS under the momentum dependent interaction as in Eq. (45), we use the u_{md}^{\prime } in the model, which is obtained by substraction the u_{md} at T=0 MeV [160], i.e.,

\begin{array}{}\text{(46)}& u_{md}^{\prime }=u_{md}-u_{md}(T=0)\mathrm{,}\end{array}

for refitting the Skyrme EOS.

3) In the version ImQMD-Sky [161, 227], a standard parametrization of the Skyrme potential energy density functional with only the spin–orbit interaction neglected is used. It is also mainly used in the beam energy ranging from 20 MeV/nucleon to 300 MeV/nucleon.

The nucleonic potential energy density functional u is

\begin{array}{}\text{(47)}& \begin{array}{ll}u=\hfill & \frac{\alpha }{2}\frac{{\rho }^2}{{\rho }_0}+\frac{\beta }{\gamma +1}\frac{{\rho }^{\gamma +1}}{{\rho }_0^{\gamma }}+\frac{g_{sur}}{2{\rho }_0}{(\nabla \rho )}^2\hfill \\ \hfill & +\frac{g_{sur\mathrm{,}iso}}{{\rho }_0}{[\nabla ({\rho }_n-{\rho }_p)]}^2\hfill \\ \hfill & +A_{sym}\left(\frac{\rho }{{\rho }_0}\right){\delta }^2\rho +B_{sym}{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }{\delta }^2\rho \hfill \\ \hfill & +u_{md}^{sky}\mathrm{.}\hfill \end{array}\end{array}

The difference between Eqs. (40, 41) and Eq. (47) is that we replace g_{\rho \tau }{(\frac{\rho }{{\rho }_0})}^{5/3}\rho and C_{sym}{(\frac{\rho }{{\rho }_0})}^{5/3}\rho {\delta }^2 in Eq. (40) with the exact Skyrme type momentum dependent interaction terms u_{sky}^{md} as in Eq. (48). The energy density of Skyrme type momentum dependent interaction we used is

\begin{array}{}\text{(48)}& \begin{array}{l}{u}_{md}^{sky}=u_{md}(\rho \tau )+u_{md}({\rho }_n{\tau }_n)+u_{md}({\rho }_p{\tau }_p)\\ =C_0{\displaystyle \int }{\text{d}}^{3}p{\text{d}}^3{p}^{\prime }f(r\mathrm{,}p)f(r\mathrm{,}{p}^{\prime })(p-p^{\prime }{)}^2\\ +D_0{\displaystyle \int }{\text{d}}^{3}p{\text{d}}^3{p}^{\prime }[f_n(r\mathrm{,}p)f_n(r\mathrm{,}{p}^{\prime })(p-p^{\prime }{)}^2\\ +f_p(r\mathrm{,}p)f_p(r\mathrm{,}{p}^{\prime })(p-p^{\prime }{)}^2]\mathrm{.}\end{array}\end{array}

This formula is derived based on the Skyrme-type momentum dependent interaction, \delta (r-r^{\prime })(p-p^{\prime }{)}^2, and f(r\mathrm{,}p) is the nucleon phase space density. In QMD approaches, f(r\mathrm{,}p)={\displaystyle {\sum }_i}\frac{1}{{(\pi \hslash )}^3}\exp [-{(r-r_i)}^2/(2{\sigma }_r^2)-{(p-p_i)}^2/(2{\sigma }_p^2)]. The coefficients C_0 and D_0 can be determined for fitting

\begin{array}{}\text{(49)}& \begin{array}{ll}{\mathcal{H}}_{eff}\hfill & =\frac{1}{8}[t_1(2+x_1)+t_2(2+x_2)]\tau \rho \hfill \\ \hfill & +\frac{1}{8}[t_2(2x_2+1)-t_1(2x_1+1)]({\tau }_n{\rho }_n+{\tau }_p{\rho }_p)\mathrm{,}\hfill \end{array}\end{array}

in nuclear matter at T=0 MeV, where f in Eq. (48) becomes the zero temperature Fermi–Dirac distributions and \tau =\frac{3}{5}{(\frac{3{\pi }^2}{2})}^{2/3}{\rho }^{2/3} in Eq. (49). Finally, we have

\begin{array}{}\text{(50)}& C_0=\frac{1}{16{\hslash }^2}[t_1(2+x_1)+t_2(2+x_2)]\mathrm{,}\end{array}

\begin{array}{}\text{(51)}& D_0=\frac{1}{16{\hslash }^2}[t_2(2x_2+1)-t_1(2x_1+1)]\mathrm{.}\end{array}

Thus, both the effects of symmetry energy potential and the neutron/proton effective mass splitting on the isospin asymmetric heavy ion collisions can be studied simultaneously with this code.

Furthermore, the same parametrization of the potential energy density makes better relation between the studies of heavy ion collisions and the nuclear structure. One should note that if the form of u_{md}^{sky} is adopted in the ImQMD model, one cannot use it to study the reaction at the beam energy above 300 MeV/nucleon. The reason is that the Skyrme type optical potential increases to infinity with the momentum, which violates the optical potential behavior obtained in experiments of nucleon–nucleus reaction.

Correspondingly, the equation of state of cold nuclear matter from the density functional used in the ImQMD can be written as

\begin{array}{}\text{(52)}& \begin{array}{ll}E/A=\hfill & \frac{3{\hslash }^2}{10m}{\left(\frac{3{\pi }^2}{2}\rho \right)}^{2/3}\hfill \\ \hfill & +\frac{\alpha }{2}\frac{\rho }{{\rho }_0}+\frac{\beta }{\gamma +1}\frac{{\rho }^{\gamma }}{{\rho }_0^{\gamma }}+g_{\rho \tau }\frac{{\rho }^{5/3}}{{\rho }_0^{5/3}}\hfill \\ \hfill & +S(\rho ){\delta }^2\mathrm{,}\hfill \end{array}\end{array}

with the gradient terms vanish in the uniform nuclear matter. The density dependence of symmetry energy S(\rho ) is as same as in Eq. (43) or Eq. (44), depending on which form of density dependence of the symmetry energy is used. The terms g_{\rho \tau } in Eq. (52) and C_{sym} in S(\rho ) come from the energy density of Skyrme-type momentum dependent interaction as in Eq. (48) when T=0 MeV in nuclear matter.

2.3.2 Collision part

In transport models for simulating the low-intermediate energy heavy ion collisions, the nucleon–nucleon collisions are determined with the concept of closest approaches. In the ImQMD model, each pair of particles within an same event is tested for a collision at every time step. More explicitly, suppose there is possible collision between particles 1 and 2, specified with (t_0\mathrm{,}{\overrightarrow{r}}_1) and (E_1\mathrm{,}{\overrightarrow{p}}_1), and (t_0\mathrm{,}{\overrightarrow{r}}_2) and (E_2\mathrm{,}{\overrightarrow{p}}_2), respectively, at the current time t₀ in the reference frame of system. In the center-of-mass frame of the two particles, their trajectories {\overrightarrow{\mathcal{R}}}_i^{*}(t^{*}) without a mean field are straight lines pointing along the constant velocities {\overrightarrow{v}}_1^{*}={\overrightarrow{p}}_1^{*}/E_1^{*} and {\overrightarrow{v}}_2^{*}={\overrightarrow{p}}_2^{*}/E_2^{*}. The asterisks represent quantities in the two-particle center-of-mass frame, while quantities without asterisk are in the calculational reference frame. The transformations from system center of mass reference to the two-particle center of mass reference frame for momentum and coordinate are

\begin{array}{}\text{(53)}& p_i^{*}=\left[(\gamma -1)p_i\cdot \frac{\beta }{{\beta }^2}-\gamma E_i\right]\beta +p_i\mathrm{,}\end{array}

\begin{array}{}\text{(54)}& r_i^{*}=(\gamma -1)r_i\cdot \frac{\beta }{{\beta }^2}\beta +r_i\mathrm{,}\end{array}

where,

\begin{array}{}\text{(55)}& \beta =\frac{p_i+p_j}{E_i+E_j}\mathrm{,}\gamma =\frac{1}{\sqrt{1-{\beta }^2}}\mathrm{.}\end{array}

As the trajectories of 1 and 2 are known exactly, the minimum distance can be calculated as

\begin{array}{}\text{(56)}& d_{\perp }^{*2}=({\overrightarrow{r}}_1^{*}-{\overrightarrow{r}}_2^{*}{)}^2-\frac{{[({\overrightarrow{r}}_1^{*}-{\overrightarrow{r}}_2^{*})\cdot {\overrightarrow{v}}_{12}^{*}]}^2}{v_{12}^{*2}}\mathrm{,}\end{array}

with {\overrightarrow{v}}_{12}^{*}={\overrightarrow{v}}_1^{*}-{\overrightarrow{v}}_2^{*}. In the ImQMD codes, the distance condition for a collision to occur is

Here, the frame-independent definition of the cross section (via the impact parameter in the two-particle rest frame) is an important factor in ensuring the approximate reference-frame independence.

Another criteria is to judge whether the collision occurs during the time interval of the current time step. In the ImQMD codes, the time of the closest approach is considered in the two-particle center-of-mass frame, where it may be written as

\begin{array}{}\text{(58)}& t_{\text{coll}}^{*}=t_0^{*}-\frac{({\overrightarrow{r}}_1^{*}-{\overrightarrow{r}}_2^{*})\cdot {\overrightarrow{v}}_{12}^{*}}{v_{12}^{*2}}\mathrm{,}\end{array}

corresponding to the minimum distance d_{\perp }^{*} given by Eq. (56). Note that {\overrightarrow{r}}_1^{*} and {\overrightarrow{r}}_2^{*} are the positions at different times t_1^{*} and t_2^{*}, respectively, in this frame.

In the Bertsch prescription [202], the condition of the closest approach for this time step is set as which is equivalent to t_{\text{coll}}^{*}\in [t_0^{*}-\frac{1}{2}\delta t\mathrm{,}{t}_0^{*}+\frac{1}{2}\delta t]. \delta t is the time interval in c.m. of two colliding particles, and it is related to the \text{Δ}t in computation reference is, \delta t=\alpha \text{Δ}t and α is defined as

\begin{array}{}\text{(59)}& \alpha =\gamma \frac{E_1^{*}{E}_2^{*}}{E_1{E}_2}\end{array}

with Lorentz factor \gamma =1/\sqrt{1-{\beta }^2} where β is the velocity of the center-of-mass of the colliding pair. We have the usual time dilation factor \alpha =1/\gamma in the limit that the two particles have a common velocity. This method can well describe the nucleon–nucleon collision ranging from low-intermediate energy heavy ion collisions. With the beam energy is high enough, the relativistic effects become more and more important, and the time order of nucleon–nucleon collision may depend on the frame. There are some discussions on it [173], but we will not extend this point in detail according to the scope of this paper.

The nucleon–nucleon cross sections and the differential cross sections in free space are taken from Ref. [228]. The in-medium total nucleon–nucleon cross section {\sigma }_{total}^{*} is taken as the form, {\sigma }_{total}^{*}=(1-\eta (E_{beam})){\sigma }_{total}^{free}, where the {\sigma }_{total}^{free} is taken as that in Ref. [228]. {\sigma }_{total}^{*}={\sigma }_{el}^{*}+{\sigma }_{inel}^{*}, and {\sigma }_{el}^{*} and {\sigma }_{inel}^{*} are the in-medium elastic and inelastic two-body cross section. The elastic and inelastic channel is determined by the possibility, P_{el}={\sigma }_{el}^{*}/{\sigma }_{total}^{*}. After generating a random number \xi, the collision for elastic channel is determined by . If the inelastic collision happens, for example NN\to N\text{Δ}, the mass of Δ resonance should be sampled according to a normalized probability distribution for the mass of produced Δ, i.e., P(m\mathrm{,}s) [229], where m is the mass of produced Δ and s is the center of mass energy. The momentum direction of outgoing nucleons is determined by the differential cross section with the Monte-Carlo method by considering the energy–momentum conservation [202].

The decay time of resonance is sampled according to its decay probability in its rest frame, i.e.,

\begin{array}{}\text{(60)}& P(t_{dec})=1-{\text{e}}^{-\frac{t_{dec}}{\tau }}=1-{\text{e}}^{-\text{Γ}{t}_{dec}}\mathrm{,}\end{array}

where τ is 1/\text{Γ}, and Γ is the decay width. Thus, we have

\begin{array}{}\text{(61)}& t_{dec}=-\frac{1}{\text{Γ}}\ln (1-\xi )\mathrm{,}\end{array}

ξ is a random number. The decay time in the computational frame (t_{dec}^{\prime }), is related to that in the rest frame of the decaying particle (t_{dec}) by a Lorentz factor, e.g.,

\begin{array}{}\text{(62)}& t_{A\to BC}^{\prime }=\gamma t_{A\to BC}=\frac{E}{m_A}{t}_{A\to BC}\mathrm{,}\end{array}

where E_A=\sqrt{m_A^2+p^2} is the energy of particle A in computation reference. The momentum of decayed particles can be calculated based on the kinematic relationship.

After the collision or decay, the corresponding momenta will be changed if there is no Pauli blocking, but the positions of them do not change.

2.3.3 Pauli blocking

The outgoing nucleon or baryon after attempted scattering will be checked whether it is Pauli-blocked. The occupation probability f_i^{\prime } at the centroid of the scattered wave packet with final momentum P_i^{\prime } is obtained from the Wigner distribution function corresponding to the QMD wave function given in Eq. (11), with the self-contribution excluded, i.e.,

\begin{array}{}\text{(63)}& \begin{array}{lll}{f}_i^{\prime }\hfill & =\hfill & f_{\tau }({\overrightarrow{R}}_i\mathrm{,}{\overrightarrow{P}}_i^{\prime })\hfill \\ \hfill & =\hfill & \frac{1}{2/(2\pi \hslash {)}^3}\frac{1}{{(\pi \hslash )}^3}{\displaystyle \sum _{k\in \tau (k\ne i)}}{\text{e}}^{-{({\overrightarrow{R}}_i-{\overrightarrow{R}}_k)}^2/(2{\sigma }_r^2)}\hfill \\ \hfill & \hfill & \times {\text{e}}^{-(2{\sigma }_r^2/{\hslash }^2)({\overrightarrow{P}}_i^{\prime }-{\overrightarrow{P}}_k{)}^2}\mathrm{,}\hfill \end{array}\end{array}

with \tau =n or p, which estimates the probability of finding nucleons in a phase space cell of dimension {(2\pi \hslash )}^3. The factor 2/(2\pi \hslash {)}^3 results from consideration of the spin in the phase-space cell. The prefactors combine into a total factor 4. The occupation probability for i is finally taken as \min \left\{\mathrm{1,}{f}_i^{\prime }\right\}.

Furthermore, the criteria \frac{4\pi }{3}{r}_{i^{\prime }k}^3\frac{4\pi }{3}{p}_{i^{\prime }k}^3⩾h^3/8 is also used in the ImQMD codes. i^{\prime } denotes the outgoing nucleon, k represents the other surrounding nucleon which should be in the range of all nucleons except itself. It means that the outgoing nucleon should not be too close to others in phase space. Similar method used in our developed UrQMD model.

2.3.4 Initialization

When solving the equations of motion for nucleons, i.e., Eq. (7), one needs to know the initial coordinate and momentum of each nucleon in projectile and target, which is the staring point of simulating heavy ion collisions. A reasonable initial condition is of vital importance for correctly describing heavy ion collision.

In the initialization, the positions and momenta of the nucleons in reaction system are sampled according to the density and momentum distributions of projectile and target nuclei. Considering the reaction geometric (such as the impact parameters) and the beam energy, the projectile and target nuclei are boosted into the center of mass frame of the reaction system. Since the width of wave packet is fixed and missing some quantum effects, the ground state of nuclei in the QMD approaches may deviate from its ground state in nature. The initial nuclei obtained in QMD approaches usually have certain excitation. In order to obtain the reasonable initial nuclei with less excitation, there are two methods to handle it. One is to reduce the excitation energy by solving the damped equation of motion to find the energy minimum of the system, and it was introduced in Ref. [216]. This method sometimes encounters a difficulty of finding the reasonable damping parameters and the evolution time for different systems. Another method is to select the pre-prepared initial nuclei that satisfactorily describe the properties of projectile and target nuclei, such as the binding energy, the root-mean-square radius ({\langle r^2\rangle }^{1/2}) and the shape for deformed nuclei. Furthermore, the selected pre-prepared initial nuclei should be stable long enough time. One should note that all these check should be under the same nucleon–nucleon interaction or the same energy density functional used in the dynamical calculations.

In the ImQMD model for low energy heavy ion reactions, the second method is adopted [155, 162]. To better describe the properties of neutron-rich nuclei, the neutron skin thickness of neutron-rich nuclei is taken into account in the initialization of the ImQMD model. Based on the 4-parameter nuclear charge radii formula proposed in Ref. [230]

\begin{array}{}\text{(64)}& \begin{array}{ll}{R}_c(\text{f}\text{m})=\hfill & 1.226A^{1/3}+2.86A^{-2/3}\hfill \\ \hfill & -1.09(I-I^2)+0.99\text{Δ}E/A\mathrm{,}\hfill \end{array}\end{array}

with which the 885 measured charge radii can be reproduced with a rms deviation of 0.022 fm (rms deviation is \sqrt{\frac{1}{m}{\displaystyle \sum }{(R_{exp}^i-R_{the}^i)}^2}), and the linear relationship between the neutron skin thickness \text{Δ}{R}_{np}={\langle r_n^2\rangle }^{1/2}-{\langle r_p^2\rangle }^{1/2} and the isospin asymmetry I=(N-Z)/A [231]

\begin{array}{}\text{(65)}& \text{Δ}{R}_{np}=0.9I-\mathrm{0.03,}\end{array}

one can obtain the proton radii R_p=\sqrt{\frac{5}{3}\left[\langle r_c^2\rangle -0.64\right]} from the charge radii {\langle r_c^2\rangle }^{1/2}=\sqrt{\frac{3}{5}}{R}_c and the neutron radii R_n=\sqrt{\frac{5}{3}}\left[{\langle r_p^2\rangle }^{1/2}+\text{Δ}{R}_{np}\right]. The nucleon positions are sampled within the hard sphere with a radius R_p-w_r for the protons and R_n-w_r for the neutrons, respectively. Here, w_r=0.8 fm is to take into account the influence of the width of the wave-packet in the coordinate space. Only those initially prepared nuclei kept good ground state properties and stable for a long enough time are finally selected as the initial nuclei applied for the simulation of the reaction.

For study of fusion reaction at near barrier energy, a more fine procedure is adopted. Figure 1 shows the average density distribution of ²⁰⁸Pb and ¹³²Sn with 500 events. The solid symbols denote the results of the ImQMD model at the initial time, in which the value of {\sigma }_r is obtained based on Eq. (39), and the curves denote the corresponding results of the Skyrme Hartree–Fock calculations with the force SkM* [232]. One can find the event averaged density distribution is in very good consistency with the Skyrme Hartree–Fock calculation at given {\sigma }_r.

In version ImQMD-05 and ImQMD-Sky, the sampling is treated less complicated than in the lower energy version. The binding energy and nuclear radius are adopted to select the reasonable initial nuclei, but the stability of initial nuclei are not forced to keep very long time due to the fact that the reaction time decreases quickly with beam energy increases.

One should note that the initial fluctuation is involved naturally in the QMD-type models. It comes from the randomness of the position and momentum of each particle at the initial time when each event is initialized under the macroscopic conditions. To illustrate the initial fluctuation in QMD type models, Fig. 2 shows the density distribution sampled in the box with periodic boundary condition. The required density profile is \rho (r)={\rho }_0+0.2{\rho }_0\sin (2\pi z/L) with L=20 fm, which is shown as black solid line. The color lines are the sampled density profile in QMD approach from ten different events. One can find that the profiles of the density distribution fluctuate obviously. This initial fluctuation propagates following the evolution of system controlled by the equation of motion Eq. (7) and scattering process. It plays important roles in many processes of heavy ion reactions, for instance, in the MNT reaction, and multiframentation process. However, one should note it is not as exactly same as the fluctuations due to many-body correlations (i.e., f_N(r_1\mathrm{,}\dots \mathrm{,}{r}_N;p_1\mathrm{,}\dots \mathrm{,}{p}_N), with N≥2), which is the main venue to go beyond dissipative mean field dynamics. The physical fluctuation contributes part of the random force in terms of macroscopic dynamics model for heavy ion reactions [126, 127, 140-142] and the other part is from the scattering process [160, 233]. It is still quantitatively unclear what is the strength of the random force, so as the strength of fluctuation, in the heavy ion reactions at different energies. How to extract the initial physical fluctuation obviously needs further study.

2.3.5 Cluster recognization

One obvious character of QMD approach is that it can describe the cluster formation due to the N-body correlations caused by the overlapping wave packets, initial fluctuation, and nucleon–nucleon scattering process fluctuation. They are identified by the cluster recognition method in the real calculations at the end of simulation time at which the dynamical process is finished. All the fragments are still in excitation, and the secondary decay is allowed. Thus, one can expect that the distance between nucleons in the same fragment should not be changed beyond the criteria of nucleon force range.

Based on this idea, the minimum spanning tree algorithm is adopted to recognize the fragments. In this algorithm, the nucleons with relative distances of coordinate and momentum satisfying the conditions and are supposed to belong to a fragment. Here, r_i and p_i are the centroid of a wave packet for ith nucleon in their coordinate and momentum space. R₀ and P₀ are phenomenological parameters determined by fitting the global experimental data, such as the intermediate mass fragments (IMF) multiplicities [150, 156]. They should roughly be in the range of nucleon–nucleon interaction. Typical values of R₀ and P₀ used in the QMD approaches are about 3.5 fm and 250 MeV/c [150, 156], respectively. This approach has been quite successful in explaining some fragmentation observables, such as the charge distributions of the emitted particles, IMF multiplicities [150, 156], yield ratios of free neutrons to protons (n/p), and the double n/p ratios in heavy ion collisions [158].

On the other hand, the MST method fails to describe other details in the production of nucleons and light charged particles [150, 156, 234]. For example, while the yields of Z=1 particles are overestimated, the yields of Z=2 particles are underestimated partly because the strong binding of α particle cannot be well described in transport models. Strong enhancements of the productions of neutron-rich isotopes observed in isoscaling [235], dynamically emitted heavy fragments [236] in neutron-rich HICs, and neutron-rich light charged particles (LCP) at mid-rapidity [237] have not been described well. Furthermore, most transport models predict more transparency than that observed experimentally in central collisions at intermediate energy [160, 238] due to insufficient production of fragments in the mid-rapidity region. Previous studies show that these problems cannot be resolved by changing only the mean field or nucleon–nucleon cross section in transport models.

There have been many attempts to improve the cluster recognition algorithm. More sophisticated algorithms such as the early cluster recognition algorithm (ECRA) [239], the simulated annealing clusterization algorithm (SACA) [240, 241], and the minimum spanning tree procedure with binding energy of fragments (MSTB) [242], have been developed to provide better description of the IMF multiplicities or the average Z_{max} of the fragments. However, these algorithms do not address the lack of isospin dependence in cluster recognition.

By considering the properties of neutron-rich nuclei, such as neutron skin or neutron halo effect and Coulomb effects, in Ref. [159] a method of cluster recognition, namely, the iso-MST was proposed to mimic the isospin dependence in the cluster formation, in which the different values of R_{nn}^0, R_{np}^0, and R_{pp}^0 are chosen. It was found that R_{nn}^0=R_{np}^0=6 fm and R_{pp}^0=3 fm give the suitable description on the isospin related observables in heavy ion collision at intermediate energies.

Our results show that the iso-MST method makes the suppression of Z=1 particles and enhancement of fragments, especially for heavier fragments with Z⩾12. Furthermore, we find enhanced production of neutron-rich isotopes at mid-rapidity. Consequently, isospin-sensitive observables, such as the double ratios, DR(t/³He), and isoscaling parameter α increase to larger values. The widths of the longitudinal and transverse rapidity distributions of Z=1–6 particles also change. In all observables examined, the effects introduced by the iso-MST algorithm are relatively small but in the direction of better agreement with data [159].

2.4 Ultra-relativistic quantum molecular dynamics model

In relativistic heavy ion collisions, kinds of baryons and mesons can be produced and the relativistic effects should be well considered. It is known that the UrQMD model inherits analogous principles as the QMD model [213] in its mean-field part and the relativistic QMD (RQMD) model [174] in the corresponding two-body collision part. It is successfully extended to describe HICs with beam energy starting from as low as several tens of MeV per nucleon (low SIS) up to the highest one available at CERN Large Hadron Collider (LHC).

In the UrQMD model [84, 173, 243, 244], it contains 55 different baryon species (including nucleon, Δ and hyperon resonances with masses up to 2.25 GeV/c²) and 32 different meson species (including strange meson resonances), which are supplemented by their corresponding antiparticle and all isospin-projected states. The collision part of UrQMD has a good treatment on the sequence of collision/decay and on the frame dependence issue by using the collision/decay time table in the computational frame, where the collision time of baryons or decay time of reasonances is calculated based on Eq. (58) and Eqs. (61), (62) and sorted according to the time order. This method is much better for handling the resonance particles’s collision and decay during the propagation time step [229]. When the beam energy is less than 200 GeV/nucleon, the above described algorithm predicts that the particle multiplicities and collision numbers are less than 3% between the laboratory frame and the center of mass frame. Further discussion on the different treatments on the attempted collision are reviewed and compared in Refs. [229, 378].

In addition, based on the cascade mode which is also constantly updated [245], it is also possible to incorporate mean field interactions in the transport calculation. Since the EOS based on the first-principle lattice QCD calculations is still not available in current model investigations, two alternative methods to consider the strong influence of EOS on the dynamics of the expanding system have been tried by our group: the (mean-field) potential updates [73, 74, 246] and the UrQMD+hydrodynamics hybrid mode [247, 248]. In recent years, some improvements we incorporated in the UrQMD especially for HICs at intermediate and low energies are, briefly, as follows: i) including and enriching the potentials such as the symmetry potential [249], the spin–orbit interaction [250], the magnetic field effect [251], the potentials for produced mesons [252, 253], as well as the Skyrme potential energy density functional [255, 254]; ii) incorporating the medium corrections effects on the NN elastic cross sections [77, 172, 256]; iii) improving the Pauli blocking by introducing the restriction of phase space for scattered particles [172].

3 Study of the phenomena and the mechanism in low-intermediate energy heavy ion reactions

3.1 Heavy ion fusion reactions

Searching for the limits of the existence of nuclei is of fundamental importance for nuclear physics study. It is known that demarcation line for existence of a nucleus is driplines beyond which protons or neutrons leak out of nuclei [25-61]. The proton dripline has been reached for many isotopic chains. However, the neutron dripline is known only up to oxygen (Z=8). The superheavy nuclide with Z=118, A=294 marks the current upper limit of nuclear charge and mass. Thus, to search for the way to produce new isotopes near driplines and superheavy nuclei becomes one of the most important tasks in nuclear research. The heavy ion fusion reactions and multi-nucleon transfer reactions are the most important methods for this aim.

In heavy ion fusion reactions, it is of great importance to explore the nucleus–nucleus potential and the formation process (neck dynamics) of compound nucleus, which are usually affected by dynamic effects, nuclear structure effects, isospin effects and so on. In this part, we try to understand it in the framework of microscopic transport model, ImQMD-v2, from the following points: i) dynamical nucleus–nucleus potential, ii) neck dynamics and iii) excitation function of fusion reaction, as they are the most important and relevant aspects relating to the reaction dynamics.

3.1.1 Dynamical nucleus–nucleus potential

The nucleus–nucleus potential is commonly described as a function of the center-to-center distance between the projectile and target nuclei, and consists of a repulsive Coulomb term and a short-ranged attractive nuclear component. Obviously, it evolves with the time during the reaction process, because the shapes of reaction partners evolve with time due to the rearrangement of particles in the system. The nucleus–nucleus potential depends on the reaction energy and the reaction system (for example, the neutron-richness, the strength of shell effects, etc.), which is named as the dynamical nucleus–nucleus potential [257].

Generally, the nucleus–nucleus potential can be calculated based on the nucleon–nucleon interaction of reaction system. By using the microscopic transport model, ImQMD, one can calculate the dynamical nucleus–nucleus potential microscopically [223] in which the densities of the system and the relative distance R between the two nuclei are functions of the evolution time.

When the projectile and target nucleus are well separated (R\gg R_1+R_2) (R₁ and R₂ are the charge radii of the projectile and the target nucleus, respectively), the collective relative motion plays a dominant role and the excitation energy of the reaction partners could be negligible. The nucleus–nucleus potential is thus expressed as

\begin{array}{}\text{(66)}& V_1=E_{c\mathrm{.}m\mathrm{.}}-T_R\mathrm{,}\end{array}

where T_R is the kinetic energy of relative motion of two colliding nuclei, which can be obtained in the ImQMD simulations since the position and momentum of each nucleon can be recorded at every time step in the time evolutions. E_{c\mathrm{.}m\mathrm{.}} is incident kinetic energy associated with the motion of center of mass of the system. It is related to the incident kinetic energy in laboratory E as

\begin{array}{}\text{(67)}& E_{c\mathrm{.}m\mathrm{.}}=\frac{m_B}{m_A+m_B}E\mathrm{,}\end{array}

with A and B are the projectile and target, respectively.

After the di-nuclear system is formed (), the nucleus–nucleus potential is described by a way like the entrance channel potential [258]

\begin{array}{}\text{(68)}& V_2=E_{\text{tot}}(R)-{\overline{E}}_1-{\overline{E}}_2\mathrm{,}\end{array}

where E_{\text{tot}}(R) is the total intrinsic energy of the composite system which is strongly dependent on the dynamical density distribution of the system. {\overline{E}}_1 and {\overline{E}}_2 are the time average of the energies of the projectile and target nuclei, respectively, which are obtained from the energies of the projectile-like and target-like nuclei in the region . R_T=R_1+R_2 denotes the touching point. In the calculations of E_{\text{tot}}(R), {\overline{E}}_1 and {\overline{E}}_2 in Eq. (68), the extended Thomas–Fermi (ETF) approximation for the intrinsic kinetic energy of the reaction system is adopted (see Refs. [223, 259] for details).

The nucleus–nucleus potential is written as a smooth function between V₁ and V₂,

\begin{array}{}\text{(69)}& V_b(R)=\frac{1}{2}\text{erfc}(s)V_2+\left[1-\frac{1}{2}\text{erfc}(s)\right]V_1\mathrm{,}\end{array}

and

\begin{array}{}\text{(70)}& s=\frac{R-R_T+\delta }{\text{Δ}R}\mathrm{,}\end{array}

with \delta =1 fm, \text{Δ}R=2 fm. The obtained nucleus-nucleus potential V_b(R) approaches to V₁ with the increase of R, and approaches to V₂ with the decrease of R.

Figure 3 shows the calculated dynamical nucleus–nucleus potentials for fusion reactions ¹⁶O+⁷⁴Ge and ¹⁶O+¹⁴⁸Nd by using the ImQMD model with the parameter set SkP*. The blue curves denote the corresponding entrance-channel potential with the Skyrme energy–density functional plus the extended Thomas–Fermi approximation including the terms up to second order (ETF2) approach in which the sudden approximation for the densities is used. The empirical barrier distribution functions for these two reactions are presented in Figs. 3(b) and (d). The dashed lines give the positions of the most probable barrier heights B_{\text{m}\mathrm{.}p\mathrm{.}}. The black squares denote the extracted most probable barrier heights from the measured barrier distributions D(E)={\text{d}}^2(E{\sigma }_{\text{fus}})/\text{d}{E}^2 based on the fusion excitation functions. For ¹⁶O+¹⁴⁸Nd, the statistics of the measured data for the fusion cross sections are not many enough to extract the most probable barrier height. It was found that the dynamical barrier height from the microscopic dynamics transport model depends on the incident energy in the fusion reactions [223]. Here, the incident energy E_{c\mathrm{.}m\mathrm{.}}=1.1B_{\text{m}\mathrm{.}\text{p}\mathrm{.}} in the ImQMD simulations. The obtained dynamical barrier height B_{\text{dyn}}\approx B_{\text{m}\mathrm{.}\text{p}\mathrm{.}} at this incident energy for the fusion events. The static potential barriers from the sudden approximation for the densities are evidently higher but relatively thinner than the dynamical ones. To reasonably reproduce the fusion excitation functions, the empirical barrier distributions [see the sub-figures (b) and (d) in Fig. 3] were proposed to take into account the nuclear structure effects and the multi-dimensional character of the realistic barrier in the ETF2 approach, and the value of the peak is lower than the corresponding barrier height from the entrance-channel potential which is based on the spherical symmetric Fermi functions for the densities of the two nuclei and the frozen-density approximation.

To understand the energy dependence of the potential barrier from the view point of the dynamical effects, the density distribution of the fusion events in ¹⁶O+¹⁸⁶W at E_{c\mathrm{.}m\mathrm{.}}=66 MeV was investigated in Refs. [223, 259], and the dynamical deformations of the reaction partners are evident.

In Ref. [257], it was the first time to explore the energy dependence of the potential barrier. It was found that the potential barrier around the Coulomb barrier increases with incident energy increasing, and its up-limits approaching the static one under the sudden approximation. The height of dynamic barrier decreases with decrease of the incident energy, and finally approaches low limit which is close to the adiabatic one.

3.1.2 Neck dynamics

The neck region is defined as connecting the projectile and target nuclei during the reaction, where the density is lowest along the centers of projectile and target nuclei, and it could be neutron-rich warmed, un-compressed region compared with the rest part of the reaction system. It is known that the neck dynamics plays important roles in fusion–fission, light charged particle emission, and so on.

Figure 4 shows the time evolution of the density at neck region in the fusion reaction ¹³²Sn+⁴⁰Ca at an incident energy of E_{c\mathrm{.}m\mathrm{.}}=115 MeV. Here, the density distribution of 328 fusion events over a total of 1000 simulation events is studied for the head-on collisions with the parameter set IQ3a. At t=400 fm/c, the reaction partners begin to touch each other, and the ratio of neutron-to-proton density {\rho }_{\text{n}}/{\rho }_{\text{p}} reaches 2.7 which is higher than N/Z of the compound nucleus by a factor of two. With the increase of the density at neck, the value of {\rho }_{\text{n}}/{\rho }_{\text{p}} deceases with some oscillations and gradually approaches the corresponding neutron-to-proton ratio of the compound nucleus (N/Z=1.37). The extremely neutron-rich density at the neck region can significantly suppress the height of the Coulomb barrier for the fusion reactions induced by neutron-rich nuclei at energies around and below the barrier. Furthermore, one can see from Fig. 4(e) that the surface diffuseness of the reaction partners at the neck side is obviously larger than that at the other side due to the transfer of nucleons.

3.1.3 Fusion cross sections

It is known that the fusion potential between two nuclei is closely related to the surface properties of the nuclei. The neutron skin thickness of neutron-rich nuclei in heavy ion fusion reactions should affect the fusion barrier and thus the fusion cross sections. To explore the influence of the symmetry potential on the fusion excitation function, the fusion cross sections of a number of fusion reactions are systematically investigated with the ImQMD model by adopting different parameter sets. Through creating certain bombarding events (hundreds to thousands) at each incident energy E_{c\mathrm{.}m\mathrm{.}} and each impact parameter b, and counting the number of fusion events, we obtain the fusion probability (or capture probability for reactions leading to super-heavy nuclei) g_{\text{fus}}(E_{c\mathrm{.}m\mathrm{.}}\mathrm{,}b) of the reaction, by which the fusion (capture) cross section can be calculated:

\begin{array}{}\text{(71)}& {\sigma }_{\text{fus}}(E_{c\mathrm{.}m\mathrm{.}})=2\pi {\displaystyle \int }b{g}_{\text{fus}}\text{d}b\simeq 2\pi {\displaystyle \sum }b{g}_{\text{fus}}\text{Δ}b\mathrm{.}\end{array}

To consider the influence of the Coulomb excitation, the initial distance d₀ between the projectile and target should be much larger than the fusion radius. The collective boost to the sampled initial nuclei is given by E_{\text{kin}}=E_{c\mathrm{.}m\mathrm{.}}-Z_1{Z}_2{e}^2/d_0 at the initial time, with the center-of-mass energy E_{c\mathrm{.}m\mathrm{.}}, the charge number Z₁ and Z₂ for the projectile and target nuclei, respectively. Here, the initial distance between the reaction partners at z-direction (beam direction) is taken to be d_0=30 fm for the intermediate reaction systems and 40 fm for the ones with stronger Coulomb repulsion such as ¹³²Sn+⁴⁰Ca. In the calculation of the fusion (capture) probability, event will be counted as a fusion (capture) event if the center-to-center distance between the two nuclei is smaller than the nuclear radius of the compound nuclei (which is much smaller than the fusion radius) and the time evolution will be terminated for those events that will not going to fusion in order to save the CPU time.

As an example, Fig. 5. shows the calculated fusion probability as a function of impact parameters. For the fusion reactions at energies above the Coulomb barrier, the fusion probability looks like a Fermi distribution, i.e., the fusion probability is about one for the central and mid-central collisions. At energies around the Coulomb barrier, the fusion probability decreases quickly with the impact parameter, which implies that the centrifugal potential due to the angular momentum affects the results significantly at this energy region.

A number of fusion reactions were studied with the ImQMD [223, 259]. As examples, the fusion excitation functions for ¹⁶O+⁷⁶Ge and ¹⁶O+¹⁵⁴Sm are shown in Fig. 6. The experimental data can be reproduced very well at energies around the barrier. We also note that in the present version of ImQMD model, the surface diffuseness of heavy nuclei is slightly over-predicted due to the approximate treatment of the Fermionic properties of nuclear system, which causes the over-predicted fusion cross sections at sub-barrier energies for the reactions with heavy target nuclei. In addition, for the fusion reactions with doubly-magic nuclei such as ²⁰⁸Pb, the fusion cross sections at sub-barrier energies are significantly over-predicted by the ImQMD calculations due to the neglecting of the shell effects. The strong shell effect of nuclei can inhibit the dynamical deformation and nucleon transfer, and therefore inhibit the lowering barrier effect. For some neutron-rich fusion systems such as ⁴⁰Ca+⁹⁶Zr and ¹³²Sn+⁴⁰Ca, the results of the ImQMD model are relatively better, which could be due to that the neutron-rich effect is more evident than the shell effect in these reactions.

3.2 Multi-nucleon transfer reactions

The synthesis of nuclei towards driplines and superheavy nuclei, has been of experimental [261-272] and theoretical interest [225, 226, 273-281]. Fusion reactions with stable beams as the traditional approach to successfully synthesize the superheavy elements (SHEs) with Z=107–118 in the last forty years [282-293]. Up to now, the number of observed extremely neutron-rich nuclides is still very limited and about 4000 masses of neutron-rich nuclides in nuclear landscape are unmeasured. For the new nuclei in the “northeast” area of the nuclear map, it is difficult to be reached in the fission reactions and thus fragmentation processes are widely used nowadays. Due to the “curvature” of the stability line, it is also difficult for reaching these new more neutron-rich nuclei in fusion reactions with stable projectiles because of the lack of neutron number.

The revival interest of MNT reaction in intermediate systems or between actinide nuclei at low-energy collisions has arisen [55, 115, 294] for synthesis of the unknown nuclei and superheavy nuclei. As one of the efficient methods to produce neutron-rich nuclei, which are important for understanding nuclear structures at the extreme isospin limit of the nuclear landscape, the MNT reaction may overcome the limits that the number of observed neutron-rich nuclides is very limited at mass region A>160, due to that neither traditional fusion reactions with stable beams nor fission of actinides easily produce new neutron-rich heavy nuclei in this region. The MNT process should compete with fusion process as well as elastic and inelastic scattering for intermediate mass nuclear systems. With the increase of system charge, fusion process is gradually suppressed and completely forbidden for actinide nuclear system due to strong Coulomb repulsion. So, it is interesting to investigate how the MNT reaction mechanism evolves with the size of nuclear system. Here we study the MNT reactions for three reaction systems corresponding to different nuclear size region: i) MNT in ⁸⁶Kr+⁶⁴Ni at 25 MeV/nucleon; ii) ¹⁵⁴Sm+¹⁶⁰Gd at E_{c\mathrm{.}m\mathrm{.}}=440 MeV (=5.6 MeV/nucleon); and iii) ²³⁸U+²³⁸U at 7 MeV/nucleon.

As a N-body and microscopic dynamics model, the ImQMD model has been used to study MNT reactions in heavy ion reactions because a large number of degrees of freedom, such as those in the excitation and deformation of projectile and target, neck formation, nucleon transfer, different types of separation of the composite system, and nucleon emission can be considered simultaneously.

Figure 7 shows the mass vs. total kinetic energy (mass-TKE) distribution in the ImQMD-v2.2 calculations for different impact parameters for the reaction ⁸⁶Kr+⁶⁴Ni at 25 MeV/nucleon [222]. The incident energy is much higher than the Coulomb barrier, which is very suitable for studying the competition among fusion, quasi-elastic scattering, deep inelastic scattering and multifragmentation. As shown in Fig. 7, at central collisions, fusion and binary scattering (quasi-elastic and deep-inelastic scattering leading to multi-nucleon transfer) play equal important roles. And with the increase of impact parameter, the fusion probability approaches to zero gradually, whereas the binary scattering events including the MNT process becomes dominate. Other processes such as ternary breakup and multifragmentation process at this energy also take place. One can see that in central collisions, the neck of the di-nuclear system can be well formed and quickly broadened at such an incident energy and many nucleons are transferred between projectile and target, and the kinetic energy of the relative motion between two colliding nuclei is significantly dissipated to the excitation energy of the composite system. As a consequence, the masses of fragments are distributed in a much broader region, in which the products come from quite different reaction processes, such as fusion, binary breakup, ternary breakup and as well as multifragmentation. With the increase of impact parameter, the mass distribution becomes narrower due to decrease of nucleon transfer between reaction partners.

The production cross sections of isotopes are calculated by using the ImQMD model (ImQMD-v2.2) together with a statistical code [295] for describing the secondary decay of fragments. The mass distributions of elements Z=30 to Z=35 are in good agreement with experimental data as shown in Fig. 8.

The second reaction system used to produce very neutron-rich nuclei by MNT reaction is the neutron-rich reaction system ¹⁵⁴Sm+¹⁶⁰Gd. It was investigated by Wang et al. [296] to study the production of new neutron-rich lanthanides at E_{c\mathrm{.}m\mathrm{.}}=440 MeV. The dynamical study by the ImQMD model (ImQMD-v2.2) shows that it is impossible to produce super-heavy nuclei in this reaction due to the disappearance of the capture pocket and the rapid increase of the potential with decreasing of the relative distance, it might produce new neutron-rich nuclide during the deep inelastic scattering process and moreover the fission barriers of lanthanides are relatively high to prevent fission of heavy fragments in the secondary decay process. Therefore this reaction is interesting for the synthesis of unmeasured lanthanides. By using the same approach used in the reaction ⁸⁶Kr+⁶⁴Ni at 25 MeV/nucleon where the measured isotope distribution of products were reasonably well reproduced [222], the productions of unknown neutron-rich nuclei for reaction ¹⁵⁴Sm+¹⁶⁰Gd at E_{c\mathrm{.}m\mathrm{.}}=440 MeV are calculated. More than 40 extremely neutron-rich unknown nuclei with Z between 58–76 are observed and the production cross sections are at the order of \text{μ}b to mb as shown in Fig. 9. The contour plots are the production probabilities of fragments in logarithmic scale for ¹⁵⁴Sm+¹⁶⁰Gd at the energy of E_{c\mathrm{.}m\mathrm{.}}=440 MeV. The curves denote the β-stability line described by Green’s formula. The circles denote the positions of known masses in AME2012. The results show that multi-nucleon transfer in the neutron-rich reaction system ¹⁵⁴Sm+¹⁶⁰Gd is an efficient way to produce new neutron-rich lanthanides. By analyzing the angular distribution of the produced heavy fragments, it is suggested that might be a suitable angular range to detect these extremely neutron-rich heavy nuclei.

The third reaction system we studied is the ²³⁸U+²³⁸U. For elements with Z>100 synthesized by “cold fusion” with lead and bismuth targets and “hot fusions” with actinide targets, only neutron-deficient isotopes were produced compared with the centers of superheavy elements. At present, the pathway to reach beyond Z=118 is not clear. It is well known that any further experimental extension of the region of SHEs to the center of the predicted first “island of stability” by means of complete fusion is limited by the neutron number of available projectiles and targets and by the very low production cross section as well. The strongly damped reactions between very heavy nuclei, such as U+U (between actinide nuclei) by MNT, could be one possible approach for those purposes [46, 52-54, 109-112, 273].

The theoretical description of the strongly damped reactions between very heavy nuclei, is one of the most difficult problems in nuclear physics due to the extremely large number of degrees of freedom involved. This motivates the use of microscopic dynamical approaches. The ImQMD model incorporating the statistical decay model (HIVAP) for describing the decay of primary fragments was used to study the reaction ²³⁸U+²³⁸U at 7 MeV/nucleon (see references Zhao et al. [226, 278, 279]). In order to calculate the isotope production cross sections of primary and residual fragments for reactions between two actinide nuclei, an empirical model for describing the mass distribution of fission fragments of fissile nuclei was introduced and incorporated with the statistical decay model (details see Ref. [226]).

Figure 10 shows the comparison between the calculated results of the mass distribution of the products in reaction of ²³⁸U around 7 MeV/nucleon and the data [269]. The same scattering angle cut as in the experimental data is selected, that is, only fragments with scattering angles of and in the center-of-mass frame are selected [269, 226]. The open triangles in the figure are the calculation results and the experimental mass spectra from Ref. [269] are indicated by solid squares, open squares, solid circles, open circles, and solid triangles for incident energies of 6.09, 6.49, 6.91, 7.10, and 7.35 MeV/nucleon, respectively. From the figure, we find that the behavior of the calculated mass distribution at 7.0 MeV/nucleon is generally in agreement with the data at the incident energy 7.10 MeV/nucleon, except that the yields at the mass region from 170 to 210 are overestimated compared with the experimental data.

The most important features of mass distribution are considered to be the following: (i) A dominant peak around uranium is observed; this can be attributed to the contribution of the reactions with large impact parameters, (ii) A steep decreasing yield above U with increasing mass number appears. The products at this mass region stem from large mass transfer in small-impact-parameter reactions. (iii) A small shoulder can be seen in the distribution of the products around Pb, compared with the products with a mass near and smaller than uranium for which the yields decrease exponentially as mass decreases. The appearance of the small shoulder is due to the very high fission barrier around Pb. The central and semicentral collisions, and even reactions with b=8–10 fm, contribute to the shoulder in the region around Pb. (iv) In the region below A\approx 190, a double hump distribution is observed. This is clearly due to the fission of actinide and transuranic nuclei, which results in the superposition of symmetric and asymmetric fission.

Figure 11 presents the isotopic production cross sections \sigma (Z\mathrm{,}A) for primary (open symbols) and residual fragments (solid symbols) with charge numbers from Z=94 to 101 in the reaction of ²³⁸U + ²³⁸U at 7.0 MeV/nucleon. An angle cut of 32.5°–44.5° in laboratory system as same as that in the experiment [262, 268] is taken. The experimental data and calculation results from Refs. [262, 268] denoted by black solid stars and blue lines, respectively, are also shown in the figure for comparison. One sees that the experimental data are generally reproduced, except the mendelevium isotopes (Z=101), which was not detected in the experiment yet. Comparing the isotope distributions for primary and residual fragments, we can find the following three features, i.e., the widths of the isotope distributions for residual fragments are much smaller than those for primary fragments for the same element; the peaks in the isotope distributions for residual fragments shift to less neutron-rich side compared with those for primary fragments; the production cross sections for the most probable residual transuranium fragments (Z=96–101) decrease almost exponentially with the increase of fragment charge number. From the isotopic cross section \sigma (Z\mathrm{,}A) for primary fragments and residual fragments, we can calculate the cross section \sigma (Z)={\displaystyle {\sum }_A}\sigma (Z\mathrm{,}A), and the most probable mass number A_Z for the isotope distribution through \sigma (Z\mathrm{,}{A}_Z)=\max [\sigma (Z\mathrm{,}A)] for each element.

Figure 12 shows the cross sections for primary fragments with Z⩾70 which are plotted by black contour lines for the reaction ²³⁸U+²³⁸U at 7.0 MeV/nucleon. It shows that a large amount of primary fragments are produced via proton and neutron transfer between projectile and target. And the most probable isotopes of primary fragments are located near the line with the isospin asymmetry close to that of ²³⁸U (the isospin asymmetry is 0.227) on the nuclear map. The superheavy primary fragment (114, 184) (the isospin asymmetry is 0.235) at the center of the first “island of stability” denoted by cross symbol in red color is not far from this line. The production cross sections for residual fragments are shown by colored rectangles. It can be found that the production cross sections for most of transactinide nuclei are smaller than 10^-8 mb because it is difficult for those primary fragments to survive against fission due to very low fission barrier. For comparison, the area of known nuclei taken from Ref. [298] is presented by the magenta thick line in the figure, and one can find that quite a few unknown neutron-rich isotopes at the “northeast” area of nuclear map can be produced through multinucleon transfer in this reaction.

For the predicted production of light uranium-like elements with , we find that they can reach the border of the proton-rich side of known nuclei in the nuclear map. Because of the high fission barrier, the light uranium-like primary fragments can survive against fission more easily and de-excite through neutron evaporation leading to the production of proton-rich nuclei.

3.3 Low-intermediate and intermediate-high energy heavy ion collisions

With the beam energy increasing, the reaction systems are compressed much more and are highly excited than low energy heavy ion reactions. It leads the system breaking up into much more fragments, three, four, five and more, i.e., from binary fission, ternary fission, to multifragmentation during its expansion phase. The phenomena appeared in this energy region are very rich. One of the characteristic phenomena is the multifragmentation, which was considered as a signal of liquid–gas phase transition. The mechanism of mutifragmentation is still an interesting topic now. Another important phenomenon is the collective flow which reflects the collective motion of emitted nucleons or light charged particles due to the pressure gradient between the participant and the spectator, and thus the collective flow can be used to obtain the stiffness of EOS. Above the threshold energy for Δ production, the resonances appear in the high density region due to the nucleon–nucleon inelastic collisions, and then followed by the production of the mesons. The yields of different kinds of mesons and the ratios between different charged mesons are supposed to be sensitive to the EOS at suprasaturation density. In this section, we will make a review of above phenomena appeared in intermediate energy heavy ion collisions and the calculations based on ImQMD and UrQMD models.

3.3.1 Fragmentation mechansim

As beam energy increases, a ternary breakup was observed experimentally for reaction ¹⁹⁷Au+¹⁹⁷Au at 15 MeV/nucleon [299, 300]. The ternary breakup observed in this reaction has two characters: (i) the masses of three fragments are in comparable size [299, 300] and (ii) the three fragments are nearly aligned along a common reseparation axis [301, 302]. The ternary breakup reaction explored in the experiment is completely different from the commonly known formation of light charged particle accompanied binary fission. The features observed in the ternary breakup reaction between two ¹⁹⁷Au nuclei indicate that the strong dissipation plays an important role in the reaction process, and the deep study of ternary breakup can help us understand the interplay between the one-body or two-body dissipation mechanism.

The ImQMD model is a suitable approach to study the reaction mechanism because in principle, the dissipation from nucleonic potential and nucleon–nucleon collisions, diffusion and correlation effects from classical N-body system are all included in the model without introducing any freely adjusting parameter, and thus, it was applied to study the microscopic mechanism of this new type ternary breakup reaction by Tian et al. [303]. In Ref. [303] the ternary breakup events were selected according to the condition given by experiment (see Refs. [299, 300]), and the microscopic mechanism of the ternary break was studied by the event tracking with the ImQMD-v2 simulation. Figure 13 shows the comparison between the calculation results and experimental data for the mass number distributions of the heaviest fragment A₁, the middle-mass fragment A₂, and the lightest fragment A₃ in the ternary breakup reactions of ¹⁹⁷Au+¹⁹⁷Au at energy of 15 MeV/nucleon. The histograms denote the experimental data of Refs. [299, 300] and the lines with open circles are the calculation results. Clearly, the calculations reproduce the experimental data rather well.

The event tracking in Ref. [303] explored that in the most of ternary breakup events, two ¹⁹⁷Au nuclei had long contact time (about 1000 fm/c or even longer) and a transient composite system was formed before separating into two parts, i.e., the projectile-like fragment (PLF) and target-like fragment (TLF) parts with large nucleon transfer. The snapshots of the time evolution of the ternary breakup process are shown in Fig. 14, which presented the composite system elongating and forming a neck led to system re-separate and then the PLF (or TLF) further quickly separated into two fragments within a time of about 100 fm/c or little short or little longer.

An analysis of the correlation between A₁, A₂, and A₃ was also performed in Ref. [303] in order to further study the mechanism of the ternary breakup process. The distribution of mass asymmetry of A₃ with respect to A₁ and A₂, i.e., {\eta }_1=(A_1-A_3)/(A_1+A_3) and {\eta }_2=(A_2-A_3)/(A_2+A_3), for different impact parameters are given in Fig. 15. It can be found that the probability distributions of {\eta }_1 for central and semi-central collisions (0–6 fm) are of the Wigner-like distributions while those of {\eta }_2 for all impact parameters are of the Poisson-type ones. From the point of view of probability theory, it implies that there exists a certain correlation between fragments A₁ and A₃ and no correlation between A₂ and A₃, which is obviously in consistence with the event tracking analysis discussed previously. From the view point of dynamics, the ternary breaking up in central collisions are in the sequence of (A_1+A_3)+A₂ and then A_1+A_3+A_2. It means the stronger correlations for A₁ and A₃ is kept during the evolution, while the correlation between A₁ and A₃ is weaken as the impact parameter increases.

More crucial information on the mechanism of the studied ternary processes comes from the features of the PLF breakup in its rest frame in the reaction plane. The basic information can be obtained from velocity vectors of the PLF fragments F1 and F2, and TLF in the laboratory reference frame (see Ref. [303]). In Fig. 16, we present the schematic view of the definition of the out-of-plane angle θ, in-plane angle φ, and the angle of the separation axis with respect to the beam direction {\theta }_{c\mathrm{.}m\mathrm{.}}. The reaction plane is defined by the beam direction and the direction of vector of {\overrightarrow{V}}_{PLF}-{\overrightarrow{V}}_{TLF}. This definition is as the same as in Refs. [302, 304].

Figure 17 shows the results of the out-of-plane angle θ and the azimuthal angle φ distributions of fragments from PLF → F1+F2 breakup, and as well as the angle {\theta }_{c\mathrm{.}m\mathrm{.}} distribution obtained from the ImQMD model simulations. The experimental results from Ref. [302] are also shown as red symbols in the figure. Figure 17 clearly tells us that the most of ternary breakup reaction events are in the reaction plane and three fragments are approximately aligned. The study of the ternary breakup reaction of ¹⁹⁷Au+¹⁹⁷Au indicates that the fusion reaction of two very heavy nuclei at energies about 15 MeV/nucleon or higher is forbidden because of the strong Coulomb repulsion but a transient composite system may be formed due to the strong dissipation (strong dumped reaction). The formed transient composite system elongates then breaks up into two parts namely, the PLF and the TLF followed by a further breakup of PLF (or TLF) after a short time, leading to a ternary breakup reaction. Similar mechanism can be extended to quaternary reactions.

Furthermore, one could expect that accompanied light particle emission should have different rapidity distributions for different kinds of breakup mode. Based on this idea, we also analyze the rapidity distributions for light charged particles in binary, ternary and multi- break-up modes.

Figure 18 presents the time evolution of the density contour plots for typical event of reaction ⁶⁴Zn+⁶⁴Zn at E_{\text{beam}}=35 MeV/nucleon which could be the onset energy of multifragmentation [305]. The results are obtained by the ImQMD05 calculations with soft symmetry energy case, i.e., {\gamma }_i=0.5, where {\gamma }_i is the symmetry potential energy coefficients in Eq. (42). As shown in Fig. 18, the system breaks via different break-up modes, such as binary, ternary and multifragmentation break-up modes, i.e., (f1), (f2) and (f3), with different probabilities due to the fluctuations which automatically appears in the QMD type models. It is clear that different break-up modes obviously lead to different emission patterns as well as the different angular and rapidity distributions of LCPs and fragments [for example, Figs. 18(f1–f3)].

To quantitatively understand the correlation between the break-up modes and LCPs emissions, we plot the rapidity distributions of the yields of LCPs corresponding to three kinds of break-up modes, binary (square symbols), ternary (circle symbols) and multifragmentation (triangle symbols) in Fig. 19. Figure 19 presents the rapidity distribution of light charged particle ³He and ⁶He for ⁶⁴Zn+⁶⁴Zn at mid-peripheral collisions. Panel (a) and (b) are for neutron-poor particle, i.e., Y(³He), (c) and (d) are for neutron-rich particles, i.e., Y(⁶He). Left panels are for symmetry potential coefficient {\gamma }_i=2.0 and right panels are for {\gamma }_i=0.5. The rapidity distributions for the yields of ^3.6He are normalized to per event. It is clear that the binary and ternary modes tend to produce more ³He and ⁶He at midrapidity than the multifragmentation breakup mode does. The difference among the yields of neutron-rich light particle Y(⁶He) in the binary, ternary modes and multifragmentation break-up modes is larger compared with that for Y(³He) due to stronger isospin migration. For example, for {\gamma }_i=2.0, the yield of ³He at y_r ≈ 0 from binary mode is 35% larger than that from multifragmentation mode, but the yield of ⁶He at y_r ≈ 0 from binary mode is 70% larger than that from multifragmentation break-up mode. Consequently, one can expect that the different break up modes lead to different shape of rapidity distribution of LCPs.

3.3.2 Fluctuation and chaoticity in liquid–gas phase transition

The theoretical studies on the properties of nuclear matter predict existing a liquid–gas phase transition at temperature of about 10–20 MeV and subnormal densities. A considerable progress has been made on the theoretical study as well as on the experimental side in order to define and collect a converging ensemble of signals connecting multifragmentation to the nuclear liquid–gas phase transition in last couple decades [3, 306-312]. The fluctuation in spinodal region was thought as one of the origins of liquid–gas phase transition in finite nuclear system. An anomalous increase of fluctuation at a phase transition and a rapid increase of chaoticity at the microscopic level stimulate further study of the dynamics of liquid–gas phase transition in nuclear systems.

A way to characterize the dynamics in the phase transition is to calculate the largest Lyapunov exponent (LLE). The largest Lyapunov exponent is defined as [313-315]

\begin{array}{}\text{(72)}& \lambda =\underset{n\to \infty }{{\displaystyle \lim }}\frac{1}{n\tau }\ln \frac{||\text{d}{\overrightarrow{X}}_n||}{||\text{d}{\overrightarrow{X}}_0||}\mathrm{.}\end{array}

The quantity ||\text{d}{\overrightarrow{X}}_n|| is the phase space distance between two trajectories corresponding to two concerned events at time t=n\tau, more details of definition can be found in Ref. [315]. The LLE is a measurement of the sensitivity of the behavior of a system to initial condition and also gives an idea of the velocity at which the system explores the available phase space, and a positive Lyapunov exponent may be taken as the defining signature of chaos. For the case of nuclear fragmentation, a nucleus at a highly excited state eventually breaks into several fragments with nucleons and light particles, where a given trajectory for the system in the phase space will never come back close to the initial state of the system. So a local-in-time LLE over an ensemble of trajectories whose initial conditions are consistent with the nuclei at a given excitation energy should be used. Thus, the time scales of the inverse LLE compared with density fluctuation become more relevant. If the time scale for density fluctuation is much longer than the inverse largest Lyapunov exponent it indicates that the dynamics during fragmentation of the nuclear system is chaotic enough. In this way the LLE calculated over an ensemble of trajectories can carry the full information of dynamics of the systems in multifragmentation.

Figure 20 shows the correlations between the LLE, i.e., λ, and the density fluctuation {\sigma }_{\rho }^2 at different temperatures from 3 to 19 MeV for systems of ¹²⁴Sn of [Fig. 20 (a)] and ²⁰⁸Pb of [Fig. 20(b)]. The density fluctuation {\sigma }_{\rho }^2 is

\begin{array}{}\text{(73)}& {\sigma }_{\rho }^2=\frac{\langle {\rho }^2(t)\rangle -{\langle \rho (t)\rangle }^2}{{\langle \rho (t)\rangle }^2}\mathrm{,}\end{array}

with

\begin{array}{}\text{(74)}& \langle A\rangle ={\displaystyle \int }A\rho (r\mathrm{,}t)\text{d}r\mathrm{,}\end{array}

and the integration is over whole space. One can see that the maximum values of both the LLE and density fluctuation are located at the same temperature, i.e., the critical temperature. There are two branches in \lambda \sim {\sigma }_{\rho }^2, one corresponding to the temperature lower than the critical temperature and another corresponding to the temperature higher than the critical temperature. For the low temperature branch, both the λ and {\sigma }_{\rho }^2 increase as the temperature increases, whereas for the high temperature branch they increase as the temperature decreases. Both branches show that the λ increases roughly linearly with {\sigma }_{\rho }^2. This correspondence between the λ and {\sigma }_{\rho }^2 is qualitatively consistent with the discussion based on Refs. [316, 317].

3.3.3 Collective flow

Collective flow is a motion characterized by space-momentum correlations of dynamic origin, and it carries the information of the pressures generating that motion, the EOS, the in-medium NN cross sections and other properties of the strongly-interacting matter.

As shown in Fig. 21, the flows that have been identified so far are radial, directed, elliptic, triangle flow and high order flow, described by the corresponding coefficients of the components of the Fourier decomposition of the azimuth angle distribution of emitted particles,

\begin{array}{}\text{(75)}& \frac{\text{d}N}{\text{d}\varphi }=p_0(1+2v_1\cos \varphi +2v_2\cos 2\varphi +2v_3\cos 3\varphi +\cdots )\mathrm{,}\end{array}

where v₁ is the directed flow, and v₂ the elliptic flow, v₃ the triangle flow.

The ImQMD05 model were applied to study the directed and elliptic flows and stopping power in ¹⁹⁷Au+¹⁹⁷Au reaction at energies lower than 0.4 GeV/nucleon by using SkP, SkM*, SLy7 and SIII Skyrme interaction parameter sets and compared with the experimental data of INDRA, FOPI and Plastic Wall taken from Ref. [318] and it was found that the SkP and SkM* best fit to the data when the energy dependent in-medium NN elastic cross sections were used [156]. As an example, Fig. 22 shows the excitation function of elliptic flow parameters at midrapidity () for Z≤2 particles for ¹⁹⁷Au+¹⁹⁷Au collisions at b=5 fm [the reduced impact parameter b/b_{max} equals 0.38 and b_{max}=1.15(A_P^{1/3}+A_T^{1/3})]. The general behavior of the excitation functions of elliptic flow parameters v₂ calculated with different Skyrme interactions is similar, i.e., the elliptic flow evolves from a preferential in-plane (rotational like) emission (v_2>0) to out-of-plane (squeeze out) emission with an increase of energies. Clearly, one can see that the harder EOS provides stronger pressure which leads to a stronger out-of-plane emission and thus to a smaller transition energy. The transition energies calculated with SkP and SkM* agree with experimental data [318], while those with SIII and SLy7 are too small compared with experimental data.

For the flow effect at further higher energy heavy ion collisions, it is more suitable to use the UrQMD model. The direct and elliptic flow for protons and light charged particles d, t, and ⁴He at a wide range of beam energies were studied and compared with experimental data [319-322]. In Fig. 23, the excitation function of elliptic flow of protons for reaction ¹⁹⁷Au+¹⁹⁷Au at energies from SIS to RHIC is presented. The UrQMD model in its cascade (UrQMD2.2) and a mean field mode with the HM-EOS and with the DBHF-like medium modification on nucleon-nulceon elastic cross sections are employed. The rapidity cut of has been used because this is the appropriate one to compare with the data at lower energies. Due to this cut, the calculation of the cascade mode (without nuclear potential) reaches negative values at low energies. With inclusion of nuclear potential, the model is in line with experimental data at both SIS and AGS energies. In recent years, for a better description of experimental data at SIS energies, the surface and the surface asymmetry energy terms from the Skyrme potential energy density functional are further incorporated into the UrQMD model [255, 323, 324]. It is found that, with a proper parameter set on the in-medium NN cross section, the published collective flow and nuclear stopping data in HICs at intermediate energies can be reproduced well [172, 252, 253, 255, 256, 323-326].

One can find that the elliptic flow calculated with ImQMD05 and UrQMD are coincidence within the common energy region. For experiments, the important issue is how to determine the reaction plane accurately, since the detectors can not cover full solid angle and can not detect all kinds of particles with the same efficiency. It blocks the experimentalists to get the flow accurately, especially for low-intermediate energy heavy ion collisions. Thus, to further improve the accuracy of the experimental measurement of flow is desirable.

3.4 Spallation reactions

The neutron/proton induced spallation reactions with the beam energy up to 1 GeV have wide applications in material science [327], biology [328], surgical therapy [329], space engineering [330] and cosmography [331]. Interest in the spallation reactions has recently been renewed because of the importance of intense neutron sources for various applications, such as spallation neutron sources for condensed matter and material science [332-334], accelerator-driven subcritical reactors for nuclear waste transmutation [335, 336] or energy production, as well as for medical therapy. So there is a growing need of nuclear data for spallation reactions at intermediate energies up to 1 GeV for targets not only the neutron production materials such as Pb, Bi, W, but also for surrounding structural materials such as Al, Fe, Ni, Zr and biologic elements such as C, O, Ca.

Experimental data are important for designing spallation sources. However, it is impossible to make measurements for all data that are of importance for the various applications [337]. When experimental data are not available, theoretical model calculations have to be employed to estimate the related data. The moving source (MS) model [338], High Energy Transport Code–Three step model (HETC-3STEP) [339], intra-nuclear cascade evaporation (INC/E) model [340] and quantum molecular dynamics (QMD) model [150, 341-344] have often been utilized in the reactions at energies higher than hundreds MeV.

It is well known that spallation reaction is usually described by three-step processes, i.e., the dynamical non-equilibrium reaction process leading to the emission of fast particles and an excited residue, followed by pre-equilibrium emission, and by the decay of the residue. The first process can be described by microscopic transport theory models, the pre-equilibrium is usually optional in different approaches, and the last one can be described by a statistical decay model. Over the past decade, by applying the ImQMD Model (ImQMD05) merged statistical decay model, a series of studies on the proton-induced spallation reactions at intermediate energies have been made [345-349]. Nuclear data including neutron double differential cross sections (DDXs), proton DDXs, DDXs of light complex particles, i.e., d, t, ³He and ⁴He, mass, charge and isotope distributions can be overall reproduced quite well.

As an example, Fig. 24 shows the results of the neutron and proton DDXs for proton-induced spallation reaction for p+²⁰⁸Pb at 256 MeV, 800 MeV, and 1 GeV by using the ImQMD model plus statistical decay model(s). In Fig. 25, the mass and charge distributions of productions in the reaction of 500 MeV p+²⁰⁸Pb [338, 339] and 800 MeV p+¹⁹⁷Au [340, 337], respectively, is shown. In Fig. 26, the isotope distributions of residues produced in the reaction of 750 MeV proton on ⁵⁶Fe is shown. And Fig. 27 shows comparisons between calculated DDXs and experimental data for light charged particles produced in the reaction of 200 and 392 MeV proton on ²⁷Al.

Except the systems shown above, more proton-induced spallation reactions have been analyzed with the ImQMD05 model merging the GEM2 and GEMINI models. The cross sections for products in proton-induced reactions on heavy targets can be reproduced quite well by both models. And it is found that the DDXs of proton and neutron are not sensitive to the parameters of ImQMD and statistical model.

However, to best reproduce the experimental data of mass and charge distributions of productions, two parameters in each model, i.e., the nucleus radius parameter r₀ and the level density parameter modification factor f_f (i.e., a_f^{*}=f_f{a}_f, a_f is the level density parameter) in GEM2; the delay time for fission t_{\text{delay}} and si{g}_{\text{delay}} in GEMINI, need to be tinily readjusted according to the incident energy. For the future, with more precise experimental data becoming available, we expect that the systematic adjustable parameters in GEM2 and GEMINI can be obtained with the present approach. However, there is still some work to be done in order to achieve a universal description for spallation reactions with arbitrary targets and arbitrary incident energy.

Another important point is the description on the LCPs produced in the spallation reaction should be refined. The yields of LCPs with high kinetic energy are underestimated in the previous ImQMD model, as the dashed curves shown in Fig. 27, because the description of the LCPs emission in pre-equilibrium process is absent. In order to overcome the limitation of the ImQMD05 model in description of LCPs emission, a phenomenological surface coalescence mechanism is introduced into the ImQMD model. The basic idea of this mechanism is: the leading nucleon ready to leave from compound nuclei can coalesce with other nucleon(s) to form a LCP, and the LCP with enough kinetic energy to overcome Coulomb barrier can be emitted. By systematic comparison between calculation results and experimental data of nucleon-induced reactions, the parameters in the surface coalescence model are fixed. Then with the fixed parameters, chosen once for all, the prediction power of the model is tested by the nucleon-induced reactions on various targets with energies from 62 to 1200 MeV. And it is found that, with surface coalescence mechanism introduced into ImQMD model, the description on the DDXs of LCPs is greatly improved. And Fig. 27 shows comparisons between calculated DDXs and experimental data for light charged particles produced in the reaction of 200 and 392 MeV proton on ²⁷Al.

As we mentioned before, the advantage of the ImQMD model is that the model parameters are directly obtained from the Skyrme interactions. Thus, one can study the impact of the different effective Skyrme interaction on proton-induced spallation reactions. Based on the best Skyrme interaction we obtained, one can also obtain the information of EOS in the mean field level.

4 Study of in-medium NN cross sections

As we discussed in previous section, collision term I_{coll} is another key ingredient of transport equation, and the in-medium effects on the collision cross sections influences the results of the simulation of HICs. Thus, the study of in-medium cross sections and extracting the medium corrections to the cross sections from experimental observables are important, which are the main contents of this section.

4.1 In-medium NN cross sections from microscopic approach

The in-medium NN cross sections can be calculated by using Bruckner G-matrix method [357-363], one-boson-exchange model (OBEM) [364-366] or the closed-time path Green’s function approach [367-370] in theory. By using the closed-time path Green’s function method, the in-medium elastic and inelastic nucleon-nucleon cross sections can be derived within the framework of self-consistent relativistic Boltzmann–Uehling–Uhlenbeck (RBUU) equations. In the early works [367-369], the in-medium elastic and in-elastic cross sections without considering the isospin dependence were computed based on quantum hadro-dynamical (QHD-1) model. Up to now, by applying the effective Lagrangian such as the QHD-II and its extension in which the couplings to vector–isovector {\rho }^{\mu } and scalor-isovector \delta (a_0(980)) mesons are considered, the isospin asymmetry α, density ρ, and center-of-mass energy \sqrt{s} dependence of the in-medium cross sections {\sigma }^{*}, such as {\sigma }_{NN\to NN}^{*}(\sqrt{s}\mathrm{,}\rho \mathrm{,}\alpha ) and {\sigma }_{NN\to N\text{Δ}}^{*}(\sqrt{s}\mathrm{,}\rho \mathrm{,}\alpha ) are computed [371].

In Fig. 28 we present both the cross sections and the suppression parameters as a function of c.m. energy \sqrt{s} at different densities, which are obtained with closed-time path Green’s function as in Ref. [369]. {\sigma }_{np}^{*(2)} and {\sigma }_{nn\mathrm{,}pp}^{*(2)} are in panels (a) and (b), and the suppression parameters, {\eta }_{np}^{(2)}={\sigma }_{np}^{*(2)}/{\sigma }_{np}^{free(2)} and {\eta }_{nn(pp)}^{(2)}={\sigma }_{nn(pp)}^{*(2)}/{\sigma }_{nn(pp)}^{free(2)} are in panels (c) and (d). The superscript “(2)” means that the cross sections are obtained with closed-time path Green’s function method as described in Ref. [157]. One can see that {\sigma }_{np}^{*(2)} changes little with density and is nearly the same as in free space. On the other hand, the cross section {\sigma }_{nn\mathrm{,}pp}^{*(2)} tends to be suppressed at lower energies, \sqrt{s}⩽2.05 GeV and enhanced at higher energies. Differences in the features of the two cross sections are associated with the differences between the isospin T=0 and the T=1 channels and, in particular, presence of a low-energy resonance in the T=1 channel and effects of {\rho }^{\mu } exchange.

Figure 29 shows the in-medium cross section {\sigma }_{NN\to N\text{Δ}}^{*} (the initial isospin averaged one) as a function of \sqrt{s} for different reduced densities and the calculation was performed with the Δ mass taken to be the pole mass [371]. It is found that at high energies, the {\sigma }_{NN\to N\text{Δ}}^{*} monotonously decreases with increasing density. This is similar to previous calculations [368] performed with another parameter set of the equation of state. However, when approaching the threshold energy, the density dependence is somewhat different from that in Ref. [368] due to the neglect of the Δ mass distribution. This should play a more important role at the lower energies. One might argue that the mass distribution ought to be influenced by the nuclear medium as well, a topic that deserves further investigation.

The results become more interesting for the isospin asymmetries \alpha =({\rho }_n-{\rho }_p)/({\rho }_n+{\rho }_p)\ne 0 to be discussed in the following. Owing to the difference in isospin matrices shown in Ref. [371], the cross-section values for the production of {\text{Δ}}^{++} and {\text{Δ}}^{-} are exactly three times larger than those of other channels ({\text{Δ}}^{+} and {\text{Δ}}^0) as long as isospin asymmetry \alpha =0. With \alpha \ne 0, this relation is destroyed because of the different effective masses for {\text{Δ}}^{++}, {\text{Δ}}^{+}, {\text{Δ}}^0, {\text{Δ}}^{-} and neutron, proton, which make the scattering amplify different as well. It is found that the influence of the mass splitting on {\sigma }_{pp\to n{\text{Δ}}^{++}}^{*} and {\sigma }_{nn\to p{\text{Δ}}^{-}}^{*} are much stronger than that on the other four channels. Figure 30 shows the effect of mass splitting in the isospin asymmetric medium on the in-medium cross sections of Δ production. The ratios R(\alpha )={\sigma }^{*}(\alpha )/{\sigma }^{*}(\alpha =0) of all channels are shown as an example in Fig. 30 as a function of the isospin asymmetry at a typical kinetic energy E_K = 1 GeV (correspondingly, s^{1/2}=2.326 GeV) where the maximum of the cross section approaches. The ratio R(\alpha ) deviates almost linearly from unity when the value of the isospin asymmetry increases from 0 to 0.5. It occurs with the sequence: R(\alpha \mathrm{,}pp\to n{\text{Δ}}^{++})>R(\alpha \mathrm{,}pp\to p{\text{Δ}}^{+})>R(\alpha \mathrm{,}pn\to n{\text{Δ}}^{+})>R(\alpha \mathrm{,}pn\to p{\text{Δ}}^0)>R(\alpha \mathrm{,}nn\to n{\text{Δ}}^0)>R(\alpha \mathrm{,}nn\to p{\text{Δ}}^{-}). It is further seen that, at \alpha =0.5, the R(\alpha ) ratio remains within the interval between 0.88 and 1.15 for the {\text{Δ}}^{+} and {\text{Δ}}^0 production channels, while it changes more rapidly to 1.74 and 0.73 for {\text{Δ}}^{++} and {\text{Δ}}^{-}, respectively.

Cui et al. also studied the influence of threshold effects on the in-medium NN\to N\text{Δ} cross sections [372] by using the one-boson-exchange model (OBEM), where the NN\to N\text{Δ} cross section was averaged over the mass distribution of Δ considering the Δ as short-living resonance. They found that the isospin splitting of medium correction factor R={\sigma }_{NN\to N\text{Δ}}^{*}/{\sigma }_{NN\to N\text{Δ}}^{free} are weaken when one consider the threshold energy effects, i.e., the changes of self-energy when a nucleon turns into a Δ in process of NN\to N\text{Δ}. In Fig. 31, the R as a function of energy above the threshold, i.e., Q=\sqrt{s}-{\sqrt{s}}_{th}=\sqrt{s}-(2m_N+m_{\pi }), for different beam energies are presented. By considering the threshold effects in the calculation of NN\to N\text{Δ} cross sections, the isospin splitting of R becomes weaker at the beam energy above 0.8 GeV because the changes of scalar and vector self-energies in the process of NN\to N\text{Δ} become smaller relative to the kinetic energy part.

Constraining those cross sections is essential for reducing the uncertainties of EOS constraints by using the transport models [157, 160, 373]. Also the in-medium cross sections are of interest for their own sake, as they underly the viscosity and other nuclear transport coefficients [374]. One of the powerful HICs observables for extracting the in-medium NN cross sections is the stopping power [375, 376] because it is constructed from the changing of momentum of the emitted particles and reflects the cross sections in the most direct manner. Furthermore, their correlation [375] with flow observables which are widely used for determining the stiffness of EOS, is also of interest. We employed the ImQMD and UrQMD models to investigate the effects of in-medium NN cross sections on HICs observables. The cross checking on the results with two models is helpful for us to deeply understand the in-medium NN cross sections and the reliability of transport models, and it will be discussed in the next part.

4.2 In-medium NN cross sections from heavy ion collisions

It is quite often in the transport model calculations that the medium corrected elastic cross section adopts a form of {\sigma }_{NN}^{*}=(1-\eta \rho /{\rho }_0){\sigma }_{NN}^{free}, with η being a parameter and usually set to be 0.2. Compared with that obtained from the theoretical calculations shown in Fig. 28 where the {\sigma }_{NN}^{*} depends also on energy in addition to density, the form of the medium corrections {\sigma }_{NN}^{*}=(1-\eta \rho /{\rho }_0){\sigma }_{NN}^{free} is too simple. To obtain the information of the in-medium elastic cross sections from HIC, Zhang et al. [157] proposed ad hoc parametrization inspired by the closed-time path Green’s function (CTPGF) results aiming at the description of the excitation function for elliptic flow in the mid-rapidity region of in Au+Au collisions. Figure 32 shows the ImQMD05 calculation results of excitation functions of elliptic flow parameter v₂ and directed flow P_{xdir}^0=p_{xdir}/u_{c\mathrm{.}m\mathrm{.}}^{beam} with different Skyrme density functionals and in-medium elastic cross sections. There are three kinds of in-medium NN cross sections, i.e., {\sigma }^{*(1)}, {\sigma }^{*(2)}, and {\sigma }^{*(3)}, where {\sigma }^{*(1)} is the one given in expression [{\sigma }_{NN}^{*}=(1-\eta \rho /{\rho }_0){\sigma }_{NN}^{free} and \eta =0.2], {\sigma }^{*(2)} is the results calculated by using CTPGF method based on QMD-II effective Lagrangian, {\sigma }^{*(3)} is the ad hoc parameterization given by Ref. [157]. It is seen that {\sigma }^{*(3)} shows the best fit to the data and {\sigma }^{*(1)} is the worst among three kinds in-medium cross section. It is found that the best fit to the experimental flow data with η=0.2 for E_{beam}⩽ 150 MeV/nucleon, η=0 for E_{beam}=150–200 MeV/nucleon, \eta =-0.2 for E_{beam} = 200–400 MeV/nucleon, can well reproduce the data of flow and stopping power simultaneously.

The impact of nucleon-nucleon cross sections on nuclear stopping is also investigated. Figure 33 shows a variety of results pertinent to vartl. The vartl is defined as the ratio of the rapidity variance in the transverse direction to the rapidity variance in the longitudinal direction,

\begin{array}{}\text{(76)}& vartl=\frac{\langle y_t^2\rangle }{\langle y_z^2\rangle }\mathrm{,}\end{array}

which has been used as a measure of the nuclear stopping. The emitted charged particles with Z=1–6 from central collisions of symmetric or near-symmetric systems are selected, with the contributions of different particles weighted with Z. The panels (a) and (b) in Fig. 33 show calculated distributions in longitudinal and transverse rapidities in 400 MeV/nucleon Au+Au collisions at . The corresponding vartl values for different calculations are quoted in those panels. The panel (c) compares the calculated vartl excitation functions for Au+Au to the data. Finally, the panel (d) compares the calculated dependence of vartl on system charge to data at 400 MeV/nucleon. From Fig. 33, one can learn that there is a good chance to restrict the in-medium cross sections using measured vartl but less chance to restrict the mean field. Our results show both flow and stopping power data favor the {\sigma }^{*(3)}, which is consistent with the recent analysis from INDRA data [377].

Very recently, Li et al. also discussed the similar things with UrQMD model [256]. They adopted a more complicated medium correction form, which reads

\begin{array}{}\text{(77)}& {\sigma }_{NN}^{*}=\mathcal{F}(\rho \mathrm{,}p){\sigma }_{NN}^{free}\mathrm{,}\end{array}

with

\begin{array}{}\text{(78)}& \mathcal{F}(\rho \mathrm{,}p)=\frac{\lambda +(1-\lambda )e^{-\rho /{\rho }_0\xi }-f_0}{1+{(p/p_0)}^{\kappa }}+f_0\mathrm{.}\end{array}

p is the momentum of nucleon in center of mass of two colliding nucleons, ξ and λ are the parameters which determine the density dependence of the cross sections, κ is the parameter which determine the momentum dependence of the cross sections.

As shown in Fig. 34, the degree of nuclear stopping \left(R_E=\frac{E_{\perp }}{2E_{/\text{}\text{}/}}\right) in central Au + Au collisions as a function of the beam energy are investigated. The results obtained with medium correction factor \mathcal{F}(\rho \mathrm{,}p)={\sigma }_{NN}^{*}/{\sigma }_{NN}^{free}=0.5 are the largest and those with \mathcal{F}=0.2 are the smallest of all. Again, this result from the nuclear stopping observables also consistently supports that the medium correction factors of about 0.2 and 0.5 are required for reasonably describing the degree of nuclear stopping in HICs at E_{lab}=40 and 150 MeV/nucleon, respectively. Furthermore, the difference between the results from the FU3FP1 [corresponding to \lambda =1/3, \xi =1/3, f_0=1, p_0=0.425 GeV/c, and \kappa =5 in Eq. (78)] parametrization and from \mathcal{F}=0.3 steadily increases with increasing beam energy. Calculations with FU3FP1 fit the experimental data quite well and reproduce the slightly increased stopping power with increasing beam energy, whereas others fail to reproduce the observed beam-energy dependence. In the future, rigorous extracted information of in-medium NN cross sections also require a refined Pauli blocking in the transport models [378].

4.3 In-medium NN cross sections from spallation reactions

Proton induced reactions also provide the information of in-medium NN cross sections below the normal density with a much clearer picture. In the following, the medium correction is taken as

\begin{array}{}\text{(79)}& {\sigma }_{tot}^{\text{*}}={\sigma }_{in}^{free}+{\sigma }_{el}^{*}={\sigma }_{in}^{free}+\mathcal{F}(u\mathrm{,}\delta \mathrm{,}p){\sigma }_{el}^{free}\mathrm{,}\end{array}

where {\sigma }_{el}^{free} and {\sigma }_{in}^{free} are the free isospin dependent elastic and inelastic cross sections, respectively. The form of the medium correction factor is as same as that proposed by Li et al. in Refs. [77, 172, 255]. \mathcal{F}(u\mathrm{,}\delta \mathrm{,}p) depends on the nuclear-reduced density u=\rho /{\rho }_0, the isospin-asymmetry \delta =({\rho }_n-{\rho }_p)/({\rho }_n+{\rho }_p) and the momentum p_{NN}. On condition of momentum independent medium correction, \mathcal{F}(u\mathrm{,}\delta \mathrm{,}p)=F_{\delta }\cdot F_u,

\begin{array}{}\text{(80)}& F_u=\lambda +(1-\lambda )\exp (-u/\zeta )\mathrm{,}\end{array}

\begin{array}{}\text{(81)}& \begin{array}{l}{F}_{\delta }=1-{\tau }_{ij}A(u)\delta \mathrm{,}A(u)=\frac{0.85}{1+3.25u}\mathrm{,}\\ i=j=n/p\mathrm{,}{\tau }_{ij}=\mp 1;i\ne j\mathrm{,}{\tau }_{ij}=0.\end{array}\end{array}

If the medium correction has an obvious momentum dependence, \mathcal{F}(u\mathrm{,}\delta \mathrm{,}p)=F_{\delta }^p\cdot F_u^p,

\begin{array}{}\text{(82)}& F_{\delta /u}^p=\{\begin{array}{ll}{f}_0\mathrm{,}\hfill & p>1\text{G}\text{e}\text{V}/c\mathrm{,}\hfill \\ \frac{F_{\delta /u}-f_0}{1+{(p/p_0)}^{\kappa }}+f_0\mathrm{,}\hfill & p⩽1\text{G}\text{e}\text{V}/c\mathrm{,}\hfill \end{array}\end{array}

with p being the momentum in the NN center-of-mass two colliding nucleons. By varying the parameters λ, ζ, f₀, p₀ and κ one can obtain different medium correction on nucleon–nucleon elastic cross section (NNECS). The parameter sets used are listed in Tables 2 and 3. Among these parameter sets, FU1–3 and FP1–5 are taken from Refs. [77, 172, 255], FU4 and FP6 are used in Ref. [379].

The various in-medium NNECS obtained from combinations by parameterizations FU1, FU2, FU3 and FP1, FP2, FP3 are tested by the excitation function of reaction cross section (RCS) for p+⁵⁶Fe in Fig. 35. By using the in-medium NNECS, the descriptions on the excitation function of RCS are great improved. Especially the combinations of FU2+FP3 and FU3+FP3 give the best two results except a little deviation at low energies. According to the momentum dependence of the in-medium NNECS F_u^p, the enhancement effect of momentum correction does not yet appear at such low energies if FP3 is adopted. FU4 provides a reasonable correction effect between ones given by FU2 and FU3.

More experimental data of the RCS of proton induced reaction from light to heavy targets ¹²C, ²⁷Al, ^40,48Ca, ⁹⁰Zr, ¹¹⁸Zn, ²⁰⁸Pb are used to test the obtained in-medium NNECS FU4FP6. As examples, the excitation functions of RCS for nucleon-induced on various targets, such as n+¹²C, n+⁶³Cu, n+²⁰⁸Pb, p+¹²C, p+⁴⁰Ca, and p+¹¹⁸Sn, calculated by the ImQMD model with FU4FP6 are shown in Fig. 36. One can see that all experimental data are reproduced quite well. In particular, not only the excitation functions of RCS for the targets along the β-stable line are well described, but also the excitation function of RCS for ⁴⁸Ca, which is far away from the β-stable line, is also well reproduced. In addition to the proton-induced reactions, the good description on the excitation functions of RCS for the neutron-induced reactions means that the medium correction of the isospin dependence of NNECS is reasonably reproduced.

5 Symmetry energy at different densities and temperatures and the constraints

The relationship between pressure, density, temperature and isospin asymmetry [\delta =({\rho }_n-{\rho }_p)/({\rho }_n+{\rho }_p)] described by the equation of state (EOS) of nuclear matter governs the compression achieved in supernova and neutron stars, as well as the structure and many other basic properties of nuclei. Thus, determination of the EOS has been one of the primary goals since first relativistic heavy ion beams started to operate in the beginning of the 80s of the last century [5-8]. The symmetry energy of nuclear matter is defined as the difference between the energy per nucleon in pure neutron matter and in symmetric nuclear matter, i.e., the energy related to the isospin asymmetry of the system. As one of the most important and unclear part of the EOS for isopin asymmetric nuclear matter, it plays an important role in the various fields, ranging from the structure of nuclei to gravitational collapse to neutron stars.

Within the parabolic approximation, the energy per nucleon for asymmetric nuclear matter E(\rho \mathrm{,}\delta ) (known as EOS of cold nuclear matter) can be written as

\begin{array}{}\text{(83)}& \begin{array}{ll}E(\rho \mathrm{,}\delta )\hfill & =E(\rho \mathrm{,}\delta =0)+E_{sym}\hfill \\ \hfill & =E(\rho \mathrm{,}\delta =0)+S(\rho ){\delta }^2+\cdots \mathrm{.}\hfill \end{array}\end{array}

The first term is the energy per nucleon for symmetric nuclear matter, E_{sym} is the symmetry energy term. Here, one should note that the S(\rho ) in the second term describes the density dependence of the symmetry energy which was also named as symmetry energy in other literatures [62]. The approaches to describe the nuclear EOS have been performed by using the variety of effective interactions within mean-field theories, relativistic [380-389], non-relativistic [390-396] as well as the effective interaction based on chiral effective field theory [397-404]. In addition, the ab initio approaches based on high precision free space nucleon–nucleon interactions and the nuclear many-body problems being treated microscopically are also applied to study the nuclear EOS [405-421]. However, large uncertainties in predictions of the density dependence of the symmetry energy away from normal nuclear matter have been found. It stimulates lots of effort to reduce the uncertainties of the predicted density and momentum dependence of symmetry potential by using very neutron-rich HICs.

Quite a few important observables in HICs, such as the yield ratio of neutron to proton, the yield ratio of {\pi }^{-} to {\pi }^{+}, and elliptic flow, were predicted to be sensitive observables for sub- and supra-saturation densities based on the IBUU approach [63, 66, 69-71, 422]. Li et al. also pointed out that comprehensive studies with multiple observables such as the yield ratios of free neutrons to protons, {\pi }^{-} to {\pi }^{+}, ³H to ³He, can provide the density dependence of the symmetry energy at a wider density region based on the UrQMD model calculations as shown in Refs. [73-75]. Tsang et al. proposed and measured the isospin sensitive observables, double neutron to proton ratios [423], isospin diffusion [424], isospin transport ratios as a function of rapidity [425]. At present, the data for the isospin diffusion [424] and the neutron to proton yield ratio data [423] from NSCL/MSU, the mass asymmetry for largest and second largest residues in CHIMERA and the flow data from TAMU [426] and GSI [236], angular distribution of neutron to proton exceed ratios from Tsinghua/IMP [427], are important observables for constraining the symmetry energy around saturation density. The data of pion yield ratio, {\pi }^{-}/{\pi }^{+} [78] and elliptic flow of neutron and proton [249, 428, 429] from GSI are available, which are supposed to be used to constrain the symmetry energy at supra-saturated densities. By comparing the data with transport model calculations, one indirectly obtains the density dependence of the symmetry energy. With the development of new generation of neutron-rich beam facilities in Institute of Modern Physics (IMP), National Superconducting Cyclotron Laboratory (NSCL), Texas A&M (TAMU), National Institute for Nuclear Physics (INFN), GSI, tight constraints of the density dependence of the symmetry energy will be possible in future.

5.1 Symmetry energy at subsaturation density

The discussions in this section focus on the neutron to proton yield ratios, isospin diffusion as observables to constrain the density dependence of the symmetry energy from normal density to subnormal density. The ImQMD code is used to extract the information of the symmetry energy from heavy ion collision data. Here, one should note that the exact input of the transport model is the effective nucleon interaction, potential energy density functional or the single particle potential as mentioned in Section 2. Consequently, one can also obtain the EOS or symmetry energy at different density or temperature based on the interaction or nucleonic potential parameters that we used in the transport models. In the community, the EOS or symmetry energy at zero temperature, which are obtained with the interaction used in the transport models, are usually used to describe the constrained results on EOS or the density dependence of the symmetry energy S(\rho ).

5.1.1 Constraints from n/p ratios and isospin diffusion

The single neutron to proton yield ratio was proposed to study the density dependence of the symmetry energy by Li et al. in Ref. [63] in 1997 where the neutron proton yield ratios are related to the strength of symmetry energy. The first experimental data was published in 2006 for ^112,124Sn+^112,124Sn by the NSCL/MSU group [423]. However, the different efficiencies for neutron and charged particles detectors cause the large errors for single n/p ratio. A double n/p ratio was proposed in order to reduce the uncertainties in the neutron detection efficiencies and the energy calibrations of neutrons and protons. The double ratio

\begin{array}{}\text{(84)}& \begin{array}{l}DR(Y(n)/Y(p))=\frac{R_{n/p}(A)}{R_{n/p}(B)}\\ =\frac{\text{d}{M}_n(A)/\text{d}{E}_{c\mathrm{.}m\mathrm{.}}}{\text{d}{M}_p(A)/\text{d}{E}_{c\mathrm{.}m\mathrm{.}}}\frac{\text{d}{M}_p(B)/\text{d}{E}_{c\mathrm{.}m\mathrm{.}}}{\text{d}{M}_n(B)/\text{d}{E}_{c\mathrm{.}m\mathrm{.}}}\mathrm{,}\end{array}\end{array}

is constructed from the ratios of energy spectra, \text{d}M/\text{d}{E}_{c\mathrm{.}m\mathrm{.}} of neutrons and protons for two systems A and B characterized by different isospin asymmetries. The comparison between the ImQMD05 calculations and the experimental data is shown in Fig. 37. The star symbols in Fig. 37(a) show the neutron–proton double ratios measured at as a function of center-of-mass (c.m.) energy of nucleons emitted from the central collisions of ¹²⁴Sn+¹²⁴Sn and ¹¹²Sn +¹¹²Sn [423]. Despite the large experimental uncertainties for higher energy data from Ref. [423], the comparisons between the data and ImQMD calculations definitely rule out both very soft ({\gamma }_i=0.35) and very stiff ({\gamma }_i=2) density-dependent symmetry terms. Here, {\gamma }_i is the symmetry potential parameter as in Eq. (42) in Section 2.

The parameterization of the symmetry energy corresponding to the interactions we used in these calculations has the form as in Eq. (44), i.e., S(\rho )=\frac{{\hslash }^2}{6m}{\left(\frac{3{\pi }^2\rho }{2}\right)}^{2/3}+\frac{C_{s\mathrm{,}p}}{2}{\left(\frac{\rho }{{\rho }_0}\right)}^{{\gamma }_i}. The potential parameter is C_{s\mathrm{,}p}= 35.2 MeV and the symmetry energy at saturation density, S₀= 30.1 MeV. The right panel shows the dependence on {\gamma }_i of {\chi }^2 computed from the difference between predicted and measured double ratios. We determine, within a 2σ uncertainty, parameter values of 0.4⩽{\gamma }_i⩽1.05 corresponding to an increase in {\chi }^2 by 4 above its minimum near {\gamma }_i~0.7.

The density dependence of the symmetry energy has also been probed in peripheral collisions between two nuclei with different isospin asymmetries by examining the diffusion of neutrons and protons across the neck that joins them. This “isospin diffusion” generally continues until the two nuclei separate or until the chemical potentials for neutrons and protons in both nuclei become equal. Thus, in theory, the changes of isospin asymmetry of projectile-like/target-like residues at their separation time reflect the isospin diffusion ability. We named this as isospin asymmetry of the emitting source. To isolate diffusion effects from other effects, such as pre-equilibrium emission, Coulomb effects and secondary decays, measurements of isospin diffusion compare “mixed” collisions, involving a neutron-rich nucleus A and a neutron-deficient nucleus B, to the “symmetric” collisions involving A+A and B+B. The degree of isospin equilibration in such collisions can be quantified by rescaling the isospin observable X according to the definition of isospin transport ratio R_i(X) [424],

\begin{array}{}\text{(85)}& R_i(X)=\frac{2X_{AB}-X_{AA}-X_{BB}}{X_{AA}-X_{BB}}\mathrm{,}\end{array}

where X is the isospin observable. In the absence of isospin diffusion, the ratios are R=1 or R=-1. If isospin equilibrium is achieved, then the ratios R=0 for the mixed systems. Eq. (85) dictates that different observables, X, provide the same results if they are linearly related [424, 425]. The agreement of experimental isospin transport ratios obtained from isoscaling parameters, α [22] and from yield ratios of A=7 mirror nuclei [425], R_7=R(X_7=\ln (Y{(}^7\text{L}\text{i})/Y{(}^7\text{B}\text{e}))) agree, i.e., R_i(\alpha )\approx R_7, reflects nearly linear relationships between α, X₇, and the asymmetry δ of the emitting source [425].

In the ImQMD calculations, the isospin asymmetry δ of emitting source is calculated from the emitted fragments and free nucleons at the same velocity or rapidity region as in experiments [227, 424, 430]. Figure 38(a) shows the calculation results of R_i(\delta ) compared with the experimental isospin diffusion transport ratios, R_i(\alpha ), plotted as shaded regions. The ImQMD calculations are performed at impact parameters of b=5, 6, 7, and 8 fm for considering the impact parameter effects. The lines in Fig. 38(a) show the predicted isospin transport ratio R_i(\delta ) as a function of impact parameter b for {\gamma }_i=0.35, 0.5, 0.75, 1, and 2. Faster equilibration occurs for smaller {\gamma }_i values which correspond to larger symmetry energies at subsaturation densities. Thus we see a monotonic decrease of the absolute value of R_i(\delta ) with decreasing {\gamma }_i. The {\chi }^2 analysis for both impact parameters b=6 and 7 fm, which corresponds to the impact parameter region in the measurements, are performed. Using the same 2σ criterion, the analysis brackets the regions 0.45⩽{\gamma }_i⩽1.0 and 0.35⩽{\gamma }_i⩽0.8 for b= 6 and 7 fm, respectively. ImQMD calculations also provide predictions for fragment yields as a function of rapidity. The star symbols in Fig. 38(b) represent measured values of R₇ obtained from the yield ratios of ⁷Li and ⁷Be [425] at b=6 fm, as shown by the lines. This first calculation of the shapes and magnitude of the rapidity dependence of the isospin transport ratios R₇ reproduces the trends accurately. The corresponding {\chi }^2 analysis with calculations at b=6 and 7 fm favors the region 0.45⩽{\gamma }_i⩽0.95. The favorite {\gamma }_i region obtained from double neutron to proton yield ratio and isospin diffusion observables are in coincidence.

Here, we should point out the constraints on the exponent {\gamma }_i depend on the symmetry energy at saturation density, S_0=S({\rho }_0). Increasing S₀ has the same effect on the isospin transport ratio as decreasing {\gamma }_i. To compare our results to constraints obtained from nuclear masses and nuclear structure, we expand S(\rho ) around the saturation density, {\rho }_0,

\begin{array}{}\text{(86)}& \begin{array}{ll}S(\rho )=\hfill & S_0+\frac{L}{3}(\rho -{\rho }_0)/{\rho }_0+\frac{K_{sym}}{18}{[(\rho -{\rho }_0)/{\rho }_0]}^2\hfill \\ \hfill & +\frac{Q_{sym}}{162}{[(\rho -{\rho }_0)/{\rho }_0]}^3+\cdots \hfill \end{array}\end{array}

where L, K_{sym}, and Q_{sym} are the slope, curvature and skewness parameters of S(\rho ) at {\rho }_0. For realistic parameterization of S(\rho ), K_{sym} is correlated to L [431]. As the second term in Eq. (86) is much important than the third term, we believe L can be determined more reliably than K_{sym}. Furthermore, the slope parameter, L=3{\rho }_0\text{d}S(\rho )/\text{d}\rho {|}_{{\rho }_0} is related to p₀, the pressure from the symmetry energy for pure-neutron matter at saturation density. The symmetry pressure, p₀, provides the dominant baryonic contribution to the pressure in neutron stars at saturation density [432-434].

More calculations have been performed at b=6 fm with different values of {\gamma }_i and S₀ to locate the approximate boundaries in the S₀ and L plane that satisfy the 2σ criterion in the {\chi }^2 analysis of the isospin diffusion data. The two diagonal lines in Fig. 39 represent estimates in such an effort. Examination of the symmetry energy functional formed along these boundaries, where the diffusion rates are similar but S₀ and L are different, reveal that diffusion rates predominantly reflect the symmetry energy at and somewhat below.

Up to now, there are a lot of efforts on the extraction of symmetry energy. They have been made by comparing measured isospin sensitive observables, such as isospin diffusion [435] at the beam energy 74 MeV/nucleon, angular distribution of neutron-excess for light charged particles at 35 MeV/nucleon [436], collective flows [249, 426] to various transport model calculations. Consensus on the symmetry energy coefficient and slope of symmetry energy has been obtained from nuclear structure and reaction studies, as partly summarized and shown in Fig. 40, where symmetry energy coefficient S_0=30–32 MeV and slope of symmetry energy L=40–65 MeV [437, 438], but the uncertainties of symmetry energy constraints are still large.

Furthermore, the symmetry energy depends not only on S₀ and L but also on the higher terms, such as K_{sym}, Q_{sym}, …, or depends on the effective mass m^{*} and neutron proton effective mass splitting \text{Δ}m=(m_n^{*}-m_p^{*})/m. It naturally requires more accumulations of the data of the isospin sensitive observables and the development of transport models to distinguish those different physics.

5.1.2 Novel probes of symmetry energy at subsaturation denisty

The new sensitive and clear observables to the symmetry potential are also welcome for obtaining accurate knowledge of the symmetry energy. For example, some types of the direct reaction, like the elastic or inelastic scattering as well as the direct projectile breakup, involve fewer degrees of freedom in the reaction process and may reduce the difficulties in modeling the collision and could be used to constrain the symmetry energy at the subsaturation density.

Ou et al. proposed a new sensitive probe that is the isovector reorientation of deuteron induced collisions on heavy nuclei based on the ImQMD model calculations [443]. It is demonstrated in Ref. [443] that in the deuteron induced peripheral collisions on heavy nuclei such as ¹²⁴Sn, the loosely bound deuteron break up into neutron and proton, which moves with different directions relative to the incoming beam direction (named as reorientation effect), due to the isovector and Coulomb (only for proton) interaction with the target. It is found that the correlation angle determined by the relative momentum vector of the proton and the neutron originating from the breakup deuteron, which is experimentally detectable, exhibits significant dependence on the isovector nuclear potential but is robust against the variation of the isoscaler sector. Figure 41 shows the distribution of correlation angle α defined as

\begin{array}{}\text{(87)}& \cos \alpha =\frac{p_z^p-p_z^n}{|p^p-p^n|}\mathrm{,}\end{array}

for 100 MeV/nucleon (unpolarized and polarized) d+¹²⁴Sn. The calculation is performed with the ImQMD05 code and the expression of density functional is given in Section 2. The Skyrme density functionals SkA [444], SkT5, SkT1 [445], SkM* [232], Skz-1, Skz1, and Skz4 [446] are applied. The MSL0-like Skyrme interaction corresponding to the density dependent parameters {\gamma }_i=0.5 and 2. From the results of ImQMD calculations, one can see that the correlation angle α is very sensitive to the density dependence of symmetry energy, and insignificantly depends on the variation of isoscalar potentials.

Figure 42(a) shows the logarithmic spectra of cos\alpha, i.e., \ln \frac{\text{d}\sigma }{\text{d}\cos \alpha }, for various {\gamma }_i with impact parameters b=7.0 fm for the reaction of \overrightarrow{d}{+}^{124}Sn. The slope of function \ln \frac{\text{d}\sigma }{\text{d}\cos \alpha }=a_0+a_1\cos \alpha, which fitting the the logarithmic spectra of \cos \alpha in \left|\cos \alpha \right|⩽0.2, are shown in Fig. 42(b). Figure 42 clearly demonstrates the sensitivity of the correlation angle distribution to the density dependence of the symmetry energy and thus can be taken as a more clear probe to the density dependence of symmetry energy at subsaturation densities. In terms of sensitivity and cleanness, the breakup reactions induced by the polarized deuteron beam at about 100 MeV/nucleon could provide a more stringent constraint to the symmetry energy at subsaturation densities.

Another effort we have done is that we proposed the angular distribution anisotropy of coalescence invariant neutron to proton yield ratio to probe the symmetry energy at subsaturation density. The mass asymmetry reaction systems, i.e., 40Ar+¹⁹⁷Au, are analyzed, since this system has a gradient of Coulomb field during the reaction and the effects of symmetry potential at forward and backward region may be different.

The angular distribution anisotropy of coalescence invariant neutron to proton yield ratio is defined as

\begin{array}{}\text{(88)}& CIR(\text{n}/\text{p})=\frac{\text{d}{M}_{n\mathrm{,}CIR}}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}/\frac{\text{d}{M}_{p\mathrm{,}CIR}}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}\mathrm{,}\end{array}

where

\begin{array}{}\text{(89)}& \frac{\text{d}{M}_{n\mathrm{,}CIR}}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}={\displaystyle \sum _{Z\mathrm{,}N}}N\frac{Y(N\mathrm{,}Z)}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}\mathrm{,}\end{array}

\begin{array}{}\text{(90)}& \frac{\text{d}{M}_{p\mathrm{,}CIR}}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}={\displaystyle \sum _{Z\mathrm{,}N}}Z\frac{Y(N\mathrm{,}Z)}{\text{d}{\theta }_{c\mathrm{.}m\mathrm{.}}}\mathrm{.}\end{array}

As shown in Fig.43, the CIR(\text{n}/\text{p}) ratios obtained with {\gamma }_i = 0.5 are greater than that obtained with {\gamma }_i = 2.0 for the beam energies we studied, because the symmetry energy for {\gamma }_i = 0.5 is stronger than that for {\gamma }_i = 2.0 at subsaturation density. Furthermore, the CIR(\text{n}/\text{p}) ratios show a different {\theta }_{c\mathrm{.}m\mathrm{.}} dependence for the different forms of symmetry potential. For {\gamma }_i = 0.5, the CIR(\text{n}/\text{p}) ratios slightly increase as a function of {\theta }_{c\mathrm{.}m\mathrm{.}}. It is the result of the isospin asymmetry changing from {\delta }_{proj} = 0.10 at the projectile region to {\delta }_{tar} = 0.19 at the target region since the single particle potentials of neutrons and protons are close. However, for {\gamma }_i = 2.0, the CIR(\text{n}/\text{p}) ratios obviously decrease with {\theta }_{c\mathrm{.}m\mathrm{.}}, and the CIR(\text{n}/\text{p}) ratios at backward regions are smaller than that at forward regions. Our calculations show that this behavior also exists at beam energy for 50 and 100 MeV/nucleon.

5.2 Symmetry energy at suprasaturation density

The symmetry energy and its constraints at suprasaturation density are very important for understanding the dense neutron-rich matter, such as the structure of neutron star and its cooling mechanism. In this section we will discuss the efforts on constraining the symmetry energy at suprasaturation densities based on the UrQMD model calculations. Some sensitive observables to the density dependence of the symmetry energy at suprasaturation densities are discussed and then the constraints on the density dependence of the symmetry energy at suprasaturation densities are deduced by comparing the UrQMD model calculations with the measurements of flow.

In Refs. [73, 74], a “hard” Skyrme-type EoS (K=300 MeV) without momentum dependence is used in the calculations, and it was demonstrated by the UrQMD model calculations that several observables, such as the transverse momentum distribution of free neutron to proton yield ratio, {}^3\text{H} to {}^3\text{He} ratio, {\pi }^{-} to {\pi }^{+}, in neutron-rich reaction systems like ²⁰⁸Pb + ²⁰⁸Pb, ¹³²Sn+¹²⁴Sn, ⁹⁶Zr +⁹⁶Zr are sensitive to the density dependence of symmetry potential energy. In the UrQMD calculations, a momentum independent symmetry potential is used, i.e., v_{sym}^q=\frac{\partial (\rho {\delta }^2{E}_{sym}^{pot})}{\partial {\rho }_q}. Here, E_{sym}^{pot} is the corresponding density dependence of the symmetry potential energy [in Eq. (43), it is the potential part of S(\rho )]. It is written as

\begin{array}{}\text{(91)}& E_{sym}^{pot}(u)=\left(S_0-\frac{{ϵ}_F}{3}\right)F(u)\mathrm{,}\end{array}

where u=\rho /{\rho }_0 is the reduced nuclear density, S₀ and {ϵ}_F are the symmetry energy coefficient and Fermi energy, respectively. The forms of F(u) adopted in the calculations are (i) u^{\gamma } with γ= 1.5 (called F15 in the following text), (ii) u\frac{(a-u)}{a-1} with a= 3 (Fa3), a is the so-called reduced critical density; (iii) and (iv) so-called DDH{\rho }^{*} and DDH3\rho {\delta }^{*} symmetry potential energies, which are inspired by the relativistic mean-field calculations of DDH{\rho }^{*} and DDH3\rho {\delta }^{*} [448]. Figure 44 illustrates the form of density dependence of symmetry potential energy used in the UrQMD. Among them, F15 and DDH3\rho {\delta }^{*} correspond to the stiff symmetry potential energy and other two correspond to the soft one.

Figure 45 shows the transverse momentum distribution of the neutron to proton yield ratios [n/p or Y(n)/Y(p)] of emitted free nucleons calculated with four different forms of the symmetry potential. The n/p ratio, especially in the low transverse momentum region, depends strongly on the choice of the symmetry potential. In the low transverse momentum region, the n/p ratio with DDH{\rho }^{*} is the largest and that with F15 is the smallest. Obviously, nucleons with low transverse momenta are mainly emitted from the low density region. It follows that the deviation between Fa3 and DDH3\rho {\delta }^{*} at densities is small. Correspondingly, the n/p ratios at low transverse momenta calculated with Fa3 and DDH3\rho {\delta }^{*} are very close. For emitted nucleons with transverse momenta larger than \sim700 MeV/c, the n/p ratios in the Fa3 and DDH3\rho {\delta }^{*} cases are close to those in the DDH{\rho }^{*} and F15 cases, respectively. This is in correspondence with the fact that E_{sym}^{pot} for Fa3 (DDH3\rho {\delta }^{*}) is close to DDH{\rho }^{*} (F15) when \rho >{\rho }_0. Free nucleons with transverse momenta larger than \sim700 MeV/c are mainly squeezed out from higher densities [449]. The calculations performed in Refs. [73, 74] showed that those nucleons with high transverse momentum are emitted at early time when the density is high.

In addition, in Refs. [73-75], the yields of {\pi }^{-} and {\pi }^{+} as well as the {\pi }^{-}/{\pi }^{+} ratios were also investigated with the F15 and Fa3. Figure 46 shows the calculation results of the yields of {\pi }^{-} and {\pi }^{+} as well as the {\pi }^{-}/{\pi }^{+} ratios for central collision in ²⁰⁸Pb +²⁰⁸Pb at 0.4 GeV/nucleon. One can find from the figure that the sequence of the relative differences between calculation results of the n/p ratios at transverse momenta larger than \sim700 MeV/c with 4 different symmetry potentials is the same as that of the {\pi }^{-}/{\pi }^{+} ratios at higher p_t. It is because that emitted nucleons with high p_t are also from high density region as pions do. But, the UrQMD calculations show that the relative difference between the calculated results of {\pi }^{-}/{\pi }^{+} ratio from different symmetry energies is weaker than n/p ratio.

Up to now a lot of efforts have been made to extract the symmetry energy at suprasaturation densities by reproducing the FOPI {\pi }^{-}/{\pi }^{+} ratio. An apparent systematic discrepancies between the extracting symmetry energy from extreme soft to extreme stiff appeared by different transport models calculations [450-456]. Also, the calculations of Hong and Danielewicz [456] showed the {\pi }^{-}/{\pi }^{+} ratios are independent of the form of density dependence of symmetry energy after including the strong pion interaction. Further work will thus be required before pion yields and yield ratios can be reliably applied to the investigation of the high-density symmetry energy.

Flow observables have been proposed by several groups as probes for the symmetry energy at high densities [69, 81]. In Ref. [429], the UrQMD model calculations were performed with the SM EOS (soft EOS and with momentum dependent interaction) in isoscalar part, and power law form of symmetry energy as in Eq. (44) in isovector part. The results showed that the elliptic-flow ratio of neutrons with respect to protons or light complex particles in reactions of neutron-rich systems at relativistic energies was sensitive to the strength of the symmetry term at suprasaturation densities. The comparison of existing data of the ratio between the elliptic flow of neutrons and hydrogen isotopes v_2^n/v_2^H from the FOPI/LAND experiment [82, 83] with calculations performed with the UrQMD model suggested a moderately soft to linear symmetry term, characterized by a coefficient {\gamma }_i=0.9\pm 0.4 in a power law form of density dependence of the symmetry energy as in Eq. (44).

In Ref. [255], the updated version of UrQMD model, where the Skyrme potential energy density functional was adopted in the mean field part, was applied to study the flow effect and to constrain the density dependence of the symmetry energy from the FOPI/LAND elliptic flow data. Figure 47 shows the calculation results of v_2^n/v_2^H for ¹⁹⁷Au+¹⁹⁷Au at 400 MeV/nucleon with impact parameter fm. Penal (a) shows the comparison between the calculations with different Skyrme potential energy density functional with FOPI/LAND data. The calculation was performed with the same range of laboratory angles accepted by LAND. The slope parameter of L=89\pm 45 MeV for the density dependence of the symmetry energy (2\sigma uncertainty) was extracted by fitting the FOPI/LAND data. Penal (b) shows the calculation results for the intermediate rapidity window 0.25⩽\left|y_0\right|⩽0.75, for the same impact parameter and rapidity interval but without the gate on laboratory angles. It is clearly seen that the differences of the various predictions steadily grow as one moves to the region of low transverse momentum. The v_2^n/v_2^H ratio in the rapidity window 0.25⩽\left|y_0\right|⩽0.75 seems considerably more sensitive to the density dependent symmetry energy than in the mid-rapidity interval \left|y_0\right|⩽0.5, thus offering interesting opportunities for future experiments.

Furthermore, the newly measured directed and elliptic flows of neutrons and light-charged particles for the reaction ¹⁹⁷Au+¹⁹⁷Au at 400 MeV/nucleon incident energy within the ASY-EOS experimental campaign at the GSI laboratory were applied to extract the density dependence of the symmetry energy at supra-saturation densities [249]. Figure 48 shows the ASY-EOS/GSI data and the UrQMD predictions for the elliptic flow ratio of neutrons over all charged particles v_2^n/v_2^{ch} for central ( fm) collisions of ¹⁹⁷Au+¹⁹⁷Au at 400 MeV/nucleon as a function of the transverse momentum per nucleon p_t/A. The deduced symmetry-term coefficient {\gamma }_i=0.75\pm 0.10.

It is interesting to see that all the deduced density dependence of the symmetry potential energy at supra-saturation densities from the comparison with the elliptic flow data of FOPI/LAND and ASY-EOS are in coincidence, which suggests a moderate and linear density dependence of the symmetry potential energy.

5.3 EOS and symmetry energy at finite temperature

For HICs, the compressed nuclear matters are excited and the temperature is not zero. How the EOS or the symmetry energy changes at finite temperature is very important for understanding the HICs observables and extracting the EOS or symmetric energy. In Ref. [457] the temperature and density dependence of the symmetry energy was studied. The energy per nucleon at density ρ and temperature T for pure neutron and symmetry matter is calculated by using Skyrme interactions within mean-field approach and the temperature dependence arises from the modification of the zero temperature step-like momentum distributions, which becomes Fermi–Dirac distributions.

The Skyrme density functional reads

\begin{array}{}\text{(92)}& \begin{array}{ll}H=\hfill & \frac{{\hslash }^2}{2m}[{\tau }_n+{\tau }_p]\hfill \\ \hfill & +\frac{1}{4}{t}_0[(2+x_0){\rho }^2-(2x_0+1)({\rho }_n^2+{\rho }_p^2)]\hfill \\ \hfill & +\frac{1}{24}{t}_3{\rho }^{\alpha }[(2+x_3){\rho }^2-(2x_3+1)]({\rho }_p^2+{\rho }_n^2)\hfill \\ \hfill & +\frac{1}{8}[t_1(2+x_1)+t_2(2+x_2)]\tau \rho \hfill \\ \hfill & +\frac{1}{8}[t_2(2x_2+1)-t_1(2x_1+1)]({\tau }_n{\rho }_n+{\tau }_p{\rho }_p)\mathrm{,}\hfill \end{array}\end{array}

where \rho ={\rho }_n+{\rho }_p, and \tau ={\tau }_n+{\tau }_p. {\rho }_q and {\tau }_q are calculated by {\rho }_q=\frac{1}{{(2\pi \hslash )}^3}2{\displaystyle {\int }_0^{\infty }}{n}_q(p){\text{d}}^{3}p and {\tau }_q=2{\displaystyle {\int }_0^{\infty }}(n_q(p)p^2)/{\hslash }^2{\text{d}}^{3}p/{\hslash }^3. The occupation number distribution for species q, obeys Fermi–Dirac distribution.

\begin{array}{}\text{(93)}& n_q(p)=\frac{1}{1+\exp [\beta ({ϵ}_q-{\mu }_q)]}\mathrm{,}\end{array}

{ϵ}_q and {\mu }_q are the single particle energy and chemical potential for species q, and \beta =1/(K_{B}T).

By solving above equations iteratively at any pair of {\mu }_n and {\mu }_p the proton and neutron density {\rho }_p and {\rho }_n at temperature T, the energy per nucleon in neutron matter (NM) and symmetric matter (SM) can be obtained. Then the symmetry energy at finite temperature can be calculated. Figure 49 shows the density and temperature dependence of the symmetry energy in nuclear matter at temperatures T = 0, 5, 10, and 20 MeV calculated with 9 different Skyrme interactions. The 9 subfigures in Fig. 49 are ordered according to the magnitude of the ratio of effective mass R_m^0 for the applied Skyrme interactions. The quantity R_m^0 is the value of R_m=m_1^{*}/m_0^{*} at normal density, i.e., R_m^0=R_m({\rho }_0), where the subscripts 0 and 1 indicate the isospin asymmetry δ = 0 (for SM) and 1 (for NM), respectively. Obviously, R_m is proportional to the neutron and proton effective mass splitting in asymmetric matter, which can characterize the strength of the splitting of the neutron and proton effective mass for the Skyrme interaction. Figure 49 shows that there is no obvious correlation between the trend of the density dependence of the symmetry energy in cold matter and the magnitude of the R_m^0 of the corresponding Skyrme interactions. Concerning the temperature dependence of the symmetry energy, one sees that the nuclear symmetry energy decreases with increasing temperature at all densities for Skyrme interactions with small R_m^0, such as SLy7, Skz3, Skt5, and SkMP, in agreement with the results given in Refs. [458, 459]. However, for SkM, SKXm, SkP, and v070 with large R_m^0, the symmetry energy decreases with increasing temperature at low density and increases with temperature when the density is higher than a certain density. We call this phenomenon as the transition of the temperature dependence of the symmetry energy (TrTDSE). It shows that the Skyrme interactions, such as SkM, SKXm, SkP, and v070 for which the TrTDSE phenomenon occurs all satisfy m_n^{*}>m_p^{*}. The density for the onset of the TrTDSE depends on the magnitude of the R_m^0 of the corresponding Skyrme interactions. The larger R_m^0 is, the lower the density for the TrTDSE onset is. For the Skyrme interactions not satisfying m_n^{*}>m_p^{*}, the TrTDSE phenomenon will not occur.

5.4 Uncertainties in symmetry energy constraints

5.4.1 Uncertainties of the symmetry energy associated parameters

With the progress of the study of the density dependence of the symmetry energy, tight constraints become urgently requisite, which need the efforts from experimental measurements, improvements of transport models, and understanding the physics parameter correlations. For example, we know that there are not only S₀ and L in the Taylor expansion of S(\rho ) as in Eq. (86), but also the high order terms, such as curvature K_{sym} and skewness Q_{sym} of S(\rho ). Obviously, the uncertainties of K_{sym} and Q_{sym} can influence the constraints on S₀ and L, or the density dependence of the symmetry energy. There are also many efforts to constrain K_{sym} and Q_{sym} from neutron skin and neutron star [460-463].

Margueron et al.’s calculations show that the simple Taylor expansion of the EOS cannot be used to reproduce the EOS well at the whole density region as well as for the symmetry energy, and they proposed a meta-EOS model to describe it [461-463]. Another method to well describe the Skyrme EOS and symmetry energy is to use the nuclear matter parameters, such as {\rho }_0, E₀, K₀, S₀, L, m_s^{*}, m_v^{*}, with two additional coefficients g_{sur} and g_{sur\mathrm{,}iso} [161, 464, 465]. Here, {\rho }_0 is the normal density, K_0=9{\rho }_0\frac{{\partial }^2ϵ/\rho }{\partial {\rho }^2}{|}_{{\rho }_0} is the incompressibility of symmetric nuclear matter, m_s^{*}/m=\left(1+\frac{2m}{{\hslash }^2}\frac{\partial }{\partial \tau }\frac{E}{A}\right){|}_{{\rho }_0} is the isoscalar effective mass, m_v^{*}=\frac{1}{1+\kappa } is the isovector effective mass where κ is the enhancement of a factor of the Thomas–Reich–Kuhn sum rule. g_{sur} and g_{sur\mathrm{,}iso} are the coefficients related to density gradient terms.

A lot of theoretical works have evidenced that all of them are related to the symmetry energy. For example, in the Skyrme–Hartree–Fock approaches, the density dependence of symmetry energy is written as

\begin{array}{}\text{(94)}& \begin{array}{ll}S(\rho )=\hfill & \frac{1}{3}\frac{{\hslash }^2}{2m}{\left(\frac{3{\pi }^2}{2}\rho \right)}^{2/3}\hfill \\ \hfill & +(A_{sym}u+B_{sym}{u}^{\eta }+C_{sym}(m_s^{*}\mathrm{,}{m}_v^{*})u^{5/3})\mathrm{,}\hfill \end{array}\end{array}

where u is the reduced density, i.e., \rho /{\rho }_0. A recent theoretical study by Mondal et al. also provided an evidence that S(\rho ) depends on the effective mass [466]. In the standard Skyrme interaction, one also observed that m_s^{*} is also related to K₀ based on the formula of Skyrme Hartree–Fock (SHF) as pointed out in Ref. [392],

\begin{array}{}\text{(95)}& K_0=B+C\sigma +D\left(1-\frac{3}{2}\sigma \right)\frac{8{\hslash }^2}{m{\rho }_0}\left(\frac{m}{m_s^{*}}-1\right)\mathrm{,}\end{array}

with B=-9E_0+\frac{3}{5}{ϵ}_F, C=-9E_0+\frac{9}{5}{ϵ}_F and D=\frac{3}{20}{\rho }_0{k}_F^2. If E₀ and {\rho }_0 are well known, K₀ depends on the m_s^{*} and σ. Focusing on the correlation between m_s^{*} and K₀, one can say K₀ is independent of m_s^{*} if \sigma =2/3, but K₀ linearly depends on the inverse of m_s^{*} if \sigma \ne \frac{2}{3}. Thus, one can expect that the constraint of S(\rho ) with less biased uncertainty should depend on the values of {\rho }_0, E₀, K₀, S₀, L, m_s^{*}, m_v^{*} rather than only on the uncertainties of S₀ and L.

Figure 50 shows the values of nuclear matter parameters, K_0, S_0, L, m_s^{*}/m, and f_I calculated from 224 effective Skyrme interactions published from the years 1970–2015 [467]. Here, f_I=\frac{1}{2\delta }\left(\frac{m}{m_n^{*}}-\frac{m}{m_p^{*}}\right)=\frac{m}{m_s^{*}}-\frac{m}{m_v^{*}}, it can be analytically incorporated into the transport model and its sign reflects m_n^{*}>m_p^{*} or . The nuclear matter incompressibility from Skyrme parameter sets converges to the region of 200–280 MeV after the year ~1990, except for the parameter from the original quark meson coupling (QMC) method [468] (red circles in upper panels of Fig. 50) which were readjusted in 2006. In Ref. [469], they show the value of K₀ is higher than generally accepted by a considerable margin, i.e., K_0=240\pm 20 MeV, based on the most precise and up-to-date data on GMR energies of Sn and Cd isotopes, together with a selected set of data from ⁵⁶Ni to ²⁰⁸Pb. This result is MeV which has been obtained without any microscopic model assumptions, except (marginally) the Coulomb effect, and revealed the essential role of surface properties in vibrating nuclei.

For other nuclear matter parameters, such as S₀, L, m_s^{*}, and f_I, most of their values fall into the regions of S_0=25–35 MeV, L=30–120 MeV, m_s^{*}/m=0.6-1.0, f_I=-0.5–0.4. The very recent results on the estimated nuclear matter parameters [461] are shown as black squares in Fig. 50.

In the simulations of neutron-rich HICs, the calculations results of the isospin sensitive observables could be influenced by those nuclear matter parameters, such as effective mass splitting and isoscalar effective mass, which makes the effect of the symmetry energy in a certain extent [161, 471-473]. Thus, the effects of different K₀, m_s^{*} and different effective mass splitting should also be investigated and make unentanglement with the effect of the symmetry potential on the isospin sensitive observables in HICs.

5.4.2 Influence of effective mass splitting on R_i, R_i(y), and n/p ratios

In the investigations with the standard Skyrme interaction in HICs, four Skyrme interaction parameter sets, SLy4, SkI2, SkM*, and Gs [232, 392, 474, 475] which have similar incompressibility (K₀), symmetry energy coefficient (S₀) i.e., K_0=230\pm 20 MeV, S_0=32\pm 2 MeV, but different isoscalar effective mass m^{*}/m=0.7\pm 0.1 and different effective mass splitting, are adopted in the ImQMD-Sky calculations. The SLy4 and SkI2 [392, 474] have similar neutron/proton effective mass splitting with , but very different slopes of symmetry energy L values, 46 MeV for SLy4, and 104 MeV for SkI2. The other two Skyrme interaction parameter sets with m_n^{*}>m_p^{*} and m^{*}/m\sim 0.78 also have different L values, 46 MeV for SkM* [232], and 93 MeV for Gs [475]. The saturation properties of nuclear matter for these four Skymre interactions are listed in Table 4. By analyzing the results calculated with these interactions, the sensitivities of the isospin observables on the different nuclear matter parameters can be investigated.

We simulated the collisions of ¹²⁴Sn+¹²⁴Sn, ¹²⁴Sn+ ¹¹²Sn, ¹¹²Sn+¹²⁴Sn, and ¹¹²Sn+¹¹²Sn reactions at beam energy of 50 MeV/nucleon using the ImQMD-Sky code. 64 000 events were performed for each reaction at each impact parameter. Previous theoretical studies [160, 476] and recent experimental studies [477] suggest that there is no strong dependence of transverse emitted neutron to proton yield ratios on the impact parameter. In Fig. 51(a), we plot the isospin transport ratios obtained with SLy4, SkI2, SkM* and Gs interactions at b=6 fm. As in previous studies [430, 160], we analyze the amount of isospin diffusion by constructing a tracer, X=\delta, from the isospin asymmetry of emitting source which includes all emitted nucleons (N) and fragments (frag) with the velocity cut (v_z^{N\mathrm{,}frag}>0.5v_{beam}^{c\mathrm{.}m\mathrm{.}}). The shaded region is experimental data obtained by constructing the isospin transport ratio using isoscaling parameter X=\alpha, near the projectile rapidity regions [424]. Our results show that the isospin transport ratio R_i [see Eq. (85)] values for SLy4 (solid circle) and SkM* (solid squares), both with L=46 MeV, lie within the experimental uncertainties while R_i values for SkI2 (open circle, L=104 MeV) and Gs (open square, L=93 MeV) are above the data range. Even though the isospin diffusion process is accelerated at subsaturation densities with the stronger Lane potential, the overall effect of mass splitting on isospin diffusion is small. This conclusion is similar with previous results even from the IBUU and SMF models [70, 471]. Since the isospin diffusion process is strongly related to the difference of local isospin concentration, the strong repulsive momentum-dependent isoscalar potential reduces the effect of isovector potential on the reaction dynamics, thus, there is no clear pattern that R_i values decrease significantly with the strength of Lane potentials.

We also compare results of the calculations to R_i as a function of the scaled rapidity y/y_{beam}^{c\mathrm{.}m\mathrm{.}} as shown in the right panel of Fig. 51. The star symbols in the right panel are experimental data, R_i(X_7)=R_7, obtained in Ref. [425]. This transport ratio was generated using the isospin tracer X=\ln [Y{(}^7\text{L}\text{i})/Y(\text{B}\text{e})], where Y{(}^7\text{L}\text{i})/Y{(}^7\text{B}\text{e}) is the yield ratio of the mirror nuclei, ⁷Li and ⁷Be [425]. For comparison, the ImQMD-Sky calculations of R_i are plotted as lines for b=6 fm. The interactions with smaller L values, SLy4 and SkM* (solid lines) agree with the data better especially in the high rapidity region. However, {\chi }^2 analysis suggests that the quality of fit with isospin diffusion data is not good enough to draw definite conclusions about mass splitting effect with confidence. We need a more sensitive and reliable observable to extract quantitative information about the nucleon effective mass splitting.

The calculated results on single ratio R(n/p) and double ratio DR(n/p) are shown in Fig. 52. We plot the Y(n)/Y(p) ratios as a function of kinetic energy of emitted nucleon in center of mass frame, E_{c\mathrm{.}m\mathrm{.}}, for ¹¹²Sn+¹¹²Sn (left panel) and ¹²⁴Sn+¹²⁴Sn (middle panel) at b=2 fm with angular gate 70°<θ_c.m.<110°. The lines connecting the circles correspond to case, and the lines connecting the squares correspond to m_n^{*}>m_p^{*} case. Not surprisingly, the Y(n)/Y(p) ratios are larger for the neutron rich system, ¹²⁴Sn+¹²⁴Sn, in the middle panel. Consistent with Refs. [471, 472], the differences in the Y(n)/Y(p) ratios between them (circles) and m_n^{*}>m_p^{*} (squares) increase with nucleons kinetic energy. At high nucleon energies, the stronger Lane potentials with enhance neutron emissions, leading to flatter Y(n)/Y(p) dependence on the nucleon kinetic energy. The calculations with SLy4 () are consistent with the double ratios data from Ref. [423] which was published in 2006, especially at high kinetic energy region.

Recently, remeasurements of the neutron to proton yield data have been finished and were published in Ref. [478]. The measured coalescence invariant spectral double ratios DR(n/p) for both beam energies, 50 and 120 MeV/nucleon, were analyzed. There are systematic uncertainties in DR(n/p) of about 10% at E_{beam} = 50 MeV/nucleon and 15% at E_{beam} = 120 MeV/nucleon stemming from the dependence of the neutron detection efficiencies on the charged particle and neutron-scattering backgrounds in LANA. The current data at E_{beam} = 50 MeV/nucleon have a factor of 4 smaller uncertainties and extend over a wider energy range than those of Ref. [423]. The two data overlap within statistical and systematic uncertainties, except for the lowest two energy data points. The new precision data show there is no agreement between the data and calculations due to the inadequate description of cluster formation mechanism at the beam energy of 50 MeV/nucleon. At higher incident energy, E_{beam} = 120 MeV/nucleon, the calculation describes the nucleon spectra fairly well. Furthermore, at high kinetic energy of emitted nucleons (E_{c\mathrm{.}m\mathrm{.}}>60 MeV), the ImQMD-Sky calculations of Ref. [227] showed the greatest sensitivity to the effective mass splitting, and calculations with SLy4 interaction lie close to the data whereas the calculations with SkM* lie below the data.

Further theoretical analyses have indicated that n/p yield ratios (as well as other observables) are also somewhat sensitive to the other nuclear matter parameters. For example, in Ref. [161], we discussed the linear correlation coefficient C_{AB} between the nuclear matter parameter and the isospin sensitive observable in the sampled parameter space. The linear-correlation coefficient C_{AB} between variable A and observable B is calculated as follows [479]:

\begin{array}{}\text{(96)}& C_{AB}=\frac{cov(A\mathrm{,}B)}{\sigma (A)\sigma (B)}\mathrm{,}\end{array}

\begin{array}{}\text{(97)}& cov(A\mathrm{,}B)=\frac{1}{N-1}{\displaystyle \sum _i}(A_i-\langle A\rangle )(B_i-\langle B\rangle )\mathrm{,}\end{array}

\begin{array}{}\text{(98)}& \sigma (X)=\sqrt{\frac{1}{N-1}{\displaystyle \sum _i}{(X_i-\langle X\rangle )}^2}\mathrm{,}X=A\mathrm{,}B\mathrm{,}\end{array}

\begin{array}{}\text{(99)}& \langle X\rangle =\frac{1}{N}{\displaystyle \sum _i}{X}_i\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}N\mathrm{.}\end{array}

cov(A\mathrm{,}B) is the covariance, \sigma (X) is the variance. C_{AB}=\pm 1 means there is a linear dependence between A and B, and C_{AB}=0 means no correlations. As shown in the panel (e) and (j) of Fig. 53, isospin diffusion also related to the isoscalar effective mass [161], which may have to be better constrained in order to accurately determine the effective mass splitting and slope of symmetry energy. Those studies also stimulate the further statistical analysis or Bayesian analysis in the multi-dimension parameter surface with respect to the neutron to proton yield data [480].

5.4.3 Influence of K₀, S₀, L, m_s^{*}, and f_I on isospin diffusion

In order to investigate the impact of other nuclear matter parameters on the isospin diffusion observable, we calculate it in five-dimensional (5D) parameter space, such as K₀, S₀, L, m_s^{*}, f_I, with ImQMD-Sky. We sampled 120 points in the range which we listed in Table 5 under the condition that \eta ⩾1.1. \eta ⩾1.1 is used for guaranteeing the reasonable three-body force in the transport model calculations. The ranges of these nuclear matter parameters are chosen based on the prior information of Skyrme parameters as shown in Fig. 50. As an example, 120 sampled points are presented as open and solid circles in two-dimensional projection in Fig. 54. The points of parameter sets uniformly distribute in two-dimensional projection except for the plots of K₀ and m_s^{*}/m due to the restriction of \eta ⩾1.1. We perform the calculations for isospin transport diffusion at 35 MeV/nucleon and 50 MeV/nucleon at b=5–8 fm with the impact parameter smearing [481] for ^112,124Sn+^112,124Sn. 10 000 events are calculated for each point in the parameter space and simulations are stopped at 400 fm/c. The calculations are performed on TianHe-1 (A), the National Supercomputer Center in Tianjin.

In Fig. 55, the lines represent the calculated results of the isospin transport ratio R_i with 120 parameter sets. Two stars are the experimental data [424, 425, 482] which is constructed from the isoscaling parameter X={\alpha }_{iso} at 50 MeV/nucleon [424] and the ratio of X=\ln (Y{(}^7\text{L}\text{i})/Y{(}^7\text{B}\text{e})) [482, 425] at the beam energy of 35 MeV/nucleon, which was assumed and evidenced to linearly related to the isospin asymmetry of emitting source [425]. And thus, one can compare R_i(\delta ) to R_i(\alpha ) or R_i(\ln (Y{(}^7\text{L}\text{i})/Y{(}^7\text{B}\text{e})). As shown in Fig. 55, the calculated results show a large spread around the experimental data. By comparing the calculations to the data, we find 22 parameter sets that can reproduce the isospin diffusion data within experimental errors. The extracted 22 parameter sets are listed in Table 6. We highlight those points that can reproduce the experimental data within experimental errors with blue colors in Fig. 54. Generally, one can observe that L increases with S₀. The constrained points distribute in the bottom-right corner in the S₀–L plot [panel (c)], and the large L with small S₀ are ruled out.

The results in panel (j) of Fig. 54 show that isospin diffusion data is not sensitive to the effective mass and its splitting. Based on Fig. 54, one can learn that there are no obvious correlations between R_i and S₀, L, m_s^{*}/m, and f_I in 5D parameter space. It is because R_i is not only correlated to L but also correlated to m_s^{*}/m [161], which broke the R_i dependence of L when we randomly chose the values of K₀, S₀, L, m_s^{*}/m, and f_I. If we fix the values of K₀, S₀, m_s^{*}/m, and f_I, the positive correlation between R_i and L can be found.

6 Discussion and prospect

The transport model has been widely and successfully applied in the study of heavy ion collisions from low energies to relativistic energies, for understanding the nuclear phenomena, such as collective flow, particle emission, fussion-fission, multifragmentation mechanism and the properties of nuclear equation of state in the laboratory.

However, the model uncertainties still need to be understood for reliable extracting the physics information of interest. In last 30 years, the transport code comparison project has been performed in the community to seek the reasons [229, 378, 483, 484] and reduce the model uncertainties. Despite a lot of efforts devoted to the improvement of the model, we still have to face some fundamental problems for the demands of more exact theoretical description and accurate calculation of heavy ion collisions.

Because of the complexity of transport equations, and in particular their dimensionality, the collision term and mean field potential term are treated separately and it causes the model uncertainties in solving the transport equation [484] and extracting the properties of nuclear system. Ideally, the final but very difficult goal for transport model is to develop a precise and accurate enough transport model, which takes enough quantum effects and treats the mean field and collision self-consistently for the quantum N-body system.

The quantum effects in transport model play more important roles for heavy ion reaction at low energy than that at high energy, but this is not well considered in the current models. One of the reason is that the coarse-grain method is adopted to derive the transport equation. It leads to a semi-classical transport equation, where the phase space distribution f will finally evolve to the Boltzmann distribution. It is mimicked by the Pauli blocking, but it is always underestimated in the current transport models, especially for low energy heavy ion reaction [378]. It seems to us that there is a long standing task to improve the treatment of the quantum effects, such as the shell effect which plays an important role in some cases as well as the Pauli blocking effect for the better description of low energy heavy ion collision by transport models. However, it is also needed to go beyond the mean field for heavy ion reaction at low energy.

Fluctuations in transport theory are the main venue to go beyond dissipative mean field dynamics, and it plays an important role in the large amplitude motions of nuclear systems, such as nuclear fission, fusion, MNT, ternary break, and multifragmentation reactions as one observed in experiments of heavy ion collisions. As an example, in Ref. [485], the models of MNT reaction were tested by reaction ¹³⁶Xe+¹⁹⁸Pt and it turned out that the ImQMD calculations provide better agreement with experimental data, and it seems to us that it is due to the stronger fluctuation and dissipation being considered dynamically in the ImQMD model. To understand and properly treat the fluctuation and dissipation in transport models is of importance for a better description of the fusion, MNT as well as the multifragmentation, etc., in heavy ion collisions from both theoretical study and practical applications.

The cluster formation mechanism is another important issue which needs to be solved in transport models, since the heavy ion collision observables used to extract the physics information of interest are obtained from the measurement of the emitted light particles or fragments. In physics, the cluster formation is related to the fluctuation when the system enters into the spinodal region as well as the N-body correlations. But both mechanisms are hard to accurately deal with and to incorporate in the calculations for describing the light particles productions in the mean field level. For heavy ion reactions at intermediate energy, one has clearly observed the multifragmentation phenomena, and there are lots of light charged particles formed. But in the transport model, except for the AMD considering the NN correlation to form light particles [486], the calculations always overpredict the Z=1 particles and underestimate the Z=2 particles. It naturally requires including a reasonable cluster formation mechanism besides the current Hamiltonian dynamics we used.

Baryons and mesons production and propagation mechanism near the threshold energy are also important issues for extracting the information of EOS at high density. At this energy region, the issue related to a covariant dynamics model in the propagation and collisions part, in-medium cross sections of NN\leftrightarrow NR/NN\leftrightarrow RR, energy conservation in the process of NN\leftrightarrow NR/NN\leftrightarrow RR/R\leftrightarrow N\pi (R is a resonance particles) and the potential for mesons, should be well investigated.

References

[1] Gaidos J. A., Gutay L. J., Hirsch A. S., Mitchell R., Ragland T. V., Scharenberg R. P., Turkot F., Willmann R. B., and Wilson C. L., Nuclear fragment emission in high-energy p–Xe and p–Kr collisions and a description of their production mechanism, Phys. Rev. Lett. 42(2), 82 (1979)

[2] Minich R. W., Agarwal S., Bujak A., Chuang J., Finn J. E., Gutay L. J., Hirsch A. S., Porile N. T., Scharenberg R. P., Stringfellow B. C., and Turkot F., Critical phenomena in hadronic matter and experimental isotopic yields in high energy proton–nucleus collisions, Phys. Lett. B 118(4–6), 458 (1982)

[3] Finn J. E., Agarwal S., Bujak A., Chuang J., Gutay L. J., Hirsch A. S., Minich R. W., Porile N. T., Scharenberg R. P., Stringfellow B. C., and Turkot F., Nuclear fragment mass yields from high-energy proton–nucleus interactions, Phys. Rev. Lett. 49(18), 1321 (1982)

[4] Landau L. D. and Lifshitz E. M., Fluid Mechanics, London: Pergamon Press, 1959

[5] Scheid W., Ligensa R., and Greiner W., Ion–ion potentials and the compressibility of nuclear matter, Phys. Rev. Lett. 21(21), 1479 (1968)

[6] Hofmann J., Stöcker H., Heinz U., Scheid W., and Greiner W., Possibility of detecting density isomers in high-density nuclear mach shock waves, Phys. Rev. Lett. 36(2), 88 (1976)

[7] Siemens P. J., Heavy ion collisions, Nucl. Phys. A 335(1–2), 491 (1980)

[8] Gustafsson H. A., Gutbrod H. H., Kolb B., Löhner H., Ludewigt B., Poskanzer A. M., Renner T., Riedesel H., Ritter H. G., Warwick A., Weik F., and Wieman H., Collective flow observed in relativistic nuclear collisions, Phys. Rev. Lett. 52(18), 1590 (1984)

[9] Stöcker H. and Greiner W., High energy heavy ion collisions — probing the equation of state of highly excited hardronic matter, Phys. Rep. 137(5–6), 277 (1986)

[10] Clare R. B. and Strottman D., Relativistic hydrodynamics and heavy ion reactions, Phys. Rep. 141(4), 177 (1986)

[11] Stöcker H., Maruhn J. A., and Greiner W., Collective sideward flow of nuclear matter in violent high-energy heavy-ion collisions, Phys. Rev. Lett. 44(11), 725 (1980)

[12] Renfordt R. E., Schall D., Bock R., Brockmann R., Harris J. W., Sandoval A., Stock R., Ströbele H., Bangert D., Rauch W., Odyniec G., Pugh H. G., and Schroeder L. S., Stopping power and collective flow of nuclear matter in the reaction Ar+Pb at 0.8 GeV/u, Phys. Rev. Lett. 53(8), 763 (1984)

[13] Ströbele H., Brockmann R., Harris J. W., Riess F., Sandoval A., Stock R., Wolf K. L., Pugh H. G., Schroeder L. S., Renfordt R. E., Tittel K., and Maier M., Charged-particle exclusive analysis of central Ar + KCl and Ar + Pb reactions at 1.8 and 0.8 GeV/nucleon, Phys. Rev. C 27(3), 1349 (1983)

[14] Baden A., Gutbrod H. H., Löhner H., Maier M. R., Poskanzer A. M., Renner T., Riedesel H., Ritter H. G., Spieler H., Warwick A., Weik F., and Wieman H., The plastic ball spectrometer: An electronic 4τ detector with particle identification, Nucl. Instrum. Methods 203(1–3), 189 (1982)

[15] Buchwald G., Graebner G., Theis J., Maruhn J., Greiner W., and Stöcker H., Kinetic energy flow in Nb(400 A MeV) + Nb: Evidence for hydrodynamic compression of nuclear matter, Phys. Rev. Lett. 52(18), 1594 (1984)

[16] Reisdorf W. and Ritter H. G., Collective flow in heavy-ion collisions, Annu. Rev. Nucl. Part. Sci. 47(1), 663 (1997)

[17] Herrmann N., Wessels J. P., and Wienold T., Collective flow in heavy-ion collisions, Annu. Rev. Nucl. Part. Sci. 49(1), 581 (1999)

[18] Wada R., Hildenbrand K. D., Lynen U., Müller W. F. J., Rabe H. J., et al., Isotopic-yield ratios of complex fragments from intermediate-energy heavy-ion reactions, Phys. Rev. Lett. 58(18), 1829 (1987)

[19] Johnston H., White T., Winger J., Rowland D., Hurst B., Gimeno-Nogues F., O’Kelly D., and Yennello S. J., Isospin equilibration in the reaction ⁴⁰Ar, ⁴⁰Ca+⁵⁸Fe, ⁵⁸Ni, Phys. Lett. B 371(3–4), 186 (1996)

[20] Veselsky M., Ibbotson R. W., Laforest R., Ramakrishnan E., Rowland D. J., Ruangma A., Winchester E. M., Martin E., and Yennello S. J., Inhomogeneous isospin distribution in the reactions of ²⁸Si+¹¹²Sn and ¹²⁴Sn at 30 and 50 MeV/nucleon, Phys. Rev. C 62(4), 041605 (2000)

[21] Mekjian A. Z., Explosive nucleosynthesis, equilibrium thermodynamics, and relativistic heavy-ion collisions, Phys. Rev. C 17(3), 1051 (1978)

[22] Tsang M. B., Friedman W. A., Gelbke C. K., Lynch W. G., Verde G., and Xu H. S., Isotopic scaling in nuclear reactions, Phys. Rev. Lett. 86(22), 5023 (2001)

[23] Xu H., Alfaro R., Davin B., Beaulieu L., Larochelle Y., et al., Fragment isospin as a probe of heavy-ion collisions, Phys. Rev. C 65, 061602(R) (2002)

[24] Danielewicz P., Lacey R., and Lynch W. G., Determination of the equation of state of dense matter, Science 298(5598), 1592 (2002)

[25] Adamian G. G., Antonenko N. V., and Scheid W., Model of competition between fusion and quasifission in reactions with heavy nuclei, Nucl. Phys. A 618(1–2), 176 (1997)

[26] Diaz-Torres A., Adamian G. G., Antonenko N. V., and Scheid W., Quasifission process in a transport model for a dinuclear system, Phys. Rev. C 64(2), 024604 (2001)

[27] Wang N., Zhao E. G., Scheid W., and Zhou S.-G., Theoretical study of the synthesis of superheavy nuclei with Z=119 and 120 in heavy-ion reactions with trans-uranium targets, Phys. Rev. C 85, 041601(R) (2012)

[28] Zagrebaev V. and Greiner W., Unified consideration of deep inelastic, quasi-fission and fusion-fission phenomena, J. Phys. G 31(7), 825 (2005)

[29] Nasirov A. K., Mandaglio G., Giardina G., Sobiczewski A., and Muminov A. I., Effects of the entrance channel and fission barrier in the synthesis of superheavy element Z=120, Phys. Rev. C 84(4), 044612 (2011)

[30] Hofmann S. and Münzenberg G., The discovery of the heaviest elements, Rev. Mod. Phys. 72(3), 733 (2000)

[31] Oganessian Yu. Ts., Abdullin F. S., Bailey P. D., Benker D. E., Bennett M. E., et al., Synthesis of a new element with atomic number Z=117, Phys. Rev. Lett. 104(14), 142502 (2010)

[32] Gupta R. K., Manhas M., Münzenberg G., and Greiner W., Theory of the compactness of the hot fusion reaction ⁴⁸Ca+²⁴⁴Pu→²⁹²114*, Phys. Rev. C 72, 014607 (2005)

[33] Adamian G. G., Antonenko N. V., Diaz-Torres A., and Scheid W., Dynamical restriction for a growing neck due to mass parameters in a dinuclear system, Nucl. Phys. A 671(1–4), 233 (2000)

[34] Lu B. N., Zhao E. G., and Zhou S. G., Potential energy surfaces of actinide nuclei from a multidimensional constrained covariant density functional theory: Barrier heights and saddle point shapes, Phys. Rev. C 85, 011301(R) (2012)

[35] Wang N., Tian J. L. and Scheid W., Systematics of fusion probability in “hot” fusion reactions, Phys. Rev. C 84, 061601(R) (2011)

[36] Sargsyan V. V., Adamian G. G., Antonenko N. V., Scheid W., and Zhang H. Q., Effects of nuclear deformation and neutron transfer in capture processes, and fusion hindrance at deep sub-barrier energies, Phys. Rev. C 84(6), 064614 (2011)

[37] Dasgupta M., Hinde D. J., Diaz-Torres A., Bouriquet B., Low C. I., Milburn G. J., and Newton J. O., Beyond the coherent coupled channels description of nuclear fusion, Phys. Rev. Lett. 99(19), 192701 (2007)

[38] Leigh J. R., Dasgupta M., Hinde D. J., Mein J. C., Morton C. R., Lemmon R. C., Lestone J. P., Newton J. O., Timmers H., Wei J. X., and Rowley N., Barrier distributions from the fusion of oxygen ions with ^144,148,154Sm and ¹⁸⁶W, Phys. Rev. C 52(6), 3151 (1995)

[39] Timmers H., Ackermann D., Beghini S., Corradi L., He J. H., Montagnoli G., Scarlassara F., Stefanini A. M., and Rowley N., A case study of collectivity, transfer and fusion enhancement, Nucl. Phys. A 633(3), 421 (1998)

[40] Zhang H. Q., Lin C. J., Yang F., Jia H. M., Xu X. X., Wu Z. D., Jia F., Zhang S. T., Liu Z. H., Richard A., and Beck C., Near-barrier fusion of ³²S+^90,96Zr: The effect of multi-neutron transfers in sub-barrier fusion reactions, Phys. Rev. C 82(5), 054609 (2010)

[41] Sobiczewski A. and Pomorski K., Description of structure and properties of superheavy nuclei, Prog. Part. Nucl. Phys. 58(1), 292 (2007)

[42] Myers W. D. and Swiatecki W. J., Nucleus–nucleus proximity potential and superheavy nuclei, Phys. Rev. C 62(4), 044610 (2000)

[43] Liu M., Wang N., Li Z., Wu X., and Zhao E., Applications of Skyrme energy-density functional to fusion reactions spanning the fusion barriers, Nucl. Phys. A 768(1–2), 80 (2006)

[44] Wang N., Liu M., and Yang Y. X., Heavy-ion fusion and scattering with Skyrme energy density functional, Sci. China G: Phys. Mech. Astron. 52(10), 1554 (2009))

[45] Wang N., Zhao K., Scheid W., and Wu X., Fusion-fission reactions with a modified Woods–Saxon potential, Phys. Rev. C 77(1), 014603 (2008)

[46] Zagrebaev V. and Greiner W., Synthesis of superheavy nuclei: A search for new production reactions, Phys. Rev. C 78(3), 034610 (2008)

[47] Cap T., Siwek-Wilczyńska K., and Wilczyński J., Nucleus–nucleus cold fusion reactions analyzed with the l-dependent “fusion by diffusion” model, Phys. Rev. C 83(5), 054602 (2011)

[48] Hagino K., Rowley N., and Kruppa A. T., A program for coupled-channel calculations with all order couplings for heavy-ion fusion reactions, Comput. Phys. Commun. 123(1–3), 143 (1999)

[49] Wong C. Y., Interaction barrier in charged-particle nuclear reactions, Phys. Rev. Lett. 31(12), 766 (1973)

[50] Puri R. K. and Gupta R. K., Fusion barriers using the energy-density formalism: Simple analytical formula and the calculation of fusion cross sections, Phys. Rev. C 45(4), 1837 (1992)

[51] Wang B., Wen K., Zhao W. J., Zhao E. G., and Zhou S. G., Systematics of capture and fusion dynamics in heavy-ion collisions, At. Data Nucl. Data Tables 114, 281 (2017)

[52] Zagrebaev V. I., Oganessian Yu. Ts., Itkis M. G., and Greiner W., Superheavy nuclei and quasi-atoms produced in collisions of transuranium ions, Phys. Rev. C 73, 031602(R) (2006)

[53] Zagrebaev V. and Greiner W., Low-energy collisions of heavy nuclei: Dynamics of sticking, mass transfer and fusion, J. Phys. G 34(1), 1 (2007)

[54] Zagrebaev V. and Greiner W., Shell effects in damped collisions: A new way to superheavies, J. Phys. G 34(11), 2265 (2007)

[55] Zhang F., Li C., Zhu L., and Wen P., Production cross sections for exotic nuclei with multinucleon transfer reactions, Front. Phys. 13(6), 132113 (2018)

[56] Umar A. S. and Oberacker V. E., Heavy-ion interaction potential deduced from density-constrained time-dependent Hartree–Fock calculation, Phys. Rev. C 74, 021601(R) (2006)

[57] Umar A. S., Oberacker V. E., Maruhn J. A., and Reinhard P. G., Microscopic composition of ion–ion interaction potentials, Phys. Rev. C 85(1), 017602 (2012)

[58] Nakatsukasa T. and Yabana K., Linear response theory in the continuum for deformed nuclei: Green’s function vs time-dependent Hartree–Fock with the absorbing boundary condition, Phys. Rev. C 71(2), 024301 (2005)

[59] Maruhn J. A., Reinhard P. G., Stevenson P. D., and Strayer M. R., Spin-excitation mechanisms in Skyrme-force time-dependent Hartree–Fock calculations, Phys. Rev. C 74(2), 027601 (2006)

[60] Guo L., Maruhn J. A., and Reinhard P. G., Boost-invariant mean field approximation and the nuclear Landau–Zener effect, Phys. Rev. C 76(1), 014601 (2007)

[61] Simenel C., Nuclear quantum many-body dynamics, Eur. Phys. J. A 48(11), 152 (2012)

[62] Li B. A., Chen L. W., and Ko C. M., Recent progress and new challenges in isospin physics with heavy-ion reactions, Phys. Rep. 464(4–6), 113 (2008)

[63] Li B. A., Ko C. M., and Ren Z. Z., Equation of state of asymmetric nuclear matter and collisions of neutron-rich nuclei, Phys. Rev. Lett. 78(9), 1644 (1997)

[64] Li B. A., Ko C. M., and Bauer W., Isospin physics in heavy-ion collisions at intermediate energies, Int. J. Mod. Phys. E 07(02), 147 (1998)

[65] Chen L. W., Ko C. M., Li B. A., Xu C., and Xu J., Probing isospin- and momentum-dependent nuclear effective interactions in neutron-rich matter, Eur. Phys. J. A 50(2), 29 (2014)

[66] Li B. A., Probing the high density behavior of the nuclear symmetry energy with high energy heavy-ion collisions, Phys. Rev. Lett. 88(19), 192701 (2002)

[67] Li B. A., Isospin dependence of the π−/π+ ratio and density dependence of the nuclear symmetry energy, Phys. Rev. C 67(1), 017601 (2003)

[68] Baran V., Colonna M., Greco V., and Di Toro M., Reaction dynamics with exotic nuclei, Phys. Rep. 410(5–6), 335 (2005)

[69] Greco V., Baran V., Colonna M., Di Toro M., Gaitanos T., and Wolter H. H., Relativistic effects in the search for high density symmetry energy, Phys. Lett. B 562(3–4), 215 (2003)

[70] Chen L. W., Ko C. M., and Li B. A., Determination of the stiffness of the nuclear symmetry energy from isospin diffusion, Phys. Rev. Lett. 94(3), 032701 (2005)

[71] Chen L. W., Greco V., Ko C. M., and Li B. A., Effects of symmetry energy on two-nucleon correlation functions in heavy-ion collisions induced by neutron-rich nuclei, Phys. Rev. Lett. 90(16), 162701 (2003)

[72] Gaitanos T., Di Toro M., Type S., Baran V., Fuchs C., Greco V., and Wolter H. H., On the Lorentz structure of the symmetry energy, Nucl. Phys. A 732, 24 (2004)

[73] Li Q. F., Li Z. X., Soff S., Gupta R. K., Bleicher M., and Stöcker H., Probing the density dependence of the symmetry potential at low and high densities, Phys. Rev. C 72(3), 034613 (2005)

[74] Li Q. F., Li Z. X., Soff S., Gupta R. K., Bleicher M., and Stöcker H., Probing the equation of state with pions, J. Phys. G 32(2), 151 (2006)

[75] Li Q., Li Z., Soff S., Gupta R. K., Bleicher M., and Stöcker H., Probing the density dependence of the symmetry potential in intermediate-energy heavy ion collisions, J. Phys. G 31(11), 1359 (2005)

[76] Li Q. F., Li Z. X., Zhao E. G., and Gupta R. K., Σ−/Σ+ ratio as a candidate for probing the density dependence of the symmetry potential at high nuclear densities, Phys. Rev. C 71(5), 054907 (2005)

[77] Li Q., Li Z., Soff S., Bleicher M., and Stöcker H., Medium modifications of the nucleon–nucleon elastic cross section in neutron-rich intermediate energy HICs, J. Phys. G 32(4), 407 (2006)

[78] Reisdorf W., et al. (FOPI Collaboration), Systematics of pion emission in heavy ion collisions in the regime, Nucl. Phys. A 781(3–4), 459 (2007)

[79] Ferini G., Colonna M., Gaitanos T., Di Toro M., and Wolter H., Isospin effects on subthreshold Kaon production at intermediate energies, Phys. Rev. Lett. 97(20), 202301 (2006)

[80] Di Toro M., Colonna M., Ferini G., Gaitanos T., Greco V., and Wolter H. H., Heavy ion collisions at relativistic energies: Testing a nuclear matter at high baryon and isospin density, Nucl. Phys. A 782(1–4), 267 (2007)

[81] Yong G. C., Li B. A., and Chen L. W., Double neutron-proton differential transverse flow as a probe for the high density behavior of the nuclear symmetry energy, Phys. Rev. C 74(6), 064617 (2006)

[82] Leifels Y., Blaich T., Elze T. W., Emling H., Freiesleben H., et al., Exclusive studies of neutron and charged particle emission in collisions of ¹⁹⁷Au+¹⁹⁷Au at 400 MeV/nucleon, Phys. Rev. Lett. 71(7), 963 (1993)

[83] Lambrecht D., Blaich T., Elze T. W., Emling H., Freiesleben H., et al., Energy dependence of collective flow of neutrons and protons in ¹⁹⁷Au+¹⁹⁷Au collisions, Z. Phys. A 350(2), 115 (1994)

[84] Petersen H., Li Q., Zhu X., and Bleicher M., Directed and elliptic flow in heavy-ion collisions from Ebeam = 90 MeV/nucleon to Ec.m. = 200 GeV/nucleon, Phys. Rev. C 74(6), 064908 (2006)

[85] Yang J. C., Xia J. W., Xiao G. Q., Xu H. S., Zhao H. W., et al., High intensity heavy ion Accelerator Facility (HIAF) in China, Nucl. Instr. Method. B 317, 263 (2013)

[86] Xiao G. Q., Xu H. S., and Wang S. C., HIAF and CiADS national research facilities: Progress and prospect, Nucl. Phys. Rev. 34, 275 (2017)

[87] Zhou X. H., Physics opportunities at the new facility HIAF, Nucl. Phys. Rev. 35, 339 (2018)

[88] https://frib.msu.edu/

[89] Moon C.-B., Nuclear physics programs for the future rare isotope beams accelerator facility in Korea, arXiv: 1601.07271 (2016)

[90] http://www.nishina.riken.jp/RIBF/

[91] https://www.ganil-spiral2.eu/

[92] https://fair-center.eu/

[93] https://web.infn.it/epics/index.php/internalprojects/spes

[94] http://nica.jinr.ru/

[95] Liu W. P., The prospects for accelerator-basednuclear physics facilities, Physics (College Park Md.) 43, 150 (2014) (in Chinese)

[96] Adamian G. G., Antonenko N. V., Scheid W., and Volkov V. V., Treatment of competition between complete fusion and quasifission in collisions of heavy nuclei, Nucl. Phys. A 627(2), 361 (1997)

[97] Huang M. H., Gan Z. G., Zhou X. H., Li J. Q., and Scheid W., Competing fusion and quasifission reaction mechanisms in the production of superheavy nuclei, Phys. Rev. C 82(4), 044614 (2010)

[98] Adamian G. G., Antonenko N. V., and Scheid W., Characteristics of quasifission products within the dinuclear system model, Phys. Rev. C 68(3), 034601 (2003)

[99] Adamian G. G., Antonenko N. V., and Scheid W., Isotopic trends in the production of superheavy nuclei in cold fusion reactions, Phys. Rev. C 69, 011601(R) (2004)

[100] Feng Z. Q., Jin G. M., Fu F., and Li J. Q., Production cross sections of superheavy nuclei based on dinuclear system model, Nucl. Phys. A 771, 50 (2006)

[101] Adamian G. G., Antonenko N. V., and Zubov A. S., Production of unknown transactinides in asymmetry-exit-channel quasifission reactions, Phys. Rev. C 71(3), 034603 (2005)

[102] Li Q., Zuo W., Li W., Wang N., Zhao E., Li J., and Scheid W., Deformation and orientation effects in the driving potential of the dinuclear model, Eur. Phys. J. A 24(2), 223 (2005)

[103] Adamian G. G., Antonenko N. V., and Lacroix D., Production of neutron-rich Ca, Sn, and Xe isotopes in transfer-type reactions with radioactive beams, Phys. Rev. C 82(6), 064611 (2010)

[104] Mun M.-H., Adamian G. G., Antonenko N. V., Oh Y., and Kim Y., Production cross section of neutron-rich isotopes with radioactive and stable beams, Phys. Rev. C 89, 034622 (2014)

[105] Zhu L., Feng Z. Q., and Zhang F. S., Production of heavy neutron-rich nuclei in transfer reactions within the dinuclear system model, J. Phys. G 42(8), 085102 (2015)

[106] Zhu L., Su J., Xie W. J., and Zhang F. S., Production of neutron-rich transcalifornium nuclei in ²³⁸U-induced transfer reactions, Phys. Rev. C 94(5), 054606 (2016)

[107] Feng Z. Q., Production of neutron-rich isotopes around N=126 in multinucleon transfer reactions, Phys. Rev. C 95(2), 024615 (2017)

[108] Zagrebaev V. and Greiner W., New way for the production of heavy neutron-rich nuclei, J. Phys. G 35(12), 125103 (2008)

[109] Zagrebaev V. I. and Greiner W., New ideas on the production of heavy and superheavy neutron rich nuclei, Nucl. Phys. A 834(1–4), 366c (2010)

[110] Zagrebaev V. I. and Greiner W., Production of heavy trans-target nuclei in multinucleon transfer reactions, Phys. Rev. C 87(3), 034608 (2013)

[111] Zagrebaev V. and Greiner W., Production of new heavy isotopes in low-energy multinucleon transfer reactions, Phys. Rev. Lett. 101(12), 122701 (2008)

[112] Zagrebaev V. I., Fornal B., Leoni S., and Greiner W., Formation of light exotic nuclei in low-energy multinucleon transfer reactions, Phys. Rev. C 89(5), 054608 (2014)

[113] Aritomo Y., Wada T., Ohta M., and Abe Y., Fluctuation-dissipation model for synthesis of superheavy elements, Phys. Rev. C 59(2), 796 (1999)

[114] Aritomo Y., Wada T., Ohta M., and Abe Y., Diffusion mechanism for synthesis of superheavy elements, Phys. Rev. C 55, R1011(R) (1997)

[115] Shen C., Kosenko G., and Abe Y., Two-step model of fusion for the synthesis of superheavy elements, Phys. Rev. C 66, 061602(R) (2002)

[116] Feng Z. Q., Jin G. M., and Li J. Q., Production of heavy isotopes in transfer reactions by collisions of ²³⁸U+²³⁸U, Phys. Rev. C 80(6), 067601 (2009)

[117] Bao X. J., Guo S. Q., Zhang H. F., and Li J. Q., Influence of entrance channel on production cross sections of superheavy nuclei, Phys. Rev. C 96(2), 024610 (2017)

[118] Tõke J., Agnihotri D. K., Baldwin S. P., Djerroud B., Lott B., Quednau B. M., Skulski W., Schröder W. U., Sobotka L. G., Charity R. J., Sarantites D. G., and de Souza R. T., Dynamical fragment production as a mode of energy dissipation in heavy-ion reactions, Phys. Rev. Lett. 77(17), 3514 (1996)

[119] Bondorf J. P., Donangelo R., Mishustin I. N., Pethick C. J., Schulz H., and Sneppen K., Statistical multifragmentation of nuclei, Nucl. Phys. A 443(2), 321 (1985)

[120] Bondorf J., Donangelo R., Mishustin I. N., and Schulz H., Statistical multifragmentation of nuclei, Nucl. Phys. A 444(3), 460 (1985)

[121] Barz H. W., Bondorf J. P., Donangelo R., and Schulz H., Connection between the thermodynamical and the percolation models of nuclear fragmentation, Phys. Lett. B 169(4), 318 (1986)

[122] Sa B. H., Zheng Y. M., and Zhang X. Z., Disassembly of hot nuclei and relevant phase transition, Int. J. Mod. Phys. A 05(05), 843 (1990)

[123] Bondorf J. P., Botvina A. S., Iljinov A. S., Mishustin I. N., and Sneppen K., Statistical multifragmentation of nuclei, Phys. Rep. 257(3), 133 (1995)

[124] Botvina A. S. and Mishustin I. N., Statistical evolution of isotope composition of nuclear fragments, Phys. Rev. C 63(6), 061601 (2001)

[125] Buyukcizmeci N., Ogul R., and Botvina A. S., Isospin and symmetry energy effects on nuclear fragment production in liquid-gas-type phase transition region, Eur. Phys. J. A 25(1), 57 (2005)

[126] Napolitani P. and Colonna M., Bifurcations in Boltzmann–Langevin one body dynamics for fermionic systems, Phys. Lett. B 726(1–3), 382 (2013)

[127] Napolitani P. and Colonna M., Frustrated fragmentation and re-aggregation in nuclei: A non-equilibrium description in spallation, Phys. Rev. C 92(3), 034607 (2015)

[128] Gaitanos T., Larionov A. B., Lenske H., and Mosel U., Breathing mode in an improved transport approach, Phys. Rev. C 81(5), 054316 (2010)

[129] Larionov A. B., Gaitanos T., and Mosel U., Kaon and hyperon production in antiproton-induced reactions on nuclei, Phys. Rev. C 85(2), 024614 (2012)

[130] Buss O., Gaitanos T., Gallmeister K., van Hees H., Kaskulov M., Lalakulich O., Larionov A. B., Leitner T., Weil J., and Mosel U., Transport-theoretical description of nuclear reactions, Phys. Rep. 512(1–2), 1 (2012); see also https://gibuu.hepforge.org

[131] Zhang F. S. and Suraud E., Analysis of multifragmentation in a Boltzmann–Langevin approach, Phys. Rev. C 51(6), 3201 (1995)

[132] Xie W.-J., Su J., Zhu L., and Zhang F.-S., Neutron-proton effective mass splitting in a Boltzmann–Langevin approach, Phys. Rev. C 88, 061601(R) (2013)

[133] Danielewicz P., Determination of the mean-field momentum-dependence using elliptic flow, Nucl. Phys. A 673(1–4), 375 (2000)

[134] Danielewicz P. and Bertsch G. F., Production of deuterons and pions in a transport model of energetic heavy-ion reactions, Nucl. Phys. A 533(4), 712 (1991)

[135] Fuchs C. and Wolter H. H., The relativistic Landau-Vlasov method in heavy-ion collisions, Nucl. Phys. A 589(4), 732 (1995)

[136] Gaitanos T., Di Toro M., Typel S., Baran V., Fuchs C., Greco V., and Wolter H. H., On the Lorentz structure of the symmetry energy, Nucl. Phys. A 732, 24 (2004)

[137] Ko C. M. and Li Q., Relativistic Vlasov–Uehling-Uhlenbeck model for heavy-ion collisions, Phys. Rev. C 37(5), 2270 (1988)

[138] Ko C. M. and Li Q., Medium effects in high energy heavy-ion collisions, J. Phys. G 22(12), 1673 (1996)

[139] Song T. and Ko C. M., Modifications of the pion-production threshold in the nuclear medium in heavy ion collisions and the nuclear symmetry energy, Phys. Rev. C 91(1), 014901 (2015)

[140] Colonna M., Di Toro M., Guarnera A., Maccarone S., Zielinska-Pfabé M., and Wolter H. H., Fluctuations and dynamical instabilities in heavy-ion reactions, Nucl. Phys. A 642(3–4), 449 (1998)

[141] Guarnera A., Colonna M., and Chomaz Ph., 3D stochastic mean-field simulations of the spinodal fragmentation of dilute nuclei, Phys. Lett. B 373(4), 267 (1996)

[142] Colonna M., Fluctuations and symmetry energy in nuclear fragmentation dynamics, Phys. Rev. Lett. 110(4), 042701 (2013)

[143] Ono A., Horiuchi H., Maruyama T., and Ohnishi A., Antisymmetrized version of molecular dynamics with two-nucleon collisions and its application to heavy ion reactions, Prog. Theor. Phys. 87(5), 1185 (1992)

[144] Ono A., Horiuchi H., Maruyama T., and Ohnishi A., Momentum distribution of fragments in heavy-ion reactions: Dependence on the stochastic collision process, Phys. Rev. C 47(6), 2652 (1993)

[145] Su J., Zhang F. S., and Bian B. A., Odd-even effect in heavy-ion collisions at intermediate energies, Phys. Rev. C 83(1), 014608 (2011)

[146] Su J. and Zhang F. S., Non-equilibrium and residual memory in momentum space of fragmenting sources in central heavy-ion collisions, Phys. Rev. C 87(1), 017602 (2013)

[147] Su J., Cherevko K., Xie W. J., and Zhang F. S., Nonisotropic and nonsingle explosion in central ¹²⁹Xe+¹²⁰Sn collisions at 50–125 MeV/nucleon, Phys. Rev. C 89(1), 014619 (2014)

[148] Aichelin J. and Bertsch G., Numerical simulation of medium energy heavy ion reactions, Phys. Rev. C 31(5), 1730 (1985)

[149] Hartnack C., Zhuxia L., Neise L., Peilert G., Rosenhauer A., Sorge H., Aichelin J., Stöcker H., and Greiner W., Quantum molecular dynamics a microscopic model from UNILAC to CERN energies, Nucl. Phys. A 495(1–2), 303 (1989)

[150] Aichelin J., “Quantum” molecular dynamics — a dynamical microscopic n-body approach to investigate fragment formation and the nuclear equation of state in heavy ion collisions, Phys. Rep. 202(5–6), 233 (1991)

[151] Hartnack C., Puri R. K., Aichelin J., Konopka J., Bass S. A., Stöcker H., and Greiner W., Modelling the many-body dynamics of heavy ion collisions: Present status and future perspective, Eur. Phys. J. A 1(2), 151 (1998)

[152] Papa M., Maruyama T., and Bonasera A., Constrained molecular dynamics approach to fermionic systems, Phys. Rev. C 64(2), 024612 (2001)

[153] Papa M., Many-body correlations in semiclassical molecular dynamics and Skyrme interaction, Phys. Rev. C 87(1), 014001 (2013)

[154] Papa M., Giuliani G., and Bonasera A., Constrained molecular dynamics II: An N-body approach to nuclear systems, J. Comput. Phys. 208(2), 403 (2005)

[155] Wang N., Li Z., and Wu X., Improved quantum molecular dynamics model and its applications to fusion reaction near barrier, Phys. Rev. C 65(6), 064608 (2002)

[156] Zhang Y. X. and Li Z. X., Elliptic flow and system size dependence of transition energies at intermediate energies, Phys. Rev. C 74(1), 014602 (2006)

[157] Zhang Y. X., Li Z. X., and Danielewicz P., In-medium NN cross sections determined from the nuclear stopping and collective flow in heavy-ion collisions at intermediate energies, Phys. Rev. C 75(3), 034615 (2007)

[158] Zhang Y. X., Danielewicz P., Famiano M., Li Z., Lynch W. G., and Tsang M. B., The influence of cluster emission and the symmetry energy on neutron–proton spectral double ratios, Phys. Lett. B 664(1–2), 145 (2008)

[159] Zhang Y. X., Li Z. X., Zhou C. S., and Tsang M. B., Effect of isospin-dependent cluster recognition on the observables in heavy ion collisions, Phys. Rev. C 85, 051602(R) (2012)

[160] Zhang Y., Coupland D. D. S., Danielewicz P., Li Z., Liu H., Lu F., Lynch W. G., and Tsang M. B., Influence of in-medium NN cross sections, symmetry potential, and impact parameter on isospin observables, Phys. Rev. C 85(2), 024602 (2012)

[161] Zhang Y., Tsang M., and Li Z., Covariance analysis of symmetry energy observables from heavy ion collision, Phys. Lett. B 749, 262 (2015)

[162] Wang N., Li Z., Wu X., Tian J., Zhang Y., and Liu M., Further development of the improved quantum molecular dynamics model and its application to fusion reactions near the barrier, Phys. Rev. C 69(3), 034608 (2004)

[163] Feng Z. Q., Momentum dependence of the symmetry potential and its influence on nuclear reactions, Phys. Rev. C 84(2), 024610 (2011)

[164] Feng Z. Q., Nuclear in-medium effects and collective flows in heavy-ion collisions at intermediate energies, Phys. Rev. C 85(1), 014604 (2012)

[165] Cao X. G., Zhang G. Q., Cai X. Z., Ma Y. G., Guo W., Chen J. G., Tian W. D., Fang D. Q., and Wang H. W., Roles of deformation and orientation in heavy-ion collisions induced by light deformed nuclei at intermediate energy, Phys. Rev. C 81(6), 061603 (2010)

[166] Zhang G. Q., Ma Y. G., Cao X. G., Zhou C. L., Cai X. Z., Fang D. Q., Tian W. D., and Wang H. W., Unified description of nuclear stopping in central heavy-ion collisions from 10A MeV to 1.2A GeV, Phys. Rev. C 84, 034612 (2011)

[167] He W. B., Ma Y. G., Cao X. G., Cai X. Z., and Zhang G. Q., Giant dipole resonance as a fingerprint of α clustering configurations in ¹²C and ¹⁶O, Phys. Rev. Lett. 113(3), 032506 (2014)

[168] Khoa D. T., Ohtsuka N., Matin M. A., Faessler A., Huang S. W., Lehmann E., and Puri R. K., In-medium effects in the description of heavy-ion collisions with realistic NN interactions, Nucl. Phys. A 548(1), 102 (1992)

[169] Uma Maheswari V., Fuchs C., Faessler A., Sehn L., Kosov D. S., and Wang Z., In-medium dependence and Coulomb effects of the pion production in heavy ion collisions, Nucl. Phys. A 628(4), 669 (1998)

[170] Shekhter K., Fuchs C., Faessler A., Krivoruchenko M., and Martemyanov B., Dilepton production in heavy-ion collisions at intermediate energies, Phys. Rev. C 68(1), 014904 (2003)

[171] Cozma M. D., Leifels Y., Trautmann W., Li Q., and Russotto P., Toward a model-independent constraint of the high-density dependence of the symmetry energy, Phys. Rev. C 88(4), 044912 (2013)

[172] Li Q., Shen C., Guo C., Wang Y., Li Z., Lukasik J., and Trautmann W., Nonequilibrium dynamics in heavy-ion collisions at low energies available at the GSI Schwerionen Synchrotron, Phys. Rev. C 83(4), 044617 (2011)

[173] Bass S. A., Belkacem M., Bleicher M., et al., Microscopic models for ultrarelativistic heavy ion collisions, Prog. Part. Nucl. Phys. 41, 255 (1998)

[174] Sorge H., Stocker H., and Greiner W., Poincaré invariant Hamiltonian dynamics: Modelling multi-hadronic interactions in a phase space approach, Ann. Phys. 192(2), 266 (1989)

[175] Wong C.-Y., Dynamics of nuclear fluid (VIII): Time-dependent Hartree–Fock approximation from a classical point of view, Phys. Rev. C 25, 1460 (1982)

[176] Chomaz Ph., Burgio G. F., and Randrup J., Inclusion of fluctuations in nuclear dynamics, Phys. Lett. B 254(3–4), 340 (1991)

[177] Chapelle F., Burgio G. F., Chomaz Ph., and Randrup J., Fluctuations in nuclear dynamics: Comparison of different methods, Nucl. Phys. A 540(1–2), 227 (1992)

[178] Zhang F. S. and Suraud E., Boltzmann-Langevin equation, dynamical instability and multifragmentation, Phys. Lett. B 319(1–3), 35 (1993)

[179] Abe Y., Ayik S., Reinhard P. G., and Suraud E., On stochastic approaches of nuclear dynamics, Phys. Rep. 275(2–3), 49 (1996)

[180] Guarnera A., Colonna M., and Chomaz Ph., 3D stochastic mean-field simulations of the spinodal fragmentation of dilute nuclei, Phys. Lett. B 373(4), 267 (1996)

[181] Colonna M., Fluctuations and symmetry energy in nuclear fragmentation dynamics, Phys. Rev. Lett. 10(4), 042701 (2013)

[182] Gavin S., Moschelli G., and Zin C., Boltzmann–Langevin approach to pre-equilibrium correlations in nuclear collisions, Phys. Rev. C 95(6), 064901 (2017)

[183] Hohenberg P. and Kohn W., Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964)

[184] Kohn W. and Sham L. J., Self-consistent equations including exchange and correlation effects, Phys. Rev. A 140, 1133 (1965)

[185] Dreizler R. M. and Gross E. K. U., Density Functional Theory: An Approach to the Quantum Many-Body Problem, Springer-Verlag, Berlin, 1990

[186] Parr R. G. and Yang W., Density-Functional Theory of Atoms and Molecules, Clarendon, Oxford, 1989

[187] Kohn W., Electronic structure of matter — wave functions and density functionals, Rev. Mod. Phys. 71(5), 1253 (1999)

[188] Wigner E. P., On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40(5), 749 (1932)

[189] Carruthers P. and Zachariasen F., Relativistic quantum transport theory approach to multiparticle production, Phys. Rev. D 13(4), 950 (1976)

[190] Kadanoff L. P. and Baym G., Quantum Statistical Mechanics, Benjamin, New York, 1962

[191] Martin P. C. and Schwinger J., Theory of many-particle systems (I), Phys. Rev. 115(6), 1342 (1959)

[192] Schwinger J., Brownian motion of a quantum oscillator, J. Math. Phys. (N.Y.) 2(3), 407 (1961)

[193] Danielewicz P., Quantum theory of nonequilibrium processes (I), Ann. Phys. 152(2), 239 (1984)

[194] Mao G. J., Li Z. X., Zhuo Y. Z., and Han Y. L., Medium effects on the NN inelastic cross section in relativistic heavy-ion collisions, Phys. Lett. B 327(3-4), 183 (1994)

[195] Mao G., Li Z., Zhuo Y., and Yu Z., Medium effects on the NN inelastic cross section in relativistic heavy-ion collisions, Phys. Lett. B 327(3–4), 183 (1994)

[196] Han Y., Mao G., Li Z., and Zhuo Y., Effective nucleon-nucleon cross sections based on Skyrme interactions, Phys. Rev. C 50(2), 961 (1994)

[197] Mao G., Li Z., and Zhuo Y., Self-consistent relativistic Boltzmann–Uehling–Uhlenbeck equation for the Δ distribution function, Phys. Rev. C 53(6), 2933 (1996)

[198] Mao G., Neise L., Stoecker H., Greiner W., and Li Z., Relativistic transport theory of N, Δ, and N* (1440) interacting through σ, ω, and π mesons, Phys. Rev. C 57(4), 1938 (1998)

[199] Li Q., Li Z., and Mao G., Isospin dependence of nucleon-nucleon elastic cross section, Phys. Rev. C 62(1), 014606 (2000)

[200] Li Q., Li Z., and Zhao E., Density and temperature dependence of nucleon-nucleon elastic cross section, Phys. Rev. C 69(1), 017601 (2004)

[201] Cassing W. and Juchem S., Semiclassical transport of particles with dynamical spectral functions, Nucl. Phys. A 665(3–4), 377 (2000)

[202] Bertsch G. F. and Das Gupta S., A guide to microscopic models for intermediate energy heavy ion collisions, Phys. Rep. 160(4), 189 (1988)

[203] Aichelin J., Rosenhauer A., Peilert G., Stoecker H., and Greiner W., Importance of momentum-dependent interactions for the extraction of the nuclear equation of state from high-energy heavy-ion collisions, Phys. Rev. Lett. 58(19), 1926 (1987)

[204] Hartnack C. and Aichelin J., New parametrization of the optical potential, Phys. Rev. C 49(5), 2801 (1994)

[205] Gale C., Bertsch G., and Das Gupta S., Heavy-ion collision theory with momentum-dependent interactions, Phys. Rev. C 35(5), 1666 (1987)

[206] Isse M., Ohnishi A., Otuka N., Sahu P. K., and Nara Y., Mean-field effects on collective flow in high-energy heavy-ion collisions at 2-158A GeV energies, Phys. Rev. C 72(6), 064908 (2005)

[207] Hama S., Clark B. C., Cooper E. D., Sherif H. S., and Mercer R. L., Global Dirac optical potentials for elastic proton scattering from heavy nuclei, Phys. Rev. C 41(6), 2737 (1990)

[208] Cugnon J., Monte Carlo calculation of high-energy heavy-ion interactions, Phys. Rev. C 22(5), 1885 (1980)

[209] Molitoris J. J., Hoffer J. B., Kruse H., and Stöcker H., Microscopic calculations of collective flow probing the short-range nature of the nuclear force, Phys. Rev. Lett. 53(9), 899 (1984)

[210] Kruse H., Jacak B. V., Molitoris J. J., Westfall G. D., and Stöcker H., Vlasov–Uehling–Uhlenbeck theory of medium energy heavy ion reactions: Role of mean field dynamics and two body collisions, Phys. Rev. C 31(5), 1770 (1985)

[211] Bertsch G. F., Kruse H., and Gupta S. D., Boltzmann equation for heavy ion collisions, Phys. Rev. C 29(2), 673 (1984)

[212] Grégoire C., Remaud B., Sebille F., Vinet L., and Raffray Y., Semi-classical dynamics of heavy-ion reactions, Nucl. Phys. A 465(2), 317 (1987)

[213] Aichelin J., Hartnack C., Bohnet A., Zhuxia L., Peilert G., Stöcker H., and Greiner W., QMD versus BUU/VUU: Same results from different theories, Phys. Lett. B 224(1–2), 34 (1989)

[214] Feldmeier H., Fermionic molecular dynamics, Nucl. Phys. A 515(1), 147 (1990)

[215] Feldmeier H. and Schnack J., Fermionic molecular dynamics SCV252SCV133 V2, Prog. Part. Nucl. Phys. 39, 393 (1997)

[216] Maruyama T., Niita K., and Iwamoto A., Extension of quantum molecular dynamics and its application to heavy-ion collisions, Phys. Rev. C 53(1), 297 (1996)

[217] Wilets L., Henley E. M., Kraft M., and Mackellar A. D., Classical many-body model for heavy-ion collisions incorporating the Pauli principle, Nucl. Phys. A 282(2), 341 (1977)

[218] Wilets L., Yariv Y., and Chestnut R., Classical many-body model for heavy-ion collisions (II), Nucl. Phys. A 301(2), 359 (1978)

[219] Dorso J. C. and Randrup J., Classical simulation of nuclear systems, Phys. Lett. B 215(4), 611 (1988)

[220] Li Z. X., Hartnack C., Stoeker H., and Greiner W., Transition from binary processes to multifragmentation in quantum molecular dynamics for intermediate energy heavy ion collisions, Phys. Rev. C 44, 824 (1990)

[221] Wang N., Ou L., Zhang Y., and Li Z., Microscopic dynamics simulations of heavy-ion fusion reactions induced by neutron-rich nuclei, Phys. Rev. C 89(6), 064601 (2014)

[222] Yao H. and Wang N., Microscopic dynamics simulations of multinucleon transfer in ⁸⁶Kr+⁶⁴Ni at 25 MeV/nucleon, Phys. Rev. C 95(1), 014607 (2017)

[223] Jiang Y. Y., Wang N., Li Z. X., and Scheid W., Dynamical nucleus–nucleus potential at short distances, Phys. Rev. C 81(4), 044602 (2010)

[224] Wang N., Zhao K., and Li Z. X., Systematic study of ¹⁶O-induced fusion with the improved quantum molecular dynamics model, Phys. Rev. C 90(5), 054610 (2014)

[225] Tian J., Wu X., Zhao K., Zhang Y., and Li Z., Properties of the composite systems formed in the reactions of ²³⁸U+²³⁸U and ²³²Th+²⁵⁰Cf, Phys. Rev. C 77(6), 064603 (2008)

[226] Zhao K., Wu X., and Li Z., Quantum molecular dynamics study of the mass distribution of products in 7.0A MeV ²³⁸U+²³⁸U collisions, Phys. Rev. C 80(5), 054607 (2009)

[227] Zhang Y. X., Tsang M. B., Li Z. X., and Liu H., Constraints on nucleon effective mass splitting with heavy ion collisions, Phys. Lett. B 732, 186 (2014)

[228] Cugnon J. J., L’Hôte D., and Vandermeulen J., Simple parametrization of cross-sections for nuclear transport studies up to the GeV range, Nucl. Instrum. Methods B 111(3–4), 215 (1996)

[229] Ono A., Xu J., Colonna M., Danielewicz P., Ko C. M., et al., Comparison of heavy-ion transport simulations: Collision integral with pions and Δ resonances in a box, Phys. Rev. C 100(4), 044617 (2019)

[230] Wang N. and Li T., Shell and isospin effects in nuclear charge radii, Phys. Rev. C 88, 011301(R) (2013)

[231] Trzcińska A., Jastrzebski J., Lubiński P., Hartmann F. J., Schmidt R., von Egidy T., and Kłos B., Neutron density distributions deduced from antiprotonic atoms, Phys. Rev. Lett. 87(8), 082501 (2001)

[232] Bartel J., Quentin P., Brack M., Guet C., and Håkansson H. B., Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force, Nucl. Phys. A 386(1), 79 (1982)

[233] Wen K., Sakata F., Li Z. X., Wu X. Z., Zhang Y. X., and Zhou S. G., Non-Gaussian fluctuations and non-Markovian effects in the nuclear fusion process: Langevin dynamics emerging from quantum molecular dynamics simulations, Phys. Rev. Lett. 111(1), 012501 (2013)

[234] Nebauer R., Aichelin J., Assenard M., Auger G., Bacri C. O., et al., Multifragmentation in Xe(50A MeV) + Sn: Confrontation of theory and data, Nucl. Phys. A 658(1), 67 (1999)

[235] Liu T. X., van Goethem M. J., Liu X. D., Lynch W. G., Shomin R., et al., Isotope yields from central ^112,124Sn+^112,124Sn collisions: Dynamical emission? Phys. Rev. C 69(1), 014603 (2004)

[236] Russotto P., De Filippo E., Pagano A., Piasecki E., Amorini F., et al., Strong enhancement of dynamical emission of heavy fragments in the neutron-rich ¹²⁴Sn+⁶⁴Ni reaction at 35A MeV, Phys. Rev. C 81(6), 064605 (2010)

[237] Kohley Z., May L. W., Wuenschel S., Colonna M., Di Toro M., et al., Transverse collective flow and midrapidity emission of isotopically identified light charged particles, Phys. Rev. C 83(4), 044601 (2011)

[238] Hudan S., Chbihi A., Frankland J. D., Mignon A., Wieleczko J. P., et al., Characteristics of the fragments produced in central collisions of ¹²⁹Xe+^natSn from 32A to 50A MeV, Phys. Rev. C 67(6), 064613 (2003)

[239] Dorso C. O. and Randrup J., Early recognition of clusters in molecular dynamics, Phys. Lett. B 301(4), 328 (1993)

[240] Puri R. K. and Aichelin J., Simulated annealing clusterization algorithm for studying the multifragmentation, J. Comput. Phys. 162(1), 245 (2000)

[241] Gossiaux P. B., Puri R. K., Hartnack C., and Aichelin J., The multifragmentation of spectator matter, Nucl. Phys. A 619(3–4), 379 (1997)

[242] Goyal S. and Puri R. K., Formation of fragments in heavy-ion collisions using a modified clusterization method, Phys. Rev. C 83(4), 047601 (2011)

[243] Bleicher M., Zabrodin E., Spieles C., Bass S. A., Ernst C., Soff S., Bravina L., Belkacem M., Weber H., Stöcker H., and Greiner W., Relativistic hadron–hadron collisions in the ultra-relativistic quantum molecular dynamics model, J. Phys. G 25(9), 1859 (1999)

[244] Li Q., Graf G., and Bleicher M., Ultrarelativistic quantum molecular dynamics calculations of two-pion Hanbury–Brown–Twiss correlations in central Pb–Pb collisions at sNN=2.76 TeV, Phys. Rev. C 85(3), 034908 (2012)

[245] https://urqmd.org/

[246] Li Q., Bleicher M., and Stöcker H., The effect of “pre-formed” hadron potentials on the dynamics of heavy ion collisions and the HBT puzzle, Phys. Lett. B 659(3), 525 (2008)

[247] Petersen H., Steinheimer J., Burau G., Bleicher M., and Stocker H., Fully integrated transport approach to heavy ion reactions with an intermediate hydrodynamic stage, Phys. Rev. C 78(4), 044901 (2008)

[248] Li Q., Steinheimer J., Petersen H., Bleicher M., and Stocker H., Effects of a phase transition on HBT correlations in an integrated Boltzmann+hydrodynamics approach, Phys. Lett. B 674(2), 111 (2009)

[249] Russotto P., Gannon S., Kupny S., Lasko P., Acosta L., et al., Results of the ASY-EOS experiment at GSI: The symmetry energy at suprasaturation density, Phys. Rev. C 94(3), 034608 (2016)

[250] Guo C., Wang Y., Li Q., and Zhang F. S., Effect of the spin–orbit interaction on flows in heavy-ion collisions at intermediate energies, Phys. Rev. C 90(3), 034606 (2014)

[251] Sun Y., Wang Y., Li Q., and Wang F., Effect of internal magnetic field on collective flow in heavy ion collisions at intermediate energies, Phys. Rev. C 99(6), 064607 (2019)

[252] Liu Y., Wang Y., Li Q., and Liu L., Collective flows of pions in Au+Au collisions at energies 1.0 and 1.5 GeV/nucleon, Phys. Rev. C 97, 034602 (2018)

[253] Du Y., Wang Y., Li Q., and Liu L., The effect of Lorentz-like force on collective flows of K⁺ in Au+Au collisions at 1.5 GeV/nucleon, Sci. China Phys. Mech. Astron. 61(6), 062011 (2018)

[254] Wang Y., Li Q., Leifels Y., and Le Fevre A., Study of the nuclear symmetry energy from the rapidity-dependent elliptic flow in heavy-ion collisions around 1 GeV/nucleon regime, Phys. Lett. B 802, 135249 (2020)

[255] Wang Y., Guo C., Li Q., Zhang H., Li Z., and Trautmann W., Collective flow of light particles in Au+Au collisions at intermediate energies, Phys. Rev. C 89(3), 034606 (2014)

[256] Li P. C., Wang Y. J., Li Q. F., and Zhang H. F., Effects of the in-medium nucleon–nucleon cross section on collective flow and nuclear stopping in heavy-ion collisions in the Fermi-energy domain, Phys. Rev. C 97(4), 044620 (2018)

[257] Wu X., Tian J., Ning W., Kai Z., and Li Z., Microscopic study on dynamic barrier in fusion reactions, Chin. Phys. C 12, 1317 (2004)

[258] Denisov V. Yu. and Nörenberg W., Entrance channel potentials in the synthesis of the heaviest nuclei, Eur. Phys. J. A 15(3), 375 (2002)

[259] Zanganeh V., Wang N., and Ghodsi O. N., Dynamical nucleus–nucleus potential and incompressibility of nuclear matter, Phys. Rev. C 85(3), 034601 (2012)

[260] Aguilera E. F., Kolata J. J., and Tighe R. J., Nuclear structure effects in the sub-barrier fusion of ¹⁶O+^{70,72,73,74,76}Ge, Phys. Rev. C 52(6), 3103 (1995)

[261] Kurtukian-Nieto T., Benlliure J., Schmidt K. H., Audouin L., Becker F., et al., Production cross sections of heavy neutron-rich nuclei approaching the nucleosynthesis r-process path around A=195, Phys. Rev. C 89(2), 024616 (2014)

[262] Hildenbrand K. D., Freiesleben H., Pühlhofer F., Schneider W. F. W., Bock R., Harrach D., and Specht H. J., Reaction between ²³⁸U and ²³⁸U at 7.42 MeV/nucleon, Phys. Rev. Lett. 39(17), 1065 (1977)

[263] Schädel M., Kratz J. V., Ahrens H., Brüchle W., Franz G., Gäggeler H., Warnecke I., Wirth G., Herrmann G., Trautmann N., and Weis M., Isotope distributions in the reaction of ²³⁸U with ²³⁸U, Phys. Rev. Lett. 41(7), 469 (1978)

[264] Essel H., Hartel K., Henning W., Kienle P., et al., Charge and mass transfer in the reaction ¹³⁶Xe+²⁰⁸Pb at energies close to the coulomb barrier, Z. Phys. A At. Nucl. 289(3), 265 (1979)

[265] Freiesleben H., Hildenbrand K. D., Pühlhofer F., Schneider W. F. W., Bock R., Harrach D., and Specht H. J., The reaction ²³⁸U+²³⁸U at 7.42 MeV/u, Z. Phys. A At. Nucl. 292(2), 171 (1979)

[266] Schädel M., Brüchle W., Gäggeler H., Kratz J. V., Sümmerer K., et al., Actinide production in collisions of ²³⁸U with ²⁴⁸Cm, Phys. Rev. Lett. 48(13), 852 (1982)

[267] Moody K. J., Lee D., Welch R. B., Gregorich K. E., Seaborg G. T., Lougheed R. W., and Hulet E. K., Actinide production in reactions of heavy ions with ²⁴⁸Cm, Phys. Rev. C 33(4), 1315 (1986)

[268] Kratz J. V., Schädel M., and Gäggeler H. W., Reexamining the heavy-ion reactions ²³⁸U+²³⁸U and ²³⁸U+²⁴⁸Cm and actinide production close to the barrier, Phys. Rev. C 88(5), 054615 (2013)

[269] Golabek C., Villari A. C. C., Heinz S., Mittig W., Bhattacharyya S., et al., Search for a long living giant system in ²³⁸U+²³⁸U collisions near the coulomb barrier, Int. J. Mod. Phys. E 17(10), 2235 (2008)

[270] Mijatović T., Szilner S., Corradi L., Montanari D., Pollarolo G., et al., Multinucleon transfer reactions in the ⁴⁰Ar+²⁰⁸Pb system, Phys. Rev. C 94, 064616 (2016)

[271] Barrett J. S., Loveland W., Yanez R., Zhu S., Ayangeakaa A. D., et al., ¹³⁶e+²⁰⁸Pb reaction: A test of models of multinucleon transfer reactions, Phys. Rev. C 91(6), 064615 (2015)

[272] Loveland W., Synthesis of transactinide nuclei using radioactive beams, Phys. Rev. C 76(1), 014612 (2007)

[273] Yanez R. and Loveland W., Predicting the production of neutron-rich heavy nuclei in multinucleon transfer reactions using a semi-classical model including evaporation and fission competition, GRAZING-F, Phys. Rev. C 91(4), 044608 (2015)

[274] Golabek C. and Simenel C., Collision dynamics of two ²³⁸U atomic nuclei, Phys. Rev. Lett. 103(4), 042701 (2009)

[275] Kedziora D. J. and Simenel C., New inverse quasifission mechanism to produce neutron-rich transfermium nuclei, Phys. Rev. C 81(4), 044613 (2010)

[276] Wang N., Li Z., Wu X., and Zhao E., Search for possible way of producing super-heavy elements: Dynamic study on damped reactions of ²⁴⁴Pu+²⁴⁴Pu, ²³⁸U+²³⁸U and ¹⁹⁷Au+¹⁹⁷Au, Mod. Phys. Lett. A 20(34), 2619 (2005)

[277] Zhao K., Li Z., Wu X., and Zhang Y., Production probability of superheavy fragments at various initial deformations and orientations in the ²³⁸U+²³⁸U reaction, Phys. Rev. C 88(4), 044605 (2013)

[278] Zhao K., Li Z., Wang N., Zhang Y., Li Q., Wang Y., and Wu X., Production mechanism of neutron-rich transuranium nuclei in ²³⁸U+²³⁸U collisions at near-barrier energies, Phys. Rev. C 92(2), 024613 (2015)

[279] Zhao K., Li Z., Zhang Y., Wang N., Li Q., Shen C., Wang Y., and Wu X., Production of unknown neutron-rich isotopes in ²³⁸U+²³⁸U collisions at near-barrier energy, Phys. Rev. C 94(2), 024601 (2016)

[280] Ayik S., Yilmaz B., Yilmaz O., Umar A. S., and Turan G., Multinucleon transfer in central collisions of ²³⁸U+²³⁸U, Phys. Rev. C 96(2), 024611 (2017)

[281] Sekizawa K., and Yabana K., Time-dependent Hartree-Fock calculations for multinucleon transfer processes in ^40,48Ca+¹²⁴Sn, ⁴⁰Ca+²⁰⁸Pb, and ⁵⁸Ni+²⁰⁸Pb reactions, Phys. Rev. C 88(1), 014614 (2013)

[282] Ghiorso A., Lee D., Somerville L. P., Loveland W., Nitschke J. M., et al., Evidence for the synthesis of 267110 produced by the ⁵⁹Co+²⁰⁹Bi reaction, Nucl. Phys. A 583, 861 (1995)

[283] Lazarev Yu. A., Lobanov Y. V., Oganessian Y. T., Utyonkov V. K., Abdullin F. S., et al., α decay of ²⁷³110: Shell closure at N=162, Phys. Rev. C 54(2), 620 (1996)

[284] Ninov V., Gregorich K. E., Loveland W., Ghiorso A., Hoffman D. C., et al., Observation of superheavy nuclei produced in the reaction of ⁸⁶Kr with ²⁰⁸Pb, Phys. Rev. Lett. 83(6), 1104 (1999)

[285] Oganessian Yu. Ts., Utyonkov V. K., Lobanov Yu. V., Abdullin F. S., Polyakov A. N., et al., Synthesis of superheavy nuclei in the ⁴⁸Ca+²⁴⁴Pu reaction, Phys. Rev. Lett. 83(16), 3154 (1999)

[286] Oganessian Yu. Ts., et al., Observation of the decay of ²⁹²116, Phys. Rev. C 63, 011301(R) (2000)

[287] Oganessian Yu. Ts., Utyonkoy V. K., Lobanov Y. V., Abdullin F. S., Polyakov A. N., et al., Experiments on the synthesis of element 115 in the reaction ²⁴³Am (⁴⁸Ca, xn)^291-x115, Phys. Rev. C 69(2), 021601 (2004)

[288] Oganessian Yu. Ts., Utyonkov V. K., Lobanov Y. V., Abdullin F. S., Polyakov A. N., et al., Measurements of cross sections and decay properties of the isotopes of elements 112, 114, and 116 produced in the fusion reactions ^233,238U, ²⁴²Pu, and ²⁴⁸Cm+⁴⁸Ca, Phys. Rev. C 70(6), 064609 (2004)

[289] Oganessian Yu. Ts., Utyonkov V. K., Lobanov Y. V., Abdullin F. S., Polyakov A. N., et al., Synthesis of the isotopes of elements 118 and 116 in the ²⁴⁹Cf and ²⁴⁵Cm+⁴⁸Ca fusion reactions, Phys. Rev. C 74(4), 044602 (2006)

[290] Oganessian Yu. Ts., Utyonkov V. K., Dmitriev S. N., Lobanov Y. V., Itkis M. G., et al., Synthesis of elements 115 and 113 in the reaction ²⁴³Am+⁴⁸Ca, Phys. Rev. C 72(3), 034611 (2005)

[291] Ellison P. A., Gregorich K. E., Berryman J. S., Bleuel D. L., Clark R. M., et al., New superheavy element isotopes: ²⁴²Pu (⁴⁸Ca, 5n)²⁸⁵114, Phys. Rev. Lett. 105(18), 182701 (2010)

[292] Oganessian Yu. Ts., Abdullin F. Sh., Dmitriev S. N., Gostic J. M., Hamilton J. H., et al., New insights into the ²⁴³Am+⁴⁸Ca reaction products previously observed in the experiments on elements 113, 115, and 117, Phys. Rev. Lett. 108(2), 022502 (2012)

[293] Oganessian Yu. Ts., Abdullin F. Sh., Alexander C., Binder J., Boll R. A., et al., Production and decay of the heaviest nuclei ^293,294117 and ²⁹⁴118, Phys. Rev. Lett. 109(16), 162501 (2012)

[294] Reisdorf W., Hessberger F. P., Hildenbrand K. D., Hofmann S., Münzenberg G., et al., Fusability and fissionability in ⁸⁶Kr-induced reactions near and below the fusion barrier, Nucl. Phys. A 444(1), 154 (1985)

[295] Charity R., McMahan M. A., Wozniak G. J., McDonald R. J., Moretto L. G., et al., Systematics of complex fragment emission in niobium-induced reactions, Nucl. Phys. A 483(2), 371 (1988)

[296] Wang N., Wu T., Zeng J., Yang Y. X., and Ou L., Improvement on fermionic properties and new isotope production in molecular dynamics simulations, J. Phys. G 43(6), 065101 (2016)

[297] Wang N. and Guo L., New neutron-rich isotope production in ¹⁵⁴Sm+¹⁶⁰Gd, Phys. Lett. B 760, 236 (2016)

[298] Audi G., Kondev F. G., Wang M., Pfeiffer B., Sun X., Blachot J., and MacCormick M., The Nubase2012 evaluation of nuclear properties, Chin. Phys. C 36(12), 1157 (2012)

[299] Skwira-Chalot I., Siwek-Wilczynska K., Wilczynski J., Amorini F., Anzalone A., et al., Dynamical evolution of the ¹⁹⁷Au+¹⁹⁷Au system at 15 MeV/nucleon, Int. J. Mod. Phys. E 15(02), 495 (2006)

[300] Skwira-Chalot I., Siwek-Wilczynska K., Wilczynski J., Amorini F., Anzalone A., et al., Ternary reactions in ¹⁹⁷Au+¹⁹⁷Au collisions revisited, Int. J. Mod. Phys. E 16(02), 511 (2007)

[301] Skwira-Chalot I., Siwek-Wilczynska K., Wilczynski J., Amorini F., Anzalone A., et al., Fast ternary and quaternary breakup of the ¹⁹⁷Au+¹⁹⁷Au system in collisions at 15 MeV/nucleon, Phys. Rev. Lett. 101(26), 262701 (2008)

[302] Wilczynski J., Skwira-Chalot I., Siwek-Wilczynska K., Wilczynski J., et al., Observation of fast collinear partitioning of the ¹⁹⁷Au+¹⁹⁷Au system into three and four fragments of comparable size, Phys. Rev. C 81(2), 024605 (2010)

[303] Tian J. L., Wu X. Z., Li Z. X., Zhao K., Zhang Y., Li X., and Yan S., Mechanism of ternary breakup in the reaction ¹⁹⁷Au+¹⁹⁷Au at 15A MeV, Phys. Rev. C 82(5), 054608 (2010)

[304] Li X., Tian J. L., Yan S. W., Cheng J. X., and Jiang X., Angular distributions of fragments produced in ternary reaction of ¹⁹⁷Au+¹⁹⁷Au at 15A MeV, Mod. Phys. Lett. A 26(06), 449 (2011)

[305] Zhang Y., Zhou C., Chen J., Wang N., Zhao K., and Li Z. X., Correlation between the fragmentation modes and light charged particles emission in heavy ion collisions, Sci. China Phys. Astro. Mech. 58(11), 112002 (2015)

[306] Natowitz J. B., Wada R., Hagel K., Keutgen T., Murray M., Makeev A., Qin L., Smith P., and Hamilton C., Caloric curves and critical behavior in nuclei, Phys. Rev. C 65(3), 034618 (2002)

[307] D’Agostino M., Bruno M., Gulminelli F., Bougault R., Cannata F., Chomaz P., Gramegna F., Le Neidre N., Moroni A., and Vannini G., Experimental signals of phase transition, Nucl. Phys. A 734, 512 (2004)

[308] Richert J. and Wagner P., Microscopic model approaches to fragmentation of nuclei and phase transitions in nuclear matter, Phys. Rep. 350(1), 1 (2001)

[309] Chomaz P., Colonna M., and Randrup J., Nuclear spinodal fragmentation, Phys. Rep. 389(5–6), 263 (2004)

[310] Poskanzer A. M., Butler G. W., and Hyde E. K., Fragment production in the interaction of 5.5 GeV protons with uranium, Phys. Rev. C 3(2), 882 (1971)

[311] Westfall G. D., Sextro R. G., Poskanzer A. M., Zebelman A. M., Butler G. W., and Hyde E. K., Energy spectra of nuclear fragments produced by high energy protons, Phys. Rev. C 17(4), l368 (1978)

[312] Pochodzalla J., Möhlenkamp T., Rubehn T., Schüttauf A., Wörner A., et al., Probing the nuclear liquid-gas phase transition, Phys. Rev. Lett. 75(6), 1040 (1995)

[313] Baldo M., Burgio G. F., and Rapisarda A., Dynamics of fragment formation in the nuclear spinodal region, Phys. Rev. C 51(1), 198 (1995)

[314] Nayak S. K., Ramaswamy R., and Chakravarty C., Maximal Lyapunov exponent in small atomic clusters, Phys. Rev. E 51(4), 3376 (1995)

[315] Zhang Y., Wu X., and Li Z., Connection between the largest Lyapunov exponent, density fluctuation, and multifragmentation in excited nuclear systems, Phys. Rev. C 69(4), 044609 (2004)

[316] Eckmann J. P. and Ruelle D., Ergodic theory of chaos and strange attractors, Rev. Mod. Phys. 57(3), 617 (1985)

[317] Gu Y., Evidences of classical and quantum chaos in the time evolution of nonequilibrium ensembles, Phys. Lett. A 149(2–3), 95 (1990)

[318] Łukasik J., Auger G., Begemann-Blaich M. L., Bellaize N., Bittiger R., et al., Directed and elliptic flow in ¹⁹⁷Au+¹⁹⁷Au at intermediate energies, Phys. Lett. B 608(3–4), 223 (2005)

[319] Pinkenburg C., et al. (E895 Collaboration), Prepared for Centennial Celebration and Meeting of the American Physical Society (Combining Annual APS General Meeting and the Joint Meeting of the APS and the AAPT), Atlanta, Georgia, Mar. 20–26 1999

[320] Chung P., et al. (E895 Collaboration), Centrality and momentum-selected elliptic flow: Tighter constraints for the nuclear equation of state, Phys. Rev. C 66(2), 021901 (2002)

[321] Alt C., et al. (NA49 Collaboration), Directed and elliptic flow of charged pions and protons in Pb+Pb collisions at 40A and 158A GeV, Phys. Rev. C 68(3), 034903 (2003)

[322] Andronic A., et al. (FOPI Collaboration), Excitation function of elliptic flow in Au+Au collisions and the nuclear matter equation of state, Phys. Lett. B 612(3–4), 173 (2005)

[323] Wang Y. J., Guo C. C., Li Q. F., Zhang H. F., Leifels Y., and Trautmann W., Constraining the high-density nuclear symmetry energy with the transverse-momentum-dependent elliptic flow, Phys. Rev. C 89(4), 044603 (2014)

[324] Wang Y. J., Guo C. C., Li Q. F., Li Z. X., Su J., and Zhang H. F., Influence of differential elastic nucleon–nucleon cross section on stopping and collective flow in heavy-ion collisions at intermediate energies, Phys. Rev. C 94(2), 024608 (2016)

[325] Li P. C., Wang Y. J., Li Q. F., and Zhang H. F., Collective flow and nuclear stopping in heavy ion collisions in Fermi energy domain, Nucl. Sci. Tech. 29(12), 177 (2018)

[326] Wang Y. J., Guo C. C., Li Q. F., Le F’evre A., Leifels Y., and Trautmann W., Determination of the nuclear incompressibility from the rapidity-dependent elliptic flow in heavy-ion collisions at beam energies 0.4A−1.0A GeV, Phys. Lett. B 778, 207 (2018)

[327] Gudowski W., Accelerator-driven transmutation projects. The importance of nuclear physics research for waste transmutation, Nucl. Phys. A 654(1–2), c436 (1999)

[328] Casolino M., Bidoli V., Morselli A., Narici L., De Pascale M. P., et al., Dual origins of light flashes seen in space, Nature 422(6933), 680 (2003)

[329] Larsson B., PhD thesis, Uppsala University, 1962

[330] Trainor J. H., Instrument and spacecraft faults associated with nuclear radiation in space, Adv. Space Res. 14(10), 685 (1994)

[331] Smith M. S. and Rehm K. E., Nuclear astrophysics measurements with radioactive beams, Annu. Rev. Nucl. Part. Sci. 51(1), 91 (2001)

[332] Bauer G. S., Proceedings of the 2nd International Conference on Accelerator Driven Transmutation Technologies, p.159, Kalmar, Uppsala University, 1996

[333] Tech. Rep. No. DOE/Er-0705, Department of Energy, 1997

[334] Carpenter J. M., Pulsed spallation neutron sources for slow neutron scattering, Nucl. Instrum. Methods 145(1), 91 (1977)

[335] Bowman C. D., Arthur E. D., Lisowski P. W., Lawrence G. P., Jensen R. J., et al., Nuclear energy generation and waste transmutation using an accelerator-driven intense thermal neutron source, Nucl. Instrum. Methods Phys. Res. Sect. A 320(1–2), 336 (1992)

[336] Takizuka T., Proceedings of the International Conference on Accelerator-driven Transmutation Technologies and Application, p. 64, AIP Press, Woodbury, NY, 1995

[337] Rejmund F., Mustapha B., Armbruster P., Benlliure J., Bernas M., et al., Measurement of isotopic cross sections of spallation residues in 800A MeV ¹⁹⁷Au+p collisions, Nucl. Phys. A 683(1–4), 540 (2001)

[338] Fernández-Dominguez B., Armbruster P., Audouin L., Benlliure J., Bernas M., et al., Nuclide cross–sections of fission fragments in the reaction ²⁰⁸Pb+p at 500A MeV, Nucl. Phys. A 747(2–4), 227 (2005)

[339] Audouin L., Tassan-Got L., Armbruster P., Benlliure J., Bernas M., et al., Evaporation residues produced in spallation of ²⁰⁸Pb by protons at 500A MeV, Nucl. Phys. A 768(1–2), 1 (2006)

[340] Benlliure J., Armbruster P., Bernas M., Boudard A., Dufour J. P., et al., Isotopic production cross sections of fission residues in ¹⁹⁷Au-on-proton collisions at 800A MeV, Nucl. Phys. A 683(1-4), 513 (2001)

[341] Iwase H., Niita K., and Nakamura T., Development of general-purpose particle and heavy ion transport Monte Carlo code, J. Nucl. Sci. Technol. 39(11), 1142 (2002)

[342] Nara Y., Otuka N., Ohnishi A., Niita K., and Chiba S., Relativistic nuclear collisions at 10A GeV energies from p+Be to Au+Au with the hadronic cascade model, Phys. Rev. C 61(2), 024901 (2000)

[343] Niita K., Chiba S., Maruyama T., Maruyama T., Takada H., Fukahori T., Nakahara Y., and Iwamoto A., Analysis of the (N,xN′) reactions by quantum molecular dynamics plus statistical decay model, Phys. Rev. C 52(5), 2620 (1995)

[344] Boudard A., Cugnon J., Leray S., and Volant C., Intranuclear cascade model for a comprehensive description of spallation reaction data, Phys. Rev. C 66(4), 044615 (2002)

[345] Ou L., Li Z., Wu X., Tian J., and Sun W., A study of proton-induced spallation reactions by the improved quantum molecular dynamics model plus statistical decay models, J. Phys. G 36(12), 125104 (2009)

[346] Ou L., Zhang Y., Tian J., and Li Z., Analysis of intermediate energy proton-induced spallation reactions by an improved quantum molecular dynamics plus statistical decay model, J. Phys. G 34(5), 827 (2007)

[347] Ou L., Zhang Y., and Li Z., Mechanism of proton-induced reactions on targets ¹⁶O, ²⁷Al, ⁵⁶Fe, ¹¹²Cd, ¹⁸⁴W and ²⁰⁸Pb at Ep=800 MeV, Chin. Phys. Lett. 24(1), 72 (2007)

[348] Wei D., Mao L., Wang N., Liu M., and Ou L., Further study on mechanism of production of light complex particles in nucleon-induced reactions, Nucl. Phys. A 933, 114 (2015)

[349] Wei D., Wang N., and Ou L., Mechanism of the production of light complex particles in nucleon-induced reactions, J. Phys. G 41(3), 035104 (2014)

[350] Meier M. M., Amian W. B., Goulding C. A., Morgan G. L., and Moss C. E., Differential neutron production cross sections for 256-MeV protons, Nucl. Sci. Eng. 110(3), 289 (1992)

[351] Trebukhovsky Y. V., Titarenko Y. E., et al., LANL Report LA-UR-03-6071; Los Alamos, 2003; ITEP Print 03-03, Moscow, 2003

[352] Trebukhovsky Y. V., Titarenko Y. E., Batyaev V. F., Mulambetov R. D., Mulambetova S. V., et al., Double-differential cross sections for the production of neutrons from Pb, W, Zr, Cu, and Al targets irradiated with 0.8-, 1.0-, and 1.6-GeV protons, Phys. At. Nucl. 68(1), 3 (2005)

[353] Chrien R. E., Krieger T. J., Sutter R. J., May M., Palevsky H., Stearns R. L., Kozlowski T., and Bauer T., Proton spectra from 800 MeV protons on selected nuclides, Phys. Rev. C 21(3), 1014 (1980)

[354] Villagrasa-Canton C., Boudard A., Ducret J. E., Fernandez B., Leray S., et al., Spallation residues in the reaction ⁵⁶Fe+p at 0.3A, 0.5A, 0.75A, 1.0A, and 1.5A GeV, Phys. Rev. C 75(4), 044603 (2007)

[355] Machner H., Aschman D. G., Baruth-Ram K., Carter J., Cowley A. A., et al., Isotopic production cross sections in proton–nucleus collisions at 200 MeV, Phys. Rev. C 73(4), 044606 (2006)

[356] Uozumi Y., Evtoukhovitch P., Fukuda H., Imamura M., Iwamoto H., et al., Magnitude factor systematics of Kalbach phenomenology for reactions emitting helium and lithium ions, Nucl. Instrum. Methods Phys. Res. A 571(3), 743 (2007)

[357] Bohnet A., Ohtsuka N., Aichelin J., Linden R., and Faessler A., Quantum molecular-dynamics approach to heavy-ion collisions with Brueckner G-matrix cross sections, Nucl. Phys. A 494(2), 349 (1989)

[358] Schulze H. J., Schnell A., Röpke G., and Lombardo U., Nucleon–nucleon cross sections in nuclear matter, Phys. Rev. C 55(6), 3006 (1997)

[359] Li G. Q. and Machleidt R., Microscopic calculation of in-medium nucleon–nucleon cross sections, Phys. Rev. C 48, 1702 (1993)

[360] Li G. Q. and Machleidt R., Microscopic calculation of in-medium proton–proton cross sections, Phys. Rev. C 49, 566 (1994)

[361] Fuchs C., Faessler A., and El-Shabshiry M., Off-shell behavior of the in-medium nucleon–nucleon cross section, Phys. Rev. C 64(2), 024003 (2001)

[362] Cai Y. H., Song H. Q., and Lombardo U., In-medium nucleon–nucleon cross section, Chin. Phys. Lett. 13(6), 420 (1996)

[363] Zhang H. F., Lombardo U., and Zuo W., Transport parameters in neutron stars from in-medium NN cross sections, Phys. Rev. C 82(1), 015805 (2010)

[364] Huber S. and Aichelin J., Production of Δ- and N*-resonances in the one-boson exchange model, Nucl. Phys. A 573(4), 587 (1994)

[365] Machleidt R., Holinde K., and Elster C., The bonn meson-exchange model for the nucleon–nucleon interaction, Phys. Rep. 149(1), 1 (1987)

[366] Larionov A. and Mosel U., The NN→NΔ cross section in nuclear matter, Nucl. Phys. A 728(1–2), 135 (2003)

[367] Mao G., Li Z., Zhuo Y., Han Y., and Yu Z., Study of in-medium NN inelastic cross section from relativistic Boltzmann–Uehling–Uhlenbeck approach, Phys. Rev. C 49(6), 3137 (1994)

[368] Mao G. J., Li Z. X., Zhuo Y. Z., and Yu Z. Q., Medium effects on the NN inelastic cross section in relativistic heavy-ion collisions, Phys. Lett. B 327(3–4), 183 (1994)

[369] Li Q., Li Z., and Mao G., Isospin dependence of nucleon-nucleon elastic cross section, Phys. Rev. C 62(1), 014606 (2000)

[370] Mao G., Relativistic Microscopic Quantum Transport Equation, NOVA, 2005

[371] Li Q. and Li Z., The isospin dependent nucleon–nucleon inelastic cross section in the nuclear medium, Phys. Lett. B 773, 557 (2017)

[372] Cui Y., Zhang Y., and Li Z., Effect of energy conservation on the in-medium NN→NΔ cross section in isospin-asymmetric nuclear matter, Phys. Rev. C 98(5), 054605 (2018)

[373] Pan Q. and Danielewicz P., From sideward flow to nuclear compressibility, Phys. Rev. Lett. 70, 2062 (1993) [Erratum Phys. Rev. Lett. 70, 3523 (1993)]

[374] Shi L. and Danielewicz P., Nuclear isospin diffusivity, Phys. Rev. C 68(6), 064604 (2003)

[375] Reisdorf W., et al. (FOPI Collaboration), Nuclear stopping from 0.09A to 1.93A GeV and its correlation to flow, Phys. Rev. Lett. 92(23), 232301 (2004)

[376] Persram D. and Gale C., Elliptic flow in intermediate energy heavy ion collisions and in-medium effects, Phys. Rev. C 65(6), 064611 (2002)

[377] Lopez O., Durand D., Lehaut G., Borderie B., Frankland J. D., et al., In-medium effects for nuclear matter in the Fermi-energy domain, Phys. Rev. C 90(6), 064602 (2014)

[378] Zhang Y. X., Wang Y. J., Colonna M., Danielewicz P., Ono A., et al., Comparison of heavy-ion transport simulations: Collision integral in a box, Phys. Rev. C 97(3), 034625 (2018)

[379] Ou L. and He X., In-medium nucleon-nucleon elastic cross-sections determined from the nucleon induced reaction cross-section data, Chin. Phys. C 43(4), 044103 (2019)

[380] Vretenar D., Niksic T., and Ring P., Relativistic nuclear energy density functionals, Int. J. Mod. Phys. E 19(04), 548 (2010)

[381] Ring P., Relativistic mean field theory in finite nuclei, Prog. Part. Nucl. Phys. 73, 193 (1996)

[382] Ring P., Covariant density functional theory and applications to finite nuclei, in: G. A. Lalazissis, P. Ring, and D. Vretenar (Eds.), Extended density functionals in nuclear structure physics, Lect. Notes Phys. 641, 175 (2004)

[383] Brockmann R., Relativistic Hartree–Fock description of nuclei, Phys. Rev. C 18(3), 1510 (1978)

[384] Horowitz C. J. and Serot B. D., Properties of nuclear and neutron matter in a relativistic Hartree–Fock theory, Nucl. Phys. A 399(2), 529 (1983)

[385] Bouyssy A., Mathiot J. F., Van Giai N., and Marcos S., Relativistic description of nuclear systems in the Hartree–Fock approximation, Phys. Rev. C 36(1), 380 (1987)

[386] Long W. H., Giai N. V., and Meng J., Density-dependent relativistic Hartree–Fock approach, Phys. Lett. B 640(4), 150 (2006)

[387] Nikšić T., Vretenar D., and Ring P., Relativistic nuclear energy density functionals: Adjusting parameters to binding energies, Phys. Rev. C 78(3), 034318 (2008)

[388] Lalazissis G. A., Nikšić T., Vretenar D., and Ring P., New relativistic mean-field interaction with density-dependent meson-nucleon couplings, Phys. Rev. C 71(2), 024312 (2005)

[389] Zhao P. W., Li Z. P., Yao J. M., and Meng J., New parametrization for the nuclear covariant energy density functional with a point-coupling interaction, Phys. Rev. C 82(5), 054319 (2010)

[390] Skyrme T. H. R., CVII. The nuclear surface, Philos. Mag. 1(11), 1043 (1956)

[391] Vautherin D. and Brink D. M., Hartree–Fock calculations with Skyrme’s interaction (I): Spherical nuclei, Phys. Rev. C 5(3), 626 (1972)

[392] Chabanat E., Bonche P., Haensel P., Meyer J., and Schaeffer R., A Skyrme parametrization from subnuclear to neutron star densities, Nucl. Phys. A 627(4), 710 (1997)

[393] Dechargé J. and Gogny D., Hartree–Fock–Bogolyubov calculations with the D1 effective interaction on spherical nuclei, Phys. Rev. C 21(4), 1568 (1980)

[394] Kleban M., Nerlo-Pomorska B., Berger J. F., Decharge J., Girod M., and Hilaire S., Global properties of spherical nuclei obtained from Hartree–Fock–Bogoliubov calculations with the Gogny force, Phys. Rev. C 65(2), 024309 (2002)

[395] Cochet B., Bennaceur K., Meyer J., Bonche P., and Duguet T., Skyrme forces with extended density dependence, Int. J. Mod. Phys. E 13(01), 187 (2004)

[396] Reinhard P. G. and Bender M., Mean field: Relativistic versus non-relativistic, in: G. A. Lalazissis, P. Ring, and D. Vretenar (Eds.), Extended density functionals in nuclear structure physics, Lect. Notes Phys. 641, 249 (2004)

[397] Serot B. D. and Walecka J. D., Recent progress in quantum hadrodynamics, Int. J. Mod. Phys. E 6(04), 515 (1997)

[398] Furnstahl R. J., Next generation relativistic models, in: A. Lalazissis, P. Ring, and D. Vretenar (Eds.), Extended density functionals in nuclear structure physics, Lect. Notes Phys. 641, 1 (2004)

[399] Lutz M., Friman B., and Appel Ch., Saturation from nuclear pion dynamics, Phys. Lett. B 474(1–2), 7 (2000)

[400] Finelli P., Kaiser N., Vretenar D., and Weise W., Nuclear many-body dynamics constrained by QCD and chiral symmetry, Eur. Phys. J. A 17(4), 573 (2003)

[401] Finelli P., Kaiser N., Vretenar D., and Weise W., Relativistic nuclear model with point-couplings constrained by QCD and chiral symmetry, Nucl. Phys. A 735(3–4), 449 (2004)

[402] Vretenar D. and Weise W., Exploring the nucleus in the context of low-energy QCD, in: G. A. Lalazissis, P. Ring, and D. Vretenar (Eds.), Extended density functionals in nuclear structure physics, Lect. Notes Phys. 641, 65 (2004)

[403] Dobaczewski J., Ab initio derivation of model energy density functionals, J. Phys. G 43(4), 04LT01 (2016)

[404] Bonnard J., Grasso M., and Lacroix D., Energy-density functionals inspired by effective-field theories: Applications to neutron drops, Phys. Rev. C 98(3), 034319 (2018)

[405] Pandharipande V. R. and Wiringa R. B., Variations on a theme of nuclear matter, Rev. Mod. Phys. 51(4), 821 (1979)

[406] Akmal A., Pandharipande V. R., and Ravenhall D. G., Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58(3), 1804 (1998)

[407] Brueckner K. A. and Gammel J. L., Properties of nuclear matter, Phys. Rev. 109(4), 1023 (1958)

[408] Jaminon M. and Mahaux C., Effective masses in relativistic approaches to the nucleon–nucleus mean field, Phys. Rev. C 40(1), 354 (1989)

[409] Zuo W., Lejeune A., Lombardo U., and Mathiot J. F., Interplay of three-body interactions in the EOS of nuclear matter, Nucl. Phys. A 706(3–4), 418 (2002)

[410] Zhou X. R., Burgio G. F., Lombardo U., Schulze H. J., and Zuo W., Three-body forces and neutron star structure, Phys. Rev. C 69(1), 018801 (2004)

[411] Haar B. and Malfliet R., Nucleons, mesons and deltas in nuclear matter a relativistic Dirac–Brueckner approach, Phys. Rep. 149(4), 207 (1987)

[412] Brockmann R. and Machleidt R., Relativistic nuclear structure (I): Nuclear matter, Phys. Rev. C 42(5), 1965 (1990)

[413] de Jong F. and Lenske H., Relativistic Brueckner–Hartree–Fock calculations with explicit intermediate negative energy states, Phys. Rev. C 58(2), 890 (1998)

[414] Gross-Boelting T., Fuchs C., and Faessler A., Covariant representations of the relativistic Brueckner T-matrix and the nuclear matter problem, Nucl. Phys. A 648(1–2), 105 (1999)

[415] Schiller E. and Müther H., Correlations and the Dirac structure of the nucleon self-energy, Eur. Phys. J. A 11(1), 15 (2001)

[416] Fuchs C., The Relativistic Dirac-Brueckner Approach to Nuclear Matter, in: G. A. Lalazissis, P. Ring, and D. Vretenar (Eds.), Extended density functionals in nuclear structure physics, Lect. Notes Phys. 641, 119 (2004)

[417] van Dalen E., Fuchs C., and Faessler A., The relativistic Dirac–Brueckner approach to asymmetric nuclear matter, Nucl. Phys. A 744, 227 (2004)

[418] Dalen E. N. E., Fuchs C., and Faessler A., Momentum, density, and isospin dependence of symmetric and asymmetric nuclear matter properties, Phys. Rev. C 72(6), 065803 (2005)

[419] Müther H. and Polls A., Two-body correlations in nuclear systems, Prog. Part. Nucl. Phys. 45(1), 243 (2000)

[420] Dickhoff W. H. and Barbieri C., Self-consistent Green’s function method for nuclei and nuclear matter, Prog. Part. Nucl. Phys. 52(2), 377 (2004)

[421] Carlson J., Morales J., Pandharipande V. R., and Ravenhall D. G., Quantum Monte Carlo calculations of neutron matter, Phys. Rev. C 68(2), 025802 (2003)

[422] Li B. A., Chen L. W., Yong G. C., and Zuo W., Double neutron/proton ratio of nucleon emissions in isotopic reaction systems as a robust probe of nuclear symmetry energy, Phys. Lett. B 634(4), 378 (2006)

[423] Famiano M. A., Liu T., Lynch W. G., Mocko M., Rogers A. M., et al., Neutron and proton transverse emission ratio measurements and the density dependence of the asymmetry term of the nuclear equation of state, Phys. Rev. Lett. 97(5), 052701 (2006)

[424] Tsang M. B., Liu T. X., Shi L., Danielewicz P., Gelbke C. K., et al., Isospin diffusion and the nuclear symmetry energy in heavy ion reactions, Phys. Rev. Lett. 92(6), 062701 (2004)

[425] Liu T. X., Lynch W. G., Tsang M. B., Liu X. D., Shomin R., et al., Isospin diffusion observables in heavy-ion reactions, Phys. Rev. C 76(3), 034603 (2007)

[426] Kohley Z., May L. W., Wuenschel S., Bonasera A., Hagel K., et al., Investigation of transverse collective flow of intermediate mass fragments, Phys. Rev. C 82(6), 064601 (2010)

[427] Zhang Y., Tian J., Cheng W., Guan F., Huang Y., et al., Long-time drift of the isospin degree of freedom in heavy ion collisions, Phys. Rev. C 95, 041602(R) (2017)

[428] Russotto P., Cozma M. D., Le Fèvre A., Leifels Y., Lemmon R., Li Q., Łukasik J., and Trautmann W., Flow probe of symmetry energy in relativistic heavy-ion reactions, Eur. Phys. J. A 50(2), 38 (2014)

[429] Russotto P., Wu P. Z., Zoric M., Chartier M., Leifels Y., Lemmon R. C., Li Q., Łukasik J., Pagano A., Pawłowski P., and Trautmann W., Symmetry energy from elliptic flow in ¹⁹⁷Au+¹⁹⁷Au, Phys. Lett. B 697(5), 471 (2011)

[430] Tsang M. B., Zhang Y., Danielewicz P., Famiano M., Li Z., Lynch W. G., and Steiner A. W., Constraints on the density dependence of the symmetry energy, Phys. Rev. Lett. 102(12), 122701 (2009)

[431] Danielewicz P. and Lee J., Symmetry energy (I): Semi-infinite matter, Nucl. Phys. A 818(1–2), 36 (2009)

[432] Horowitz C. J. and Piekarewicz J., Neutron star structure and the neutron radius of ²⁰⁸Pb, Phys. Rev. Lett. 86(25), 5647 (2001)

[433] Piekarewicz J., Unmasking the nuclear matter equation of state, Phys. Rev. C 69(4), 041301 (2004)

[434] Yoshida S. and Sagawa H., Isovector nuclear matter properties and neutron skin thickness, Phys. Rev. C 73(4), 044320 (2006)

[435] Galichet E., et al. (INDRA Collaboration), Isospin diffusion in ⁵⁸Ni-induced reactions at intermediate energies (I): Experimental results, Phys. Rev. C 79(6), 064614 (2009)

[436] Wang R. S., Zhang Y., Xiao Z. G., Tian J. L., Zhang Y. X., et al., Time-dependent isospin composition of particles emitted in fission events following ⁴⁰Ar+¹⁹⁷Au at 35 MeV/u, Phys. Rev. C 89(6), 064613 (2014)

[437] Tsang M. B., Stone J. R., Camera F., Danielewicz P., Gandolfi S., et al., Constraints on the symmetry energy and neutron skins from experiments and theory, Phys. Rev. C 86(1), 015803 (2012)

[438] Lattimer J. M. and Steiner A. W., Constraints on the symmetry energy using the mass-radius relation of neutron stars, Eur. Phys. J. A 50(2), 40 (2014)

[439] Li B. A. and Chen L. W., Nucleon–nucleon cross sections in neutron-rich matter and isospin transport in heavy-ion reactions at intermediate energies, Phys. Rev. C 72(6), 064611 (2005)

[440] Klimkiewicz A., Paar N., Adrich P., Fallot M., Boretzky K., et al., Nuclear symmetry energy and neutron skins derived from pygmy dipole resonances, Phys. Rev. C 76(5), 051603 (2007)

[441] Colò G. and Danielewicz P., Constraints, limits and extensions for nuclear energy functionals, AIP Conf. Proc. 1128, 59 (2009)

[442] Trippa L., Colò G., and Vigezzi E., Giant dipole resonance as a quantitative constraint on the symmetry energy, Phys. Rev. C 77, 061304(R)

[443] Ou L., Xiao Z. G., Yi H., Wang N., Liu M., and Tian J., Dynamic isovector reorientation of deuteron as a probe to nuclear symmetry energy, Phys. Rev. Lett. 115(21), 212501 (2015)

[444] Köhler S., Skyrme force and the mass formula, Nucl. Phys. A 258(2), 301 (1976)

[445] Tondeur F., Brack M., Farine M., and Pearson J. M., Static nuclear properties and the parametrisation of Skyrme forces, Nucl. Phys. A 420(2), 297 (1984)

[446] Margueron J., Navarro J., and Van Giai N., Instabilities of infinite matter with effective Skyrme-type interactions, Phys. Rev. C 66(1), 014303 (2002)

[447] Wu Q., Zhang Y., Xiao Z., Wang R., Zhang Y., Li Z., Wang N., and Showalter R. H., Competition between Coulomb and symmetry potential in semi-peripheral heavy ion collisions, Phys. Rev. C 91(1), 014617 (2015)

[448] Gaitanos T., Di Toro M., Ferini G., Colonna M., and Wolter H. H., Isospin effects in intermediate energy heavy ion collisions, arXiv: nucl-th/0402041 (2004)

[449] Bass S. A., Hartnack C., Stöcker H., and Greiner W., High p t pions as probes of the dense phase of relativistic heavy ion collisions, Phys. Rev. C 50(4), 2167 (1994)

[450] Xiao Z., Li B. A., Chen L. W., Yong G. C., and Zhang M., Circumstantial evidence for a soft nuclear symmetry energy at suprasaturation densities, Phys. Rev. Lett. 102(6), 062502 (2009)

[451] Xie W. J., Su J., Zhu L., and Zhang F. S., Symmetry energy and pion production in the Boltzmann–Langevin approach, Phys. Lett. B 718(4–5), 1510 (2013)

[452] Feng Z. Q. and Jin G. M., Probing high-density behavior of symmetry energy from pion emission in heavy-ion collisions, Phys. Lett. B 683(2–3), 140 (2010)

[453] Song T. and Ko C. M., Modifications of the pion-production threshold in the nuclear medium in heavy ion collisions and the nuclear symmetry energy, Phys. Rev. C 91(1), 014901 (2015)

[454] Cozma M. D., Constraining the density dependence of the symmetry energy using the multiplicity and average p_T ratios of charged pions, Phys. Rev. C 95(1), 014601 (2017)

[455] Zhang Z. and Ko C. M., Medium effects on pion production in heavy ion collisions, Phys. Rev. C 95(6), 064604 (2017)

[456] Hong J. and Danielewicz P., Subthreshold pion production within a transport description of central Au+Au collisions, Phys. Rev. C 90(2), 024605 (2014)

[457] Ou L., Li Z., Zhang Y., and Liu M., Effect of the splitting of the neutron and proton effective masses on the nuclear symmetry energy at finite temperatures, Phys. Lett. B 697(3), 246 (2011)

[458] Xu J., Chen L. W., Li B. A., and Ma H. R., Temperature effects on the nuclear symmetry energy and symmetry free energy with an isospin and momentum dependent interaction, Phys. Rev. C 75(1), 014607 (2007)

[459] Xu J., Chen L. W., Li B. A., and Ma H. R., Effects of isospin and momentum dependent interactions on thermal properties of asymmetric nuclear matter, Phys. Rev. C 77(1), 014302 (2008)

[460] Xie W. J. and Li B. A., Bayesian inference of high-density nuclear symmetry energy from radii of canonical neutron stars, Astrophys. J. 883(2), 174 (2019)

[461] Margueron J., Hoffmann Casali R., and Gulminelli F., Equation of state for dense nucleonic matter from metamodeling (I): Foundational aspects, Phys. Rev. C 97(2), 025805 (2018)

[462] Margueron J. and Gulminelli F., Effect of high-order empirical parameters on the nuclear equation of state, Phys. Rev. C 99(2), 025806 (2019)

[463] Zhang N. B. and Li B. A., Delineating effects of nuclear symmetry energy on the radii and tidal polarizabilities of neutron stars, J. Phys. G 46(1), 014002 (2019)

[464] Agrawal B. K., Shlomo S., and Au V. K., Determination of the parameters of a Skyrme type effective interaction using the simulated annealing approach, Phys. Rev. C 72(1), 014310 (2005)

[465] Chen L. W., Cai B. J., Ko C. M., Li B. A., Shen C., and Xu J., Higher-order effects on the incompressibility of isospin asymmetric nuclear matter, Phys. Rev. C 80(1), 014322 (2009)

[466] Mondal C., Agrawal B. K., De J. N., Samaddar S. K., Centelles M., and Viñas X., Interdependence of different symmetry energy elements, Phys. Rev. C 96, 021302(R) (2017)

[467] Dutra M., Lourenco O., Sa Martins J. S., Delfino A., Stone J. R., and Stevenson P. D., Skyrme interaction and nuclear matter constraints, Phys. Rev. C 85(3), 035201 (2012)

[468] Guichon P. A. M. and Thomas A. W., Quark structure and nuclear effective forces, Phys. Rev. Lett. 93(13), 132502 (2004)

[469] Stone J. R., Stone N. J., and Moszkowski S. A., Incompressibility in finite nuclei and nuclear matter, Phys. Rev. C 89(4), 044316 (2014)

[470] Zhang Y., Liu M., Xia C. J., Li Z., and Biswal S. K., Constraints on the symmetry energy and its associated parameters from nuclei to neutron stars, Phys. Rev. C 101(3), 034303 (2020)

[471] Rizzo J., Colonna M., Baran V., Di Toro M., Wolter H. H., and Zielinska-Pfabe M., Isospin dynamics in peripheral heavy ion collisions at Fermi energies, Nucl. Phys. A 806(1–4), 79 (2008)

[472] Li B. A., Das C. B., Das Gupta S., and Gale C., Momentum dependence of the symmetry potential and nuclear reactions induced by neutron-rich nuclei at RIA, Phys. Rev. C 69(1), 011603 (2004)

[473] Giordano V., Colonna M., DiToro M., Greco V., and Rizzo J., Isospin emission and flow at high baryon density: A test of the symmetry potential, Phys. Rev. C 81(4), 044611 (2010)

[474] Reinhard P. G. and Flocard H., Nuclear effective forces and isotope shifts, Nucl. Phys. A 584(3), 467 (1995)

[475] Friedrich J. and Reinhard P. G., Skyrme-force parametrization: Least-squares fit to nuclear ground-state properties, Phys. Rev. C 33(1), 335 (1986)

[476] Coupland D. D. S., Lynch W. G., Tsang M. B., Danielewicz P., and Zhang Y., Influence of transport variables on isospin transport ratios, Phys. Rev. C 84(5), 054603 (2011)

[477] Coupland D. D. S., PhD thesis, Michigan State University, 2013

[478] Coupland D. D. S., Youngs M., Chajecki Z., Lynch W. G., Tsang M. B., et al., Probing effective nucleon masses with heavy-ion collisions, Phys. Rev. C 94(1), 011601 (2016)

[479] Brandt S., Data Analysis: Statistical and Computational Methods for Scientists and Engineers, 4th Ed., Springer, 2014

[480] Morfouace P., Tsang C. Y., Zhang Y., Lynch W. G., Tsang M. B., et al., Constraining the symmetry energy with heavy-ion collisions and Bayesian analyses, Phys. Lett. B 799, 135045 (2019)

[481] Li L., Zhang Y., Li Z., Wang N., Cui Y., and Winkelbauer J., Impact parameter smearing effects on isospin sensitive observables in heavy ion collisions, Phys. Rev. C 97(4), 044606 (2018)

[482] Sun Z. Y., Tsang M. B., Lynch W. G., Verde G., Amorini F., et al., Isospin diffusion and equilibration for Sn+Sn collisions at E/A=35 MeV, Phys. Rev. C 82, 051603(R) (2010)

[483] Kolomeitsev E. E., Hartnack C., Barz H. W., Bleicher M., Bratkovskaya E., et al., Transport theories for heavy-ion collisions in the 1A GeV regime, J. Phys. G 31(6), S741 (2005)

[484] Xu J., Chen L. W., Tsang M. B., Wolter H., Zhang Y. X., et al., Understanding transport simulations of heavy-ion collisions at 100A and 400A MeV: Comparison of heavy-ion transport codes under controlled conditions, Phys. Rev. C 93(4), 044609 (2016)

[485] Desai V. V., Loveland W., McCaleb K., Yanez R., Lane G., et al., The ¹³⁶Xe+¹⁹⁸Pt reaction: A test of models of multi-nucleon transfer reactions, Phys. Rev. C 99(4), 044604 (2019)

[486] Ono A., Dynamics of clusters and fragments in heavy-ion collisions, Prog. Part. Nucl. Phys. 105, 139 (2019)

Rights & permissions

Higher Education Press 2020