# Influence of the treatment of initialization and mean-field potential on the
neutron to proton yield ratios

Junping Yang  China Institute of Atomic Energy, Beijing 102413, China
Department of Physics, Guangxi Normal University, Guilin, 541004, China
Yingxun Zhang<EMAIL_ADDRESS>China Institute of Atomic Energy, Beijing
102413, China  Guangxi Key Laboratory of Nuclear Physics and Technology,
Guangxi Normal University, Guilin, 541004, China  Ning Wang  Guangxi Key
Laboratory of Nuclear Physics and Technology, Guangxi Normal University,
Guilin, 541004, China  Department of Physics, Guangxi Normal University,
Guilin, 541004, China  Zhuxia Li  China Institute of Atomic Energy, Beijing
102413, China

###### Abstract

In this work, we firstly investigate how to reproduce and how well one can
reproduce the Woods-Saxon density distribution of initial nuclei in the
framework of the improved quantum molecular dynamics model. Then, we propose a
new treatment for the initialization of nuclei which is correlated with the
nucleonic mean-field potential by using the same potential energy density
functional. In the mean field potential, the three-body force term is
accurately calculated. Based on the new version of the model, the influences
of precise calculations of the three-body force term, the slope of symmetry
energy, the neutron-proton effective mass splitting, and the width of the wave
packet on heavy ion collision observables, such as the neutron to proton yield
ratios for emitted free nucleons [$R(n/p)$] and for coalescence invariant
nucleons [$R_{ci}(n/p)$] for 124Sn+112Sn at the beam energy of 200 MeV per
nucleon, are discussed. Our calculations show that the spectra of neutron to
proton yield ratios [$R(n/p)$] can be used to probe the slope of symmetry
energy ($L$) and the neutron-proton effective mass splitting. In detail, the
$R(n/p)$ in the low kinetic energy region can be used to probe the slope of
symmetry energy ($L$). With a given $L$, the inclination of $R(n/p)$ to
kinetic energy ($E_{k}$) can be used to probe the effective mass splitting. In
the case where the neutron-proton effective mass splitting is fixed, $R(n/p)$
at high kinetic energy can also be used to learn the symmetry energy at
suprasaturation density.

###### pacs:

21.60.Jz, 21.65.Ef, 24.10.Lx, 25.70.-z

##  I Introduction

The isospin asymmetric nuclear equation of state is fundamental for
understanding the objectives of both nuclear physics, such as the properties
of neutron-rich nuclei[1, 2, 3, 4, 5, 6, 7, 8] and the mechanism of heavy ion
collisions[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 9], and astrophysics,
such as the properties of neutron star masses, radii, and tidal
deformability[21, 22, 23, 24, 25, 26, 27, 28, 29]. However, the theoretical
predictions present large uncertainties on the isospin asymmetric nuclear
equation of state away from normal density, especially the density dependence
of the symmetry energy. In the laboratory, low-intermediate energy heavy ion
collisions (HICs) can provide the constraints of symmetry energy from
subnormal density to twice the saturation density by comparing the data of HIC
observables, such as neutron to proton yield ratios[10, 14, 15], triton to He3
yield ratios[30, 31], isospin diffusion[11, 12, 9], neutron excess[32],
$\pi^{-}/\pi^{+}$ ratios[33, 34], collective flow[35, 36], and so on, with the
transport model calculations.

Among those transport model calculations of isospin sensitive HIC observables,
the initialization of nuclei plays an important role, as does the isospin
dependent mean-field potential. For example, consideration of the neutron skin
in the initialization can influence the prediction of $\pi^{-}/\pi^{+}$
ratios[37, 38] in peripheral HICs. In these pioneer calculations, the slope of
symmetry energy and the thickness of neutron skin were treated separately.
However, the theoretical calculations show that the neutron skin thickness of
a heavy nucleus, which is given by the difference between the root-mean-square
radii of neutrons and protons, i.e., $\Delta r_{np}=\langle
r^{2}\rangle^{1/2}_{n}-\langle r^{2}\rangle^{1/2}_{p}$, is strongly correlated
to the slope of symmetry energy $L$[3, 4, 5, 9]. Thus, a consistent treatment
of neutron skin in the initialization and isospin dependent mean-field
potential in nucleon propagation is highly desired in the development of
transport models to reduce the uncertainties of symmetry energy constraint
caused by separately treating $L$ and $\Delta r_{np}$.

In the QMD type models, each nucleon is represented by a Gaussian wave packet,

$\phi_{i}(\mathbf{r})=\frac{1}{(2\pi\sigma_{r}^{2})^{3/4}}\exp{\\{-\frac{[\mathbf{r}-\mathbf{r}_{i}(t)]^{2}}{4\sigma_{r}^{2}}+i\mathbf{r}_{i}\cdot\mathbf{p}_{i}(t)\\}}.$
(1)

The initial nuclei are prepared by sampling he centroids of wave packets
$\mathbf{r}_{i}(t=0)$ in a hard sphere[42, 43, 44, 39, 40, 41]. As a result,
the sampled density profile of the nucleus has a larger tail than the required
density due to the finite width of the Gaussian wave packet. This is the main
difficulty in reproducing the density distribution calculated by microscopic
nuclear models. Even for the simple Woods-Saxon (WS) density distribution,
which is wildly used in transport model simulations, it is still hard to
reproduce and has been discussed in the transport model comparison
project[42]. There were some efforts to improve the description of the density
profiles in the initialization. For example, in the ultrarelativistic quantum
molecular dynamics (UrQMD) model, the centroids of wave packets are sampled
within the radius which is reduced by half a layer of nucleons from the
original nuclear radius[43]. The sampled density profile is similar to a
Woods-Saxon density distribution, but they have different diffuseness values.
In the parton-hadron quantum molecular dynamics (PHQMD)[44] and Tübingen
quantum molecular dynamics (TuQMD) models[45], a small width of Gaussian wave
packet is used to well reproduce the Woods-Saxon density distribution. Thus,
an important theoretical question in the QMD type models is how to reproduce
and how well one can reproduce the required density profile of the initial
nucleus.

In this work, we investigate how to reproduce the initial density distribution
of the nucleus by doing the inverse Weierstrass transformation, and then redo
the treatment of the initialization and mean-field potential in the improved
quantum molecular dynamics (ImQMD) model with the same Skyrme energy density
functional. The symmetry energy and neutron-proton effective mass splitting
effects on the isospin sensitive observables, such as neutron to proton yield
ratios, for the reaction of 124Sn+112Sn at the beam energy of 200 MeV/u are
examined and discussed. The new calculations show that the spectra of neutron
to proton yield ratios can be used to probe $L$ and the neutron-proton
effective mass splitting. In the case where the neutron-proton effective mass
splitting is fixed, neutron to proton yield ratios at kinetic energy greater
than 100 MeV can be used to probe the symmetry energy at suprasaturation
density.

##  II Theoretical approaches

For describing HICs, the ImQMD model is adopted. The treatments of nucleon-
nucleon collision, Pauli blocking, etc. are as the same as in the previous
version of ImQMD. More details about them can be found in the review paper
[46]. The improvements we made in this work are in the initialization and
nucleonic mean-field potential.

###  II.1 Nucleonic Mean-field potential

In the ImQMD model [9], the Skyrme type nucleonic potential energy density
without the spin-orbit term is used:

$\displaystyle u(\mathbf{r})_{sky}=$
$\displaystyle\frac{\alpha}{2}\frac{\rho^{2}}{\rho_{0}}+\frac{\beta}{\eta+1}\frac{\rho^{\eta+1}}{\rho_{0}^{\eta}}+\frac{g_{sur}}{2\rho_{0}}(\nabla\rho)^{2}$
$\displaystyle+\frac{g_{sur,iso}}{\rho_{0}}[\nabla(\rho_{n}-\rho_{p})]^{2}$
$\displaystyle+A_{sym}\frac{\rho^{2}}{\rho_{0}}\delta^{2}+B_{sym}\frac{\rho^{\eta+1}}{\rho_{0}^{\eta}}\delta^{2}+u_{md}.$

The $\alpha$ is the parameter related to the two-body term, $\beta$ and $\eta$
are related to the three-body term, and $A_{sym}$ and $B_{sym}$ are the
coefficients in the symmetry potential and come from the two- and the three-
body interaction terms. The density $\rho$ is obtained by integrating the
Wigner phase space density
$f_{i}(\mathbf{r},\mathbf{p})=\frac{1}{(\pi\hbar)^{3}}e^{-\frac{(\mathbf{r}-\mathbf{r}_{i})}{2\sigma_{r}^{2}}-\frac{(\mathbf{p}-\mathbf{p}_{i})^{2}}{2\sigma_{p}^{2}}}$
in momentum space for all nucleons, i.e.,

$\displaystyle\rho(\mathbf{r})=\sum_{i=1}^{A}\rho_{i}(\mathbf{r})$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{A}\int
f_{i}(\mathbf{r},\mathbf{p})d^{3}\mathbf{p}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{A}\frac{1}{(2\pi\sigma_{r}^{2})^{3/2}}e^{-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{2\sigma_{r}^{2}}},$

and $\delta=(\rho_{n}-\rho_{p})/(\rho_{n}+\rho_{p})$. Here, $\mathbf{r}_{i}$
and $\mathbf{p}_{i}$ are centroids of distribution, and also the variational
parameters of Gaussian single-particle wave functions [46].

The Skyrme-type momentum dependent energy density functional $u_{md}$ is
written based on its interaction form
$\delta(\mathbf{r}_{1}-\mathbf{r}_{2})(\mathbf{p}_{1}-\mathbf{p}_{2})^{2}$
[47]. In the ImQMD model,

$\displaystyle u_{md}(\mathbf{r},\\{\mathbf{p}_{i}-\mathbf{p}_{j}\\})$
$\displaystyle=$ $\displaystyle C_{0}\sum_{ij}\int
d^{3}pd^{3}p^{\prime}f_{i}(\mathbf{r},\mathbf{p})f_{j}(\mathbf{r},\mathbf{p}^{\prime})(\mathbf{p}-\mathbf{p}^{\prime})^{2}$
$\displaystyle+D_{0}\sum_{ij\in n}\int
d^{3}pd^{3}p^{\prime}f_{i}(\mathbf{r},\mathbf{p})f_{j}(\mathbf{r},\mathbf{p}^{\prime})(\mathbf{p}-\mathbf{p}^{\prime})^{2}$
$\displaystyle+D_{0}\sum_{ij\in p}\int
d^{3}pd^{3}p^{\prime}f_{i}(\mathbf{r},\mathbf{p})f_{j}(\mathbf{r},\mathbf{p}^{\prime})(\mathbf{p}-\mathbf{p}^{\prime})^{2}.$

$C_{0}$, $D_{0}$ are parameters related to the momentum dependent interaction.

The parameters in Eqs.(II.1) and (II.1) can be obtained from the standard
Skyrme interaction parameters as in Refs.[48, 9]. The connection between the
seven parameters $\alpha$, $\beta$, $\eta$, $A_{sym}$, $B_{sym}$, $C_{0}$,
$D_{0}$ used in the ImQMD model and the seven nuclear matter parameters-the
saturation density $\rho_{0}$, binding energy at saturation density $E_{0}$,
incompressibility $K_{0}$, symmetry energy coefficient $S_{0}$, the slope of
symmetry energy $L$, isocalar effective mass $m_{s}^{*}$, and isovector
effective mass $m_{v}^{*}$-are given in Ref. [9]. Thus, one can alternatively
use $\rho_{0}$, $E_{0}$, $K_{0}$, $S_{0}$, $L$, $m_{s}^{*}$, $m_{v}^{*}$, as
input to study the influence of different nuclear matter parameters.

Based on recent constraints of nuclear matter parameters related to symmetry
energy from nuclei to neutron stars[9], we choose $K_{0}=240$ MeV, $S_{0}=30$
MeV, $m_{v}^{\ast}/m=0.7$, and $m_{s}^{\ast}/m=0.8$, which are also in the
reasonable region for the Skyrme type force[49]. The $g_{sur}$ and
$g_{sur,iso}$ values are taken as 24.5 and -4.99 MeVfm2, respectively. $L$ is
varied from 30 to 110 MeV to analyze the influence of $L$ on the spectra of
n/p ratios. These values are the default parameters we used, and are listed in
Table 1. In addition, the influence of $S_{0}$, $m_{s}^{\ast}$, and the
neutron-proton effective mass splitting $\Delta
m_{np}^{\ast}=m_{n}^{\ast}-m_{p}^{\ast}$ on the n/p ratios are also discussed.
One point we want to mention is that $\Delta
m_{np}^{\ast}=m_{n}^{\ast}-m_{p}^{\ast}$ is alternatively described by
$f_{I}=(m/m_{s}^{\ast}-m/m_{v}^{\ast})=\frac{1}{2\delta}(m/m_{n}^{\ast}-m/m_{p}^{\ast})$
in the following discussions, because it can be analytically incorporated into
the transport code and is independent of isospin asymmetry.

Table 1: Values of the nuclear matter parameter used in the ImQMD-L model. $\rho_{0}$ is in fm-3, $E_{0}$, $K_{0}$, $S_{0}$, $L$ are in MeV, $g_{sur}$ and $g_{sur,iso}$ are in MeVfm2. $f_{I}=\frac{m}{m_{s}^{\ast}}-\frac{m}{m_{v}^{\ast}}$=-0.178. $\ K_{0}$ | $S_{0}$ | $E_{0}$ | $\rho_{0}$ | $m_{v}^{\ast}/m$ | $m_{s}^{\ast}/m$ | $g_{sur}$ | $g_{sur,iso}$ | $L$
---|---|---|---|---|---|---|---|---
240 | 30 | -16 | 0.16 | 0.7 | 0.8 | 24.5 | -4.99 | 30,50,70,90,110

The nucleonic force acting on the $i$th nucleon is

$F_{i}=\dot{\mathbf{p}}_{i}=-\frac{\partial
U(\\{\mathbf{r}_{i},\mathbf{p}_{i}\\})}{\partial\mathbf{r}_{i}}.$ (5)

In the above formula, the potential energy
$U(\\{\mathbf{r}_{i},\mathbf{p}_{i}\\})$ is obtained by integrating the
potential energy density in coordinate space, i.e.,
$U(\\{\mathbf{r}_{i},\mathbf{p}_{i}\\})=\int u(\mathbf{r})d^{3}\mathbf{r}$,
and it is a function of $\mathbf{r}_{i}$ and $\mathbf{p}_{i}$. Because the
density has a Gaussian form, all the integrals in potential energy
calculations can be done analytically except for the three-body related terms
with $\eta\neq 2$. When $\eta\neq 2$, one has to solve it with an
approximation or numerical method. Usually, the three-body related potential
energy is approximately calculated as

$\displaystyle\frac{\beta}{\eta+1}\int\frac{\rho^{\eta}}{\rho_{0}^{\eta}}\rho
d^{3}\mathbf{r}$ $\displaystyle=$
$\displaystyle\frac{\beta}{\eta+1}\sum_{i=1}^{N}<\frac{\rho^{\eta}}{\rho_{0}^{\eta}}>_{i}$
$\displaystyle\approx$
$\displaystyle\frac{\beta}{\eta+1}\sum_{i=1}^{N}<\frac{\rho}{\rho_{0}}>_{i}^{\eta}+O(\rho)$

by omitting the higher order term[50, 51, 52, 46] or adjusting the value of
$\sigma_{r}$ in the three-body term[41, 44]. The three-body force acting on
particle $i$ is calculated according to the followings expression:

$\dot{\mathbf{p}}_{i}=-\frac{\partial
U_{3}}{\partial\mathbf{r}_{i}}=-\frac{\beta}{\eta+1}\frac{\partial}{\partial\mathbf{r_{i}}}\sum_{j=1}^{N}\frac{\langle\rho\rangle^{\eta}_{j}}{\rho_{0}^{\eta}}.$
(7)

In the above formula,
$\langle\rho\rangle_{j}=\sum_{k}\frac{1}{(4\pi\sigma_{r}^{2})^{3/2}}e^{-\frac{(\mathbf{r}_{j}-\mathbf{r}_{k})^{2}}{4\sigma_{r}^{2}}}$.
For uniform matter, this approximation is good enough since the density
fluctuation is zero.

However, the density variance in intermediate energy HICs is strong and time
dependent. Thus, the influence of the higher order term in Eq. (II.1) should
not be neglected. In this work we exactly calculate the three-body term by
using the numerical quadrature method. Further, the force acting on particle
$i$ related to the three-body term is calculated as

$\dot{\mathbf{p}}_{i}=-\frac{\partial
U_{3}}{\partial\mathbf{r}_{i}}=-\beta\rho_{0}\int\frac{\rho^{\eta}}{\rho_{0}^{\eta}}\frac{\rho_{i}}{\rho_{0}}\frac{\mathbf{r}-\mathbf{r}_{i}}{\sigma_{r}^{2}}d^{3}\mathbf{r}.$
(8)

The integral in Eq. (8) is solved by using the 11-point Gauss-Legendre
quadrature method, and it gives a stronger three-body force in ImQMD-L than in
ImQMD. To distinguish with the previous version of the ImQMD model, we named
it ImQMD-L (L, means the lattice method) in the following discussions.

###  II.2 Initialization of nucleus in ImQMD-L

In this section, we will introduce how to obtain the neutron and proton Woods-
Saxon density distribution of the initial nucleus with the same Skyrme energy
density functional as in the mean-field propagation in the ImQMD-L model.
Then, we investigate the criteria for reproducing the Woods-Saxon density
distribution with a Gaussian wave packet. Then, the influences of different
widths of wave packet, which lead to the different shapes of Woods-Saxon
density distributions, on the stability and binding energy of the initial
nucleus are discussed. Finally, we describe the method of initialization used
in the ImQMD-L model.

####  II.2.1 Density distribution of the initial nucleus with restricted
density variational method

In the calculations, we take the density distribution as a Woods-Saxon density
function with

$\rho_{i}=\rho_{0i}\frac{1}{1+\exp(\frac{r-R_{i}}{a_{i}})},i={n,p}.$ (9)

Here, $\rho_{0n}$, $\rho_{0p}$, $R_{p}$, $R_{n}$, $a_{p}$, and $a_{n}$ are the
saturation density, radius, and diffuseness values of proton and neutron
density distributions, and they are obtained by minimizing the total energy of
the system given by,

$E=\int\mathcal{H}dr=\int\\{\frac{\hbar^{2}}{2m}[\tau_{n}(\mathbf{r})+\tau_{p}(\mathbf{r})]+u_{sky}+u_{coul}\\}dr,$
(10)

under the condition of the conservation of particle number in the system. This
method is named as the restricted density variational method (RDV) method
[53].

The same semiclassical expression of the Skyrme energy density functional as
in ImQMD, i.e., $u_{sky}$, is applied. One should note, that $u_{md}$ is
reduced to

$u_{md}=\frac{C_{0}}{2\hbar^{2}}\rho\tau+\frac{D_{0}}{2\hbar^{2}}(\rho_{n}\tau_{n}+\rho_{p}\tau_{p})$
(11)

in the calculations of the nucleus in its ground state. The kinetic energy
density $\tau_{i}$ in the RDV method is given by

$\displaystyle\tau_{i}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{3}{5}(3\pi^{2})^{2/3}\rho_{i}^{5/3}+\frac{1}{36}\frac{(\nabla\rho_{i})^{2}}{\rho_{i}}+\frac{1}{3}\triangle\rho_{i}$
$\displaystyle+\frac{1}{6}\frac{\nabla\rho_{i}\nabla f_{i}+\rho_{i}\triangle
f_{i}}{f_{i}}-\frac{1}{12}\rho_{i}(\frac{\nabla f_{i}}{f_{i}})^{2}$
$\displaystyle+\frac{1}{2}\rho_{i}(\frac{2m}{\hbar^{2}}\frac{W_{0}}{2}\frac{\nabla(\rho+\rho_{i})}{f_{i}})^{2},$

where we use the extended Thomas-Fermi (ETF) approach including all terms up
to second order (ETF2) and fourth-order (ETF4) as in Ref. [54]. $\rho_{i}$
denotes the proton and neutron densities of the nucleus, and
$\rho=\rho_{n}+\rho_{p}$. $W_{0}$ is the strength of the spin-orbit
interaction, and we set it to zero in order to use the same form of Skyrme
energy density functional as in ImQMD-L; the parameter $f_{i}(\mathbf{r})$ is
the same as in Ref. [53]. The calculated results of $a_{p}$, $R_{p}$, $a_{n}$,
$R_{n}$, binding energy $B$, and rms radii for neutron and proton obtained by
RDV with the Skyrme energy density functional are listed in Table 2.

Table 2: $a_{p}$, $R_{p}$, $a_{n}$, $R_{n}$, binding energy $B$, and rms radius for neutron and proton for 124Sn obtained with the RDV method. $L$ and $B$ are in MeV, $a_{p}$, $R_{p}$, $a_{n}$, $R_{n}$, ${<r_{p}^{2}>}^{1/2}$, and ${<r_{n}^{2}>}^{1/2}$ are in fm. The values in the bracket are the results obtained with given $a_{p}=a_{n}=f^{-1}(\sigma_{r}=1.287)$, for details, see Secs. II.2.2 and II.2.3. $L$ | $a_{p}(a^{\prime})$ | $R_{p}$ | $a_{n}(a^{\prime})$ | $R_{n}$ | $B(B^{*})$ | ${<r_{p}^{2}>}^{1/2}$ | ${<r_{n}^{2}>}^{1/2}$
---|---|---|---|---|---|---|---
30 | 0.414 | 5.733 | 0.514 | 5.777 | -7.971 | 4.700 | 4.865
| (0.743) |  | (0.743) |  | (-7.306) | (5.230) | (5.259)
50 | 0.415 | 5.729 | 0.507 | 5.811 | -8.021 | 4.698 | 4.881
| (0.743) |  | (0.743) |  | (-7.357) | (5.227) | (5.281)
70 | 0.419 | 5.707 | 0.503 | 5.838 | -8.073 | 4.687 | 4.893
| (0.743) |  | (0.743) |  | (-7.420) | (5.213) | (5.299)
90 | 0.422 | 5.686 | 0.496 | 5.872 | -8.129 | 4.676 | 4.907
| (0.743) |  | (0.743) |  | (-7.479) | (5.199) | (5.321)
110 | 0.426 | 5.656 | 0.487 | 5.909 | -8.191 | 4.659 | 4.922
| (0.743) |  | (0.743) |  | (-7.539) | (5.179) | (5.346)

####  II.2.2 Criteria for reproducing the Woods-Saxon density distribution

Following the wave function of a nucleon used in the quantum molecular
dynamics model, the nuclear density can be written as,

$\rho(\mathbf{r})=\sum_{i}^{A}\rho_{i}(\mathbf{r}),$ (13)

where

$\rho_{i}(\mathbf{r})=\frac{1}{(2\pi\sigma_{r}^{2})^{3/2}}\exp\\{-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{2\sigma_{r}^{2}}\\}.$
(14)

To reproduce the Woods-Saxon density profile of a nucleus, such as
$\rho^{ws}(\mathbf{r})$,

$\rho^{ws}(\mathbf{r})=\rho_{0}\frac{1}{1+\exp(\frac{|r|-R}{a})},$ (15)

where $\rho_{0}$ is the saturation density, $R$ is the half density radius,
and $a$ is the diffuseness parameter, by the method of Monte-Carlo sampling,
one has to know the distribution of the centroid of the wave packet,
$\tilde{\rho}(\mathbf{r}_{i})$, which satisfies

$\rho^{ws}(\mathbf{r})=\sum_{i}^{A}\int
d^{3}\mathbf{r}_{i}\tilde{\rho}(\mathbf{r}_{i})\frac{1}{(2\pi\sigma_{r}^{2})^{3/2}}e^{-\frac{(\mathbf{r}-\mathbf{r}_{i})^{2}}{2\sigma_{r}^{2}}}.$
(16)

For a spherical nucleus, the form of $\rho^{ws}(\mathbf{r})$ has spherical
symmetry, and thus the density distribution along the $x$ direction can be
written as

$\displaystyle\rho^{ws}(x)$ $\displaystyle=$
$\displaystyle\sum_{i}^{A}\int_{-\infty}^{\infty}dx_{i}\tilde{\rho}(x_{i})\frac{1}{(2\pi\sigma_{r}^{2})^{1/2}}e^{-\frac{(x-x_{i})^{2}}{2\sigma_{r}^{2}}}$
$\displaystyle=$
$\displaystyle\sum_{i}^{A}\int_{-\infty}^{\infty}dx_{i}\tilde{\rho}(x_{i})G(x-x_{i}).$

Equation (II.2.2) is a generalized Weierstrass transform[55] with a kernel
$G(x-x_{i})=\frac{1}{(2\pi\sigma_{r}^{2})^{1/2}}\exp\\{-\frac{(x-x_{i})^{2}}{2\sigma_{r}^{2}}\\}$,
where $x_{i}$ is the centroid of a Gaussian kernel. Thus, the key point for
the initialization of a nucleus with a Gaussian wave packet is to do an
inverse Weierstrass transform to get the solution of $\tilde{\rho}(x_{i})$ at
a given Gaussian kernel $G(x-x_{i})$.

Since the centroids of nucleons are sampled randomly within the density
distribution of the nucleus, $\tilde{\rho}(x_{i})$ can be thought as being the
same for all nucleons. Thus, Eq. (II.2.2) can be simplified as

$\rho^{ws}(x)=A\int_{-\infty}^{\infty}dX\tilde{\rho}(X)G(x-X),$ (18)

with one-dimensional distribution $\tilde{\rho}(x_{i})=\tilde{\rho}(X)$, and
where $A$ is the number of nucleons. Principally, the distribution of the
centroid of a Gaussian, i.e., $\tilde{\rho}(X)$, should be obtained by doing
the inverse Weierstrass transformation. However, it cannot be obtained
analytically. We approximate the solution by using the form of
$\tilde{\rho}(X)=\rho_{0}\Theta(R^{WT}-X)$, where $\Theta(R^{WT}-X)$ is a
Heaviside step function. By best fitting
$\rho^{WT}(x)=\int\tilde{\rho}(X)G(x-X)dX$ to the required density profile
$\rho^{ws}(x)$, the parameter $R^{WT}$ and the $\sigma_{r}$ in $G(x-X)$ are
obtained.

In Fig. 1 (a), a Woods-Saxon density distribution with $R=7$ fm and $a=0.5$ fm
is plotted as the red line. The black dashed line and dash-dotted line in
panel (a) are the results obtained with different $\sigma_{r}$ at $R^{WT}=R=7$
fm. The results with $\sigma_{r}\approx 0.87$ fm are close to the Woods-Saxon
distribution. Figure 1 (b) shows the $\chi^{2}$ of fitting $\rho^{WT}$ to
$\rho^{ws}$ on $R^{WT}$ and the $\sigma_{r}$ plane. At $\chi^{2}=2.5\times
10^{-5}$, we obtain that $R^{WT}=R$, $\sigma_{r}=0.87$ fm. More generally, for
reproducing the required Woods-Saxon distribution, two criteria should be
fulfilled: a) the radius of the hard sphere is equal to half the density
radius $R$, i.e., $R^{WT}=R$, and b)$\sigma_{r}$ is related to $a$ according
to the relationship

$\sigma_{r}=f(a)=ka+c,R\geq 4.4a,$ (19)

where $k=1.71217\pm 0.01548$ and $c=0.01564\pm 0.01047$ fm should be
fulfilled. For the Woods-Saxon distributions in three-dimensional space, the
above criteria are also approximately available.

Figure 1: (Color online) (a) Woods-Saxon density distribution with $R$=7 fm
and $a$=0.5 fm (red line), and Weierstrass transformation of
$\Theta(R^{WT}-r)$ with different $\sigma_{r}$ at $R$=7 fm (black lines); (b)
$\chi^{2}=\frac{1}{N}\sum_{i=1}^{N}{(\rho^{WT}_{i}-\rho^{WS}_{i})}^{2}$ as a
function of $\sigma_{r}$ and $R^{WT}$.

####  II.2.3 Stability of initial nucleus with different $\sigma_{r}$

In the QMD approach, the widths of the wave packets are time independent and
is the same for the proton and neutron. Thus, one cannot exactly reproduce the
density profile within the current framework of ImQMD due to the following
points. First, the diffuseness parameters of the density distribution of
sampled nuclei have $a_{n}=a_{p}$, which is not the same as in the
distribution obtained with RDV where $a_{n}\neq a_{p}$. Second, the commonly
used values of $\sigma_{r}$ in the QMD type models are larger than the values
extracted based on Eq. (19). For example, $\sigma_{r}=1.414$ fm is usually
used in the QMD calculations for Au+Au, but in the case of fitting the density
distribution obtained by RDV the values of $\sigma_{r}$ for proton and neutron
are $\sigma_{r}^{p}=$0.73 fm and $\sigma_{r}^{n}=0.88$ fm, respectively. It is
found that the small $\sigma_{r}$ which can reproduce the density distribution
obtained with RDV is worse in the stability of initial nuclei.

To understand the influence of $\sigma_{r}$ on the stability of sampled
nucleus, we present the root-mean-square ($rms$) values of sampled 124Sn in
ImQMD-L as a function of time with three values of $\sigma_{r}$=1.1, 1.29, 1.4
fm in Fig. LABEL:fig2-rms-t. In these calculations, the positions of neutrons
and protons are sampled within $R_{n}$ and $R_{p}$, whose values are obtained
by the RDV method with the same Skyrme energy density functional as that used
in the mean-field propagation. To quantify the stability of the sampled
initial nucleus, a variable named the probability of stability is used, and is
defined as

$P_{stab}=\frac{N_{stab}}{N_{total}},$ (20)

where $N_{stab}$ is the number of events which keeps the $rms$ variation
within 20% at 200 fm/$c$ and $N_{tot}$ is the number of total events. For
$\sigma_{r}$=1.1 fm, in panel (a), the $rms$ values obviously increase with
time and all of them are larger than 20% of $rms$ at t=0 fm/$c$, i.e.,
$P_{stab}=0.0$. The bad stability of the initial nucleus with smaller
$\sigma_{r}$ is caused by a large initial fluctuation of the density
distribution which may produce a stronger repulsive force on the particle.
With the $\sigma_{r}$ increasing, the stability of the initial nucleus becomes
better. $P_{stab}$ reaches about 87% with $\sigma_{r}$=1.29 fm, and reaches
about 100% with $\sigma_{r}$=1.4 fm.

Consequently, the stability of the initial nucleus and reproducing the initial
density profile cannot be achieved simultaneously in a QMD type model. To
balance the accuracy of the reproduction of initial density distributions and
stability of the initial nucleus, we select the width of the wave packet as
$\sigma_{r}=1.287$ fm in the following 124Sn+112Sn calculations.

####  II.2.4 Initialization in ImQMD-L with $R_{n}$ and $R_{p}$

In the initialization used in this work, the centroids of the wave packet for
neutrons and protons are sampled within the half-density radii $R_{n}$ and
$R_{p}$, and the binding energy of the sampled nucleus is in the range of
$B\pm 0.2$ MeV. The values of $B$, $R_{n}$, and $R_{p}$ are calculated based
on the RDV with the same Skyrme energy density functional as that used in the
mean-field propagation of the ImQMD-L model. In Table 2, we present the
$R_{n}$, $R_{p}$, and the binding energy $B$ of 124Sn obtained based on RDV
with five parameter sets characterized by $L$=30, 50, 70, 90, 110 MeV,
respectively. Thus, the neutron skin effect is correlated to the energy
density functional used in the mean-field potential in ImQMD-L. It is
different from the method used in many other QMD codes, in which the $R_{n/p}$
is obtained by using the formula $r_{0}A^{1/3}$[42].

For 124Sn+112Sn, $\sigma_{r}=1.287$ fm is used and the sampled density
distribution has a corresponding diffuseness with
$a_{n}=a_{p}=a=f^{-1}(\sigma_{r}=1.287)=0.743$ fm. If there is no
specification, $\sigma_{r}=1.287$ fm is the default value in the calculations.
This $a$ value is larger than that obtained with the RDV method. Thus, the
binding energy $B$ of the sampled nucleus deviates from the ground state
energy. How much deviation of the binding energy is caused by the sampled
density distribution with $a=0.743$ fm? In Table 2, the values in the brackets
of the sixth column are the binding energies of 124Sn obtained with the RDV
method for the Woods-Saxon density distribution with $R_{n}$, $R_{p}$, and
$a=a_{n}=a_{p}=0.743$ fm. The binding energy per nucleon of 124Sn in the
initialization with $a=0.743$ fm (or $\sigma_{r}=1.287$ fm) is about $\approx
0.65$ MeV larger than its ground state energy obtained with RDV, i.e., the
binding energy deviation is about 8.2%. It means that there is a spurious
excitation in the QMD initialization. ion.

##  III Results and discussions

In the following studies, we perform the calculation of ${}^{124}Sn+^{112}Sn$
at the beam energy of 200 MeV/u, and impact parameter b=2 fm. Two isospin
sensitive HIC observables will be discussed in this section. One is the
neutron to proton yield ratio, i.e., $R(n/p)$, for emitted free nucleons,

$R(n/p)=\frac{dY_{n}}{dE_{k}}/\frac{dY_{p}}{dE_{k}}.$ (21)

This observable is sensitive to the strength of the symmetry potential[10,
16], because the symmetry potential is opposite in sign for protons and
neutrons. Another observable is the coalescence invariant neutron to proton
yield ratio proposed in Ref.[16, 14], which is defined as
$R_{ci}(n/p)=\frac{dY^{CI}_{n}}{d(E_{k}/A)}/\frac{dY^{CI}_{p}}{d(E_{k}/A)}$,
with

$\displaystyle\frac{dY^{CI}_{n}}{d(E_{k}/A)}$ $\displaystyle=$
$\displaystyle\sum_{A\leq 16,Z\leq 6}\frac{dY(Z,N)}{d(E_{k}/A)}\times N,$ (22)
$\displaystyle\frac{dY^{CI}_{p}}{d(E_{k}/A)}$ $\displaystyle=$
$\displaystyle\sum_{A\leq 16,Z\leq 6}\frac{dY(Z,N)}{d(E_{k}/A)}\times Z.$

The $R_{ci}(n/p)$ ratio still retain sensitivity to the symmetry energy, and
it could eliminate the problem related to the absolute yield of light charged
particles in the transport model simulations. Both $R(n/p)$ and $R_{ci}(n/p)$
are constructed from the transverse emitted nucleons and particles with an
angular cut of $70^{\circ}\leq\theta_{c.m.}\leq 110^{\circ}$.

In practical calculations, the exact values of $R(n/p)$ and $R_{ci}(n/p)$ also
depend on the treatment of the three-body force in the mean-field potential.
In this section, we reexamine the sensitivity of isospin sensitive observables
to the slope of the symmetry energy ($L$) with the new treatment of the
initialization and mean-field potential. In addition, the influences of the
physical variables, such as $S_{0}$, $f_{I}$, $m_{s}^{\ast}$, and model
parameter $\sigma_{r}$, are also discussed.

###  III.1 Influences of the treatment of three-body force term and $L$ on the
n/p ratios

Figure LABEL:fig3-np-ek (a) shows the results of $R(n/p)$ obtained with ImQMD
[gray lines, using Eq.(7) to calculate the three-body force term], and ImQMD-L
[color shaded region, using Eq.(8) to calculate the three-body force]. Three
different symmetry energy cases, i.e., $L$=30 MeV (gray solid lines, blue
shaded region), 70 MeV (gray dashed lines, green shaded region), 110 MeV (gray
dash-dotted lines, red shaded region), are presented.

As shown in panel (a), the $R(n/p)$ values of emitted free nucleons with
$E_{k}<$80 MeV are greater for the soft symmetry energy case than for the
stiff symmetry energy case, and $R(n/p)_{L=30}>R(n/p)_{L=70}>R(n/p)_{L=110}$.
This is because the emitted nucleons with lower kinetic energy mainly come
from the overlap region during the expansion phase. This region is below
saturation density, where the symmetry energy is larger for smaller $L$.
Larger symmetry energy in this region results in enhanced neutron emissions
from the neutron-rich system and, thus, larger values of the $R(n/p)$ ratio.
This behavior has been observed in other transport model predictions,
especially for the widely studied Sn+Sn at the beam energy of 50 MeV/u[56].

At $E_{k}>120$ MeV, there is $R(n/p)_{L=30}<R(n/p)_{L=70}<R(n/p)_{L=110}$. It
is different than we observed at the beam energy around 50 MeV/u. This is
because the emitted nucleons with high kinetic energy mainly come from the
overlap region at the early stage of the expansion. During this stage, the
density in the overlapped region is above the saturation density, where the
stiff symmetry energy has larger values than the soft symmetry energy. It
results in the $R(n/p)$ values being greater in the stiff symmetry energy case
than in the soft symmetry energy case. The important finding is that the
sensitivity of $R_{n/p}$ to $L$ becomes clearer within the calculations of
ImQMD-L compared with the calculations of ImQMD. The reason is that the
calculations in ImQMD-L provide stronger strength of the three-body force at
high density, and thus of the strength of the symmetry potential.

In Fig. LABEL:fig3-np-ek (b), the $R_{ci}(n/p)$ are presented. The sensitivity
of $R_{ci}(n/p)$ to $L$ at low kinetic energy vanishes due to the
contributions from the light particles and clusters. At high kinetic energy,
the sensitivity of $R_{ci}(n/p)$ to $L$ is retained in the calculations with
ImQMD-L.

###  III.2 Influences of $S_{0}$, $f_{I}$ and $m_{s}^{*}$ on the n/p ratios

As mentioned in Ref.[9], the strength of the symmetry energy not only depends
on $L$ but also depends on the symmetry energy coefficient $S_{0}$, isoscalar
effective mass $m^{\ast}_{s}$, and $f_{I}$. However, the parameters $S_{0}$,
$m_{s}^{\ast}$ and $f_{I}$ also have some uncertainties[9, 49], and it is
worthwhile to understand their influence on $R(n/p)$.

The results of $R(n/p)$ obtained with $S_{0}=34$ MeV are presented as solid
shaded regions in panel (a) of Fig.LABEL:fig4-np-s0fims. All the other
parameters are kept the same as in the default parameter sets. The different
colors correspond to the results obtained with $L$=30, 70, 110 MeV. When a
larger $S_{0}$ is adopted in the calculations, the curves of $R(n/p)$ shift up
by $<$8% owing to the enhanced symmetry energy in all density regions.
However, the sensitivity of $R(n/p)$ to $L$ is not dramatically changed.

Another important ingredient that can influence $R(n/p)$ is the neutron-proton
effective mass splitting, i.e., $m_{n}^{\ast}>m_{p}^{\ast}$ or
$m_{n}^{\ast}<m_{p}^{\ast}$[57]. In panel (b) of the Fig.LABEL:fig4-np-s0fims,
the solid shaded regions are the results obtained with
$m_{n}^{\ast}<m_{p}^{\ast}$, i.e., $f_{I}=0.178$, and different colors
represent different $L$. For the given sets with $m_{n}^{\ast}<m_{p}^{\ast}$,
i.e., $f_{I}=0.178$, the sensitivity of $R(n/p)$ to $L$ at lower kinetic
energy region is similar to that obtained with default parameter sets where
$f_{I}=-0.178$ ($m_{n}^{\ast}>m_{p}^{\ast}$). Thus, the $R(n/p)$ at $E_{k}\leq
80$ MeV can help us to determine the values of $L$.

In the high kinetic energy region, the sensitivity of $R(n/p)$ to $L$ becomes
complicated if the neutron-proton effective mass splitting is not well fixed.
For example, the values of $R(n/p)$ at $E_{k}>100$ MeV obtained with
$f_{I}$=0.178 and $L$=30 MeV are close to that obtained with $f_{I}$=-0.178
and $L$=70 MeV. This is because the sets with $f_{I}$=0.178
($m_{n}^{\ast}<m_{p}^{\ast}$) have the strong lane potential as described in
Ref. [57], and it enhances the values of $R(n/p)$. To distinguish it, the
inclination of $R(n/p)$ to $E_{k}$, which is calculated with the values of
$R(n/p)$ obtained at $E^{(1)}_{k}$=50 MeV and $E^{(2)}_{k}$=150 MeV, i.e.,

$RE_{k}=\frac{\Delta R(n/p)}{\Delta E_{k}},$ (23)

where $\Delta R(n/p)=R(n/p)(E^{(1)}_{k})-R(n/p)(E^{(2)}_{k})$ and $\Delta
E_{k}=E^{(1)}_{k}-E^{(2)}_{k}$, can be adopted. The results of $\Delta
R/\Delta E_{k}$ for different $f_{I}$ at $K_{0}$=240 MeV, $S_{0}$=30 MeV, and
$m_{s}^{*}=0.8m$ are presented Fig. LABEL:fig4-ms-fis, and the calculations
with larger $f_{I}$ predict larger inclination.

Based on the above discussions, one can expect that the spectra of $R(n/p)$
can be used to obtain information on $L$ and $f_{I}$ simultaneously, if one
does the analysis on the $L$ and $f_{I}$ parameter space.

The influence of $m_{s}^{\ast}$ on $R(n/p)$ is also investigated in the case
where other parameters are kept the same as in Table 1. The values of $R(n/p)$
obtained with $m_{s}^{\ast}=0.65m$ are presented as solid shaded region in
panel (c) of Fig. LABEL:fig4-np-s0fims. Different colors are the results
obtained with different $L$. For $L$=70 and 110 MeV, reducing the isoscalar
effective mass $m_{s}^{\ast}/m$ to 0.65 changes the $R(n/p)$ weakly. For
$L$=30 MeV, obvious differences of $R(n/p)$ are observed between the set with
$m_{s}^{\ast}=0.65m$ and the default parameter set. The $R(n/p)$ obtained with
$m_{s}^{\ast}/m=0.65$ are largely enhanced in the lower kinetic energy region
and are obviously suppressed in the high kinetic energy region. This is
because the set with $L$=30 MeV and $m_{s}^{\ast}/m=0.65$ have strong symmetry
energy at low density and strong momentum dependent interaction, both
providing a stronger repulsive force than other parameter sets and worsening
the stability of the initial nucleus in the ImQMD-L simulations.

###  III.3 Influence of $\sigma_{r}$ on the n/p ratios

The width of wave packet $\sigma_{r}$ is an important model parameter in the
QMD type models, and one can expect that the different values of $\sigma_{r}$
may influence the results. Thus, investigating the influence of $\sigma_{r}$
on the $R(n/p)$ within the ImQMD-L model can help us to understand the
robustness of model calculations. In this work, the test calculations are
performed with the value of $\sigma_{r}=1.0$ fm, which is not really used in
low-intermediate energy HIC simulations.

Figure 6 (a) shows the $R(n/p)$ obtained with different $L$ at
$\sigma_{r}=1.0$ fm. In the calculations with $\sigma_{r}=1.0$ fm, the values
of $R(n/p)$ still depend on the slope of symmetry energy, and the behavior is
similar to that with $\sigma_{r}=1.287$ fm. However, the sensitivity becomes
weaker in both low and high kinetic energy regions than that with
$\sigma_{r}=1.287$ fm. The reason is that the initial nuclei obtained with
$\sigma_{r}=1.0$ fm in the ImQMD-L model have stronger initial fluctuation and
worse stability than that with $\sigma_{r}=1.287$ fm. It increases spurious
emission of nucleons with low kinetic energy. For example, for $L=30$ MeV, 22%
(18%) more protons (neutrons) are emitted for the case of $\sigma_{r}=1.0$ fm
than that with $\sigma_{r}=1.287$ fm. Consequently, the $R(n/p)$ obtained with
$\sigma_{r}=1.0$ fm is smaller than that with $\sigma_{r}=1.287$ fm at lower
kinetic energy. For $L=110$ MeV, 6% (6.7%) more protons (neutrons) are emitted
for the case with $\sigma_{r}=1.0$ fm than that with $\sigma_{r}=1.287$ fm,
and thus the $R(n/p)$ obtained with $\sigma_{r}=1.0$ fm is enhanced in the low
kinetic energy region. Due to the conservation of nucleon number in the
reaction system, $R(n/p)$ obtained with $\sigma_{r}=1.0$ fm is enhanced at
high kinetic energy region for $L=30$ MeV, and suppressed for $L=110$ MeV.
Thus, the effects of symmetry energy become relatively weak in the case of
$\sigma_{r}=1.0$ fm. However, the sensitivity of $R(n/p)$ to $L$ is still
large enough to distinguish the stiffness of symmetry energy.

Figure 6: (Color online) Panels (a)-(d) are the neutron to proton yield ratios
for emitted free nucleons obtained with different $L$ and different $f_{I}$ at
$\sigma_{r}=1.0$ fm.

Figure 6(b), 6(c) and 6(d) show the results of $R(n/p)$ obtained with
$f_{I}=-0.178$ ($m_{n}^{*}>m^{*}_{p}$) and $f_{I}=0.178$
($m_{n}^{*}<m^{*}_{p}$) at $L$=30, 70, 110 MeV, respectively. As we discussed
previously, the values of $R(n/p)$ in the high kinetic energy region obtained
with $f_{I}=0.178$ ($m_{n}^{*}<m^{*}_{p}$) are greater than that obtained with
$f_{I}=-0.178$ ($m_{n}^{*}>m^{*}_{p}$) at given $L$. The finding is similar to
the results obtained with $\sigma_{r}=1.287$ fm, and it demonstrates that the
small value of $\sigma_{r}=1$ fm do not dramatically change the sensitivity of
$R(n/p)$ to the effective mass splitting.

##  IV Summary and outlook

In summary, we have investigated the issues of how to reproduce and how well
to reproduce the Woods-Saxon density distribution with a Gaussian wave
function. Based on the inverse Weierstrass transformation, we obtain two
criteria for reproducing the Woods-Saxon density profile: one is that the
centroid of the Gaussian wave packet should be sampled within a hard sphere
with a radius approximately equal to the half density radius of the Woods-
Saxon density distribution. Another is that the width of the Gaussian wave
packet should be taken based on the parameter $a$, i.e.
$\sigma_{r}=0.01564+1.71217a$. The second condition requires a smaller width
of Gaussian wave packet than the commonly used values in the quantum molecular
dynamics model, and it causes initial nuclei to have bad stability due to the
strong initial fluctuation. To balance the stability of initial nuclei and the
requirement of initial density profile, the width of the Gaussian wave packet
is taken as the one that can give the probability of stability greater than
85% in this work.

To correlate the initialization and mean-field potential with the same energy
density functional in the improved quantum molecular dynamics, we incorporate
the RDV method into the ImQMD model. First, we calculate the half density
radius of $R_{n}$, $R_{p}$, and the binding energy $B$ of the initial nuclei
with the RDV method by using the same energy density functional as in the
mean-field propagation. Then, the initial nuclei are sampled by using the
obtained values of $R_{n}$, $R_{p}$, and binding energy $B$. In the mean-field
part, we precisely calculate the three-body related term by using the 11-point
Gauss-Legendre quadrature method at the expense of 50 times longer CPU times.
These methods improve the theoretical reliability of the transport model and
reduce the uncertainties of theoretical predictions owing to the consideration
of the correlation between the initialization and the mean-field potential.

Based on the new version of the model in this work, we study the influence of
the slope of symmetry energy, effective mass splitting on the heavy ion
collision observables, such as the neutron to proton yield ratios for emitted
free nucleons and for coalescence invariant nucleon yield, for 124Sn+112Sn at
the beam energy of 200 MeV per nucleon. Our calculations show that the
$R(n/p)$ at low kinetic energy region is sensitive to the slope of symmetry
energy, and the $R(n/p)$ at high kinetic energy depends on the stiffness of
symmetry energy and neutron-proton effective mass splitting. With a given $L$,
the inclination of $R(n/p)$ to kinetic energy ($E_{k}$) can be used to probe
the effective mass splitting. On the other hand, if the neutron-proton
effective mass splitting is fixed, the values of $R(n/p)$ in the high kinetic
energy region increase with $L$, which also implies that the $R(n/p)$ in the
high kinetic energy region can be used to probe the symmetry energy above the
saturation density. Varying the width of wave packet parameter $\sigma_{r}$ in
a reasonable region does not dramatically change our conclusions. The current
calculations show that the analysis of spectra of $R(n/p)$ on the $L$ and
$f_{I}$ parameter space is necessary in future to tightly constrain the
symmetry energy.

## Acknowledgements

This work was partly inspired by the transport code comparison project, and it
was supported by the National Natural Science Foundation of China No.
11875323, 11705163, 11790320, 11790323, U1867212, 114222548, and 11961141003,
the National Key R&D Program of China under Grant No. 2018YFA0404404, the
Continuous Basic Scientific Research Project (No. WDJC-2019-13, BJ20002501),
Innovation Project of Guangxi Graduate Education No. XYCSZ2018060, and the
funding of China Institute of Atomic Energy. The work was carried out at
National Supercomputer Center in Tianjin, and the calculations were performed
on TianHe-1 (A).

##

## References

  * [1] P. Moller, W.D. Myers, H. Sagawa, and S. Yoshida, Phys. Rev. Lett. 108, 052501 (2012).
  * [2] P. Danielewicz and J. Lee, Nucl. Phys. A 922, 1 (2014) .
  * [3] M. Centelles, X. Roca-Maza, X. Vinas, M. Warda, Phys. Rev. Lett. 102, 122502 (2009).
  * [4] B.A. Brown, Phys. Rev. Lett. 85, 5296 (2000).
  * [5] L.W. Chen, C.M. Ko, B.A. Li, J. Xu, Phys. Rev. C 82, 024321 (2010).
  * [6] M Liu, N Wang, Z.-X. Li, and F.-S. Zhang, Phys. Rev. C 82, 064306 (2010).
  * [7] H. Jiang, G. J. Fu, Y. M. Zhao, and A. Arima, Phys. Rev. C 85, 024301 (2012).
  * [8] N. Wang, L. Ou, and M. Liu, Phys. Rev. C 87, 034327 (2013).
  * [9] Y. Zhang, M. Liu, C.-J. Xia, Z. Li, and S.K. Biswal, Phys. Rev. C 101, 034303 (2020).
  * [10] B.A. Li, C.M. Ko, and Z. Ren, Phys. Rev. Lett. 78, 1644 (1997).
  * [11] M.B. Tsang, T.X. Liu, L.Shi, P. Danielewicz, C.K. Gelbke, X.D. Liu, W.G. Lynch, W.P. Tan, G. Verde et al., Phys. Rev. Lett. 92, 062701 (2004).
  * [12] L.W.Chen, C.M. Ko, B.A. Li, Phys. Rev. Lett. 94, 032701 (2005).
  * [13] V. Baran, M. Colonna, V. Greco, and M. Di Toro, Phys. Rep. 410, 335 (2005).
  * [14] M.A. Famiano, T. Liu, W.G. Lynch, M. Mocko, A.M. Rogers, M.B. Tsang et al., Phys. Rev. Lett. 97, 052701 (2006).
  * [15] Y. Zhang, P. Danielewicz, M. Famiano, Zhuxia Li, W.G. Lynch, M.B. Tsang, Phys. Lett. B 664, 145 (2008).
  * [16] M.B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W.G. Lynch, and A. W. Steiner, Phys. Rev. Lett. 102, 122701 (2009).
  * [17] B.A. Li, L.W. Chen, and C.M. Ko, Phys. Rep. 464, 113 (2008).
  * [18] M.B. Tsang, J.R. Stone, F. Camera, P. Danielewicz, S. Gandolfi, K. Hebeler et al., Phys. Rev. C 86, 015803 (2012).
  * [19] B.A. Li, A. Ramos, G. Verde et al. (editors), Eur. Phys. J. A 50, 9 (2014), Topical Issue on Nuclear Symmetry Energy.
  * [20] Q. Wu, Y. Zhang, Z. Xiao, R. Wang, Y. Zhang, Z. Li, N. Wang, and R.H. Showalter, Phys. Rev. C 91, 014617 (2015).
  * [21] B.P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * [22] B.P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018).
  * [23] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018).
  * [24] F.J. Fattoyev, J. Piekarewicz, and C.J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * [25] B.P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 9, 011001 (2019).
  * [26] T. Malik, B.K. Agrawal, J.N. De, S.K. Samaddar, C. Providencia, C. Mondal, and T.K. Jha, Phys. Rev. C 99, 052801(R) (2019).
  * [27] N.B. Zhang and B.A. Li, J. Phys. G: Nucl. Part. Phys. 46, 014002 (2019).
  * [28] C.Y. Tsang, M.B. Tsang, P. Danielewicz, F.J. Fattoyev, W.G. Lynch, Phys. Lett. B 796, 1 (2019).
  * [29] M.B. Tsang, W.G. Lynch, P. Danielewicz, C.Y. Tsang, Phys. Lett. B 795, 533 (2019).
  * [30] Q. Li, Z. Li, S. Soff, M. Bleicher, and H. Stocker, Phys. Rev. C.72, 034613 (2005).
  * [31] G.-C. Yong, B.-A. Li, L.-W. Chen, and X.-C. Zhang, Phys. Rev. C 80, 044608 (2009).
  * [32] Y. Zhang, J. Tian, W. Cheng, F. Guan, Y. Huang et al., Phys. Rev.C 95, 041602(R) (2017).
  * [33] B.A. Li, Phys. Rev. Lett. 88, 192701 (2002).
  * [34] Y. Liu, Y. Wang, Y. Cui, C.-J. Xia, Z. Li, Y. Chen, Q. Li, and Y. Zhang, Phys. Rev. C 103, 014616 (2021).
  * [35] P. Russotto, S. Gannon, S. Kupny, P. Lasko, L. Acosta et. al., Phys. Rev. C 94, 034608 (2016).
  * [36] Y. J. Wang, Q. Li, Y. Leifels, and A. Fevre, Phys. Lett. B 802, 135249 (2020).
  * [37] C. Hartnack, A. Le Fevre, and Y. Leifels, J. Aichelin, arXiv:1808.09868.
  * [38] G.-F. Wei, B.A. Li, J. Xu, L.W. Chen, Phys. Rev. C 90, 014610 (2014).
  * [39] S.R. Souza, L.de Paula, S. Leray, J. Nemeth, C. Ngo, H. Ngo, Nuclear Physics A 571, 159-184 (1994).
  * [40] K. Niita, S. Chiba, T. Maruyama, T. Maruyama, H. Takada, T. Fukahori, Y. Nakahara, A. Iwamoto, Phys. Rev. C. 52, 2620 (1995).
  * [41] T. Maruyama, A. Ohnishi, and H. Horiuchi, Phys. Rev. C 42, 386 (1990).
  * [42] J. Xu, L.-W. Chen, ManYee B. Tsang, H. Wolter, Y.-X. Zhang, J. Aichelin,M. Colonna, D. Cozma, P. Danielewicz, Z.-Q. Feng, A.L. Fevre, T. Gaitanos, C. Hartnack, K. Kim, Y. Kim, C.-M. Ko, B.-A. Li, Q.-F. Li, Z.-X. Li, P. Napolitani et al. Phys. Rev. C. 93, 044609 (2016).
  * [43] S.A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L.Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann, J.Konopka et al., Prog. Part. Nucl. Phys. 41, 255 (1998).
  * [44] J. Aichelin, E. Bratkovskaya, A. Le Fevre, V. Kireyeu, V. Kolesnikov, Y. Leifels and V. Voronyuk, and G. Coci, Phys. Rev. C 101, 044905 (2020).
  * [45] M.D. Cozma, Eur. Phys. J. A 54, 40 (2018).
  * [46] Y. Zhang, N. Wang, Q. Li, L. Ou, J. Tian, M. Liu, K. Zhao, X. Wu, Z. Li, Front. Phys. 15, 54301 (2020).
  * [47] T.H.R. Skyrme, Philos. Mag. 1 ,1043 (1956).
  * [48] Y. Zhang and Z. Li, Phys. Rev. C 74, 014602 (2006).
  * [49] M. Dutra, O. Lourenco, J. S. SaMartins, A. Delfino, J. R. Stone, and P. D. Stevenson, Phys. Rev. C 85, 035201 (2012).
  * [50] Z. Q. Feng, Phys. Rev. C 84, 024610 (2011).
  * [51] M. D. Cozma, Y. Leifels, W. Trautmann, Q. Li, and P. Russotto, Phys. Rev. C 88, 044912 (2013).
  * [52] J. Su, F. S. Zhang, and B. A. Bian, Phys. Rev. C 83, 014608 (2011).
  * [53] M. Liu, N. Wang, Z.-X. Li, and X.-Z. Wu, Chin. Phys. Lett. 23, 804 (2006).
  * [54] M. Brack, C. Guet, and H.B. Hakanson, Phys. Rep. 123, 275 (1985).
  * [55] G.G. Bilodeau, Duke Math. J. 29, 293 (1962).
  * [56] H.-Y. Kong, Y. Xia, J. Xu, L.-W. Chen, B.-A. Li, and Y.-G. Ma, Phys. Rev. C 91, 047601 (2015).
  * [57] Y. Zhang, M.B. Tsang, Z. Li, and H. Liu, Phys. Lett. B 732, 186 (2014).