Nuclear collective dynamics in the lattice Hamiltonian Vlasov method

TL;DR

This paper introduces a new lattice Hamiltonian method for solving the Vlasov equation in nuclear physics, enabling stable simulation of nuclear ground states and collective excitations with results consistent with experiments.

Contribution

The paper develops a novel lattice Hamiltonian Vlasov approach incorporating Skyrme pseudopotential up to N3LO, improving stability and accuracy in nuclear dynamic simulations.

Findings

Stable nuclear ground state evolution demonstrated.

Results for giant monopole and dipole modes agree with RPA and experiments.

Method effectively models long-time nuclear dynamics.

Abstract

The lattice Hamiltonian method is developed for solving the Vlasov equation with nuclear mean-field based on the Skyrme pseudopotential up to next-to-next-to-next-to leading order. The ground states of nuclei are obtained through varying the total energy with respect to the density distribution of nucleons. Owing to the self-consistent treatment of initial nuclear ground state and the exact energy conservation in the lattice Hamiltonian method, the present framework of solving the Vlasov equation exhibits very stable nuclear ground state evolution. As a first application of the new lattice Hamiltonian Vlasov method, we explore the iso-scalar giant monopole and iso-vector giant dipole modes of finite nuclei. The obtained results are shown to be comparable to that from random-phase approximation and consistent with the experimental data, indicating the capability of the present method in dealing with the long-time near-equilibrium nuclear dynamics.

Tables

Table 1. Table 1: Parameters related to the nuclear matter properties of three N 3 3 3 LO Skyrme pseudopotentials, SP 6 6 6 s, SP 6 6 6 m, and SP 6 6 6 h, and one conventional Skyrme interaction MSL 1 1 1 , where the recombination of Skyrme parameters defined in Eqs. ( 11 ) and ( 12 ) are used.

	SP$6$s	SP$6$m	SP$6$h	MSL$1$
$t_0$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^3$)	-1814.64	-1956.75	-1675.52	1963.23
$x_0$	0.5400	0.2306	-0.0902	0.3208
$t_3$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^{3+3\alpha }$)	10796.2	11402.9	9873.1	12174.9
$x_3$	0.8257	0.1996	-0.4990	0.3219
$\alpha$	0.2923	0.2523	0.3168	0.2694
$C^{\left[2\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^5$)	597.877	637.195	677.884	435.519
$D^{\left[2\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^5$)	-446.695	-524.373	-601.990	-367.583
$C^{\left[4\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^7$)	-26.2027	-28.5209	-31.2026	0.0
$D^{\left[4\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^7$)	23.2525	27.6873	32.4607	0.0
$C^{\left[6\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^9$)	0.0903	0.1000	0.1121	0.0
$D^{\left[6\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^9$)	-0.0896	-0.1080	-0.1292	0.0

Table 2. Table 2: Macroscopic characteristic quantities for SP 6 6 6 s, SP 6 6 6 m, SP 6 6 6 h, and MSL 1 1 1 , where ρ sc subscript 𝜌 sc \rho_{\rm sc} = = 0.11 / 0.16 ​ ρ 0 0.11 0.16 subscript 𝜌 0 0.11/0.16\rho_{0} and ρ h subscript 𝜌 h \rho_{\rm h} = = 0.5 ​ fm − 3 0.5 superscript fm 3 0.5~{}{\rm fm}^{-3} .

	SP$6$s	SP$6$m	SP$6$h	MSL$1$
${\rho }_0$ (${\mathrm{fm}}^{-3}$)	$0.1614$	$0.1630$	$0.1647$	$0.1586$
$E_0$ ($\mathrm{MeV}$)	$-16.04$	$-15.94$	$-15.61$	$-16.00$
$K_0$ ($\mathrm{MeV}$)	$240.9$	$233.4$	$240.8$	$235.1$
$J_0$ ($\mathrm{MeV}$)	$-377.0$	$-384.2$	$-358.2$	$-372.7$
$E_{\mathrm{sym}}\left({\rho }_{\mathrm{sc}}\right)$ ($\mathrm{MeV}$)	$25.43$	$25.83$	$25.98$	$26.67$
$L\left({\rho }_{\mathrm{sc}}\right)$ ($\mathrm{MeV}$)	$32.47$	$46.75$	$62.19$	$46.19$
$E_{\mathrm{sym}}\left({\rho }_0\right)$ ($\mathrm{MeV}$)	$28.84$	$31.93$	$34.97$	$32.33$
$L\left({\rho }_0\right)$ ($\mathrm{MeV}$)	$18.20$	$49.10$	$82.17$	$45.25$
$K_{\mathrm{sym}}$ ($\mathrm{MeV}$)	$-242.7$	$-158.0$	$-70.5$	$-183.3$
$E_{\mathrm{sym}}\left(2{\rho }_0\right)$ ($\mathrm{MeV}$)	$24.06$	$41.31$	$61.62$	$39.00$
$E_{\mathrm{sym}}\left({\rho }_h\right)$ ($\mathrm{MeV}$)	$0.03$	$41.32$	$79.82$	$31.01$
$m_{s,0}^{\ast }/m$	$0.759$	$0.758$	$0.755$	$0.806$
$m_{v,0}^{\ast }/m$	$0.678$	$0.663$	$0.648$	$0.706$

Table 3. Table 3: The gradient parameters, both E [ 2 ] superscript 𝐸 delimited-[] 2 E^{[2]} and E ~ [ 2 ] superscript ~ 𝐸 delimited-[] 2 \tilde{E}^{[2]} for SP 6 6 6 s, SP 6 6 6 m, SP 6 6 6 h, and MSL 1 1 1 used in the present work. The obtained binding energy and proton rms radius of \isotope [208]Pb of TF initialization and LHV calculations are also shown.

	SP$6$s	SP$6$m	SP$6$h	MSL$1$	Exp.
$E^{\left[2\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^5$)	-250.0	-200.0	-150.0	-250.0	-
${\tilde{E}}^{\left[2\right]}$ ($\mathrm{MeV}\cdot {\mathrm{fm}}^5$)	-15.0	-10.0	-10.0	-20.0	-
BE ($\mathrm{MeV}$)	$1637.2$	$1669.7$	$1654.5$	$1632.7$	$1636.4$
$\sqrt{{⟨r^2⟩}_p}$ ($\mathrm{fm}$)	$5.48$	$5.44$	$5.40$	$5.51$	$5.45$
$\mathrm{BE}$ ($\mathrm{MeV}$) in LHV	$1557.2$	$1585.1$	$1565.1$	$1553.5$	-
$\sqrt{{⟨r^2⟩}_p}$ ($\mathrm{fm}$) in LHV	$5.52$	$5.49$	$5.44$	$5.56$	-

Equations

$$v_{Sk}=V_{\text{N3LO}}^C+V_{\text{N1LO}}^{DD},$$

$$V_{\text{N3LO}}^C=$$

$$\displaystyle+t^{[4]}_{2}(1+x^{[4]}_{2}\hat{P}_{\sigma})(\hat{\vec{k}}^{\prime}\cdot\hat{\vec{k}})(\hat{\vec{k}}^{\prime 2}+\hat{\vec{k}}^{2})+t^{[6]}_{1}(1+x^{[6]}_{1}\hat{P}_{\sigma})(\hat{\vec{k}}^{\prime 2}+\hat{\vec{k}}^{2})\Big{[}\frac{1}{2}(\hat{\vec{k}}^{\prime 2}+\hat{\vec{k}}^{2})^{2}+6(\hat{\vec{k}}^{\prime}\cdot\hat{\vec{k}})^{2}\Big{]}$$

$$+t_2^{[6]}(1+x_2^{[6]}{\hat{P}}_{\sigma })({\hat{\vec{k}}}^{\prime }\cdot \hat{\vec{k}})[3({\hat{\vec{k}}}^{\prime 2}+{\hat{\vec{k}}}^2{)}^2+4({\hat{\vec{k}}}^{\prime }\cdot \hat{\vec{k}}{)}^2],$$

$$V^{DD}_{\text{N1LO}}=\frac{1}{6}t_{3}(1+x_{3}\hat{P}_{\sigma})\rho^{\alpha}\Big{(}\frac{\vec{r}_{1}+\vec{r}_{2}}{2}\Big{)}.$$

$$\begin{array}{rl}\mathcal{H}(\vec{r})& ={\mathcal{H}}^{\mathrm{kin}}(\vec{r})+{\mathcal{H}}^{\mathrm{loc}}(\vec{r})\\ & +{\mathcal{H}}^{\mathrm{MD}}(\vec{r})+{\mathcal{H}}^{\mathrm{grad}}(\vec{r})+{\mathcal{H}}^{\mathrm{DD}}(\vec{r}),\end{array}$$

$${\mathcal{H}}^{\mathrm{kin}}(\vec{r})=\sum _{\tau =n,p}\int d^{3}p\frac{p^2}{2m_{\tau }}{f}_{\tau }(\vec{r},\vec{p})$$

$${\cal H}^{\rm loc}(\vec{r})=\frac{t_{0}}{4}\bigg{[}(2+x_{0})\rho(\vec{r})^{2}-(2x_{0}+1)\sum_{\tau=n,p}\rho_{\tau}(\vec{r})^{2}\bigg{]}$$

$$\begin{array}{rl}{\mathcal{H}}^{\mathrm{MD}}(\vec{r})& =\int d^{3}p{d}^3{p}^{\prime }{\mathcal{K}}_s(\vec{p},{\vec{p}}^{\prime })f(\vec{r},\vec{p})f(\vec{r},{\vec{p}}^{\prime })\\ & +\sum _{\tau =n,p}\int d^{3}p{d}^3{p}^{\prime }{\mathcal{K}}_v(\vec{p},{\vec{p}}^{\prime })f_{\tau }(\vec{r},\vec{p})f_{\tau }(\vec{r},{\vec{p}}^{\prime }),\end{array}$$

$${\mathcal{K}}_{\mathrm{s}}(\vec{p},{\vec{p}}^{\prime })$$

$$+\frac{C^{[6]}}{16{\hslash }^2}(\vec{p}-{\vec{p}}^{\prime }{)}^6,$$

$${\mathcal{K}}_{\mathrm{v}}(\vec{p},{\vec{p}}^{\prime })$$

$$+\frac{D^{[6]}}{16{\hslash }^2}(\vec{p}-{\vec{p}}^{\prime }{)}^6.$$

$${\cal H}^{\rm grad}(\vec{r})=\frac{1}{16}E^{[2]}\Big{\{}2\rho(\vec{r})\nabla^{2}\rho(\vec{r})-2\big{[}\nabla\rho(\vec{r})\big{]}^{2}\Big{\}}.$$

$$C^{[n]}$$

$$D^{[n]}$$

$$E^{[2]}=-\frac{1}{4}\big{[}t_{1}^{[2]}(2+x_{1}^{[2]})-t_{2}^{[2]}(2+x_{2}^{[2]})\big{]}.$$

$$\begin{split}{\cal H}^{\rm DD}(\vec{r})=\frac{t_{3}}{24}\bigg{[}(2+x_{3})\rho^{2}-(2x_{3}+1)\sum_{\tau=n,p}\rho_{\tau}^{2}\bigg{]}\rho^{\alpha}.\end{split}$$

$$\frac{\partial f}{\partial t}+\frac{\vec{p}}{E}{\nabla }_{\vec{r}}f+{\nabla }_{\vec{p}}U(\vec{r},\vec{p})\cdot {\nabla }_{\vec{r}}f-{\nabla }_{\vec{r}}U(\vec{r},\vec{p})\cdot {\nabla }_{\vec{p}}f=0,$$

$$f(\vec{r},\vec{p})=\frac{1}{(2\pi\hbar)^{3}}\int{\rm exp}\Big{(}-i\frac{\vec{p}}{\hbar}\cdot\vec{s}\Big{)}\rho(\vec{r}+\vec{s}/2,\vec{r}-\vec{s}/2)d\vec{s}.$$

$$\begin{array}{rl}{I}_{\mathrm{c}}& =-\int \frac{d{\vec{p}}_2}{(2\pi \hslash {)}^3}\frac{d{\vec{p}}_3}{(2\pi \hslash {)}^3}\frac{d{\vec{p}}_4}{(2\pi \hslash {)}^3}|{\mathcal{M}}_{12\to 34}{|}^2\\ & \times (2\pi {)}^4{\delta }^4(p_1+p_2-p_3-p_4)\\ & \times [f_1{f}_2(1-f_3)(1-f_4)-f_3{f}_4(1-f_1)(1-f_2)],\end{array}$$

$$f_{\tau}(\vec{r},\vec{p},t)\propto\sum_{i}\delta\big{[}\vec{r}_{i}(t)-\vec{r}\big{]}\delta\big{[}\vec{p}_{i}(t)-\vec{p}\big{]},$$

$$H=\int \mathcal{H}(\vec{r})d\vec{r}\approx l_x{l}_y{l}_z\sum _{\alpha }\mathcal{H}({\vec{r}}_{\alpha })\equiv H_L,$$

$$\tilde{f}_{\tau}(\vec{r}_{\alpha},\vec{p},t)=\frac{1}{2}\frac{(2\pi\hbar)^{3}}{N_{\rm E}}\sum_{i}^{\alpha,\tau}S\big{[}\vec{r}_{i}(t)-\vec{r}_{\alpha}\big{]}\delta\big{[}\vec{p}_{i}(t)-\vec{p}\big{]},$$

$$\tilde{\rho}_{\tau}(\vec{r}_{\alpha},t)=2\int\tilde{f}_{\tau}\frac{d\vec{p}}{(2\pi\hbar)^{3}}=\frac{1}{N_{\rm E}}\sum_{i}^{\alpha,\tau}S\big{[}\vec{r}_{i}(t)-\vec{r}_{\alpha}\big{]}.$$

$$\frac{d{\vec{r}}_i}{dt}$$

$$\frac{d{\vec{p}}_i}{dt}$$

$$\displaystyle=-N_{\rm E}l_{x}l_{y}l_{z}\sum_{\alpha\in V_{i}}\bigg{[}\sum_{\tau}^{n,p}\Big{(}\frac{\partial{\tilde{\cal H}}^{\rm loc}_{\alpha}}{\partial\tilde{\rho}_{\tau,\alpha}}+\frac{\partial{\tilde{\cal H}}^{\rm Cou}_{\alpha}}{\partial\tilde{\rho}_{\tau,\alpha}}+\frac{\partial{\tilde{\cal H}}^{\rm DD}_{\alpha}}{\partial\tilde{\rho}_{\tau,\alpha}}+\frac{\partial{\tilde{\cal H}}^{\rm grad}_{\alpha}}{\partial\tilde{\rho}_{\tau,\alpha}}+\sum_{n}(-1)^{n}\nabla^{n}\frac{\partial{\tilde{\cal H}}^{\rm grad}_{\alpha}}{\partial\nabla^{n}\tilde{\rho}_{\tau,\alpha}}\Big{)}\frac{\partial\tilde{\rho}_{\tau,\alpha}}{\partial\vec{r}_{i}}+\frac{\partial{\tilde{\cal H}}^{\rm MD}_{\alpha}}{\partial\vec{r}_{i}}\bigg{]}.$$

$${\cal H}^{\rm Cou}(\vec{r}_{\alpha})=e^{2}\rho_{p}(\vec{r}_{\alpha})\bigg{\{}\frac{1}{2}\int\frac{\rho_{p}(\vec{r}^{\prime})}{|\vec{r}_{\alpha}-\vec{r}^{\prime}|}d\vec{r}^{\prime}-\frac{3}{4}\Big{[}\frac{3\rho_{p}(\vec{r}_{\alpha})}{\pi}\Big{]}^{1/3}\bigg{\}}\approx\tilde{\cal H}^{\rm Cou}(\vec{r}_{\alpha})=e^{2}\tilde{\rho}_{p}(\vec{r}_{\alpha})\bigg{\{}\frac{1}{2}\sum_{\alpha^{\prime}\neq\alpha}\frac{\tilde{\rho}_{p}(\vec{r}_{\alpha^{\prime}})l_{x}l_{y}l_{z}}{|\vec{r}_{\alpha}-\vec{r}_{\alpha^{\prime}}|}-\frac{3}{4}\Big{[}\frac{3\tilde{\rho}_{p}(\vec{r}_{\alpha})}{\pi}\Big{]}^{1/3}\bigg{\}},$$

$$\frac{\partial {\tilde{\rho }}_{\tau ,\alpha }}{\partial {\vec{r}}_i}=\frac{\partial }{\partial {\vec{r}}_i}\sum _{{\vec{r}}_j\in V_{\alpha }}^{{\tau }_j=\tau }S({\vec{r}}_j-{\vec{r}}_{\alpha })=\{\begin{array}{ll}& \frac{\partial S({\vec{r}}_i-{\vec{r}}_{\alpha })}{\partial {\vec{r}}_i},{\tau }_i=\tau ,\\ & 0,{\tau }_i\ne \tau .\end{array}$$

$$\frac{\partial {\tilde{\mathcal{H}}}^{\mathrm{MD}}({\vec{r}}_{\alpha })}{\partial {\vec{r}}_i}$$

$$\frac{\partial {\tilde{\mathcal{H}}}^{\mathrm{MD}}({\vec{r}}_{\alpha })}{\partial {\vec{p}}_i}$$

$$S({\vec{r}}_i-\vec{r})$$

$$g(q)$$

$$\sum _{\alpha }\tilde{\rho }({\vec{r}}_{\alpha })l_x{l}_y{l}_z=\frac{1}{N_{\mathrm{E}}}\sum _{\alpha }\sum _{i}S({\vec{r}}_i-{\vec{r}}_{\alpha })l_x{l}_y{l}_z=A.$$

Full text

I Introduction

One-body transport models like the Boltzmann-Uehling-Uhlenbeck (BUU) equation (see, e.g., Ref. Bertsch and Das Gupta (1988)) provide a powerful tool to describe heavy-ion collisions. Based on these models, valuable information on the details of nuclear reaction dynamics and nuclear matter properties has been obtained by analyzing data in heavy-ion collisions Li et al. (1998); Danielewicz et al. (2002); Baran et al. (2005); Li et al. (2008). Among the obtained information, those concerning the equation of state (EOS) of asymmetric nuclear matter at both subsaturation and suprasaturation densities and the in-medium nuclear effective interaction are of fundamental importance, since they are crucial in investigating the properties of various nuclear systems or astrophysical objects Li et al. (1998); Danielewicz et al. (2002); Lattimer and Prakash (2004); Baran et al. (2005); Steiner et al. (2005); Lattimer and Prakash (2007); Li et al. (2008); Trautmann and Wolter (2012); Horowitz et al. (2014); Li et al. (2014); Hebeler et al. (2015); Baldo and Burgio (2016); Oertel et al. (2017); Li (2017); Li et al. (2018). Heavy-ion collisions not only serve as an alternative method in extracting information on subsaturation EOS through collective flows and particle production (at intermediate beam energy) Li (2000); Baran et al. (2005); Li et al. (2008); Chen et al. (2005); Tsang et al. (2009), but also turn out to be the unique tool in terrestrial labs for the exploration of the suprasaturation density behaviors of the nuclear matter EOS (at intermediate to high beam energy) Li (2002a); Xiao et al. (2009); Feng and Jin (2010); Russotto et al. (2011); Xie et al. (2013); Cozma et al. (2013); Russotto et al. (2016); Cozma (2018). In this sense, it is important to improve and upgrade the one-body transport models so that they can help to provide more reliable information on the EOS of asymmetric nuclear matter and the in-medium nuclear effective interaction. To this end, the community of heavy-ion collisions at intermediate energies have actually made great efforts in recent years to improve the transport models for heavy-ion collisions by the program of transport model code comparisons Xu et al. (2016); Zhang et al. (2018).

One of the basic input in one-body transport models is the single nucleon potential (nuclear mean-field potential) under non-equilibrium conditions, which is generally dependent on nucleon momentum. The momentum dependence of single nucleon potential, which has various origins like the exchange term of finite-range nuclear interaction and nuclear short-range correlations, is evident from the observed momentum or energy dependence of nucleon optical model potential Hama et al. (1990); Cooper et al. (1993). This has driven the construction and development of a lot of momentum-dependent mean-field potentials, and they have been extensively employed to study both nuclear matter and heavy-ion collisions Gale et al. (1987); Prakash et al. (1988); Welke et al. (1988); Greco et al. (1999); Danielewicz (2000); Bombaci (2001); Persram and Gale (2002); Das et al. (2003); Chen et al. (2005); Xu and Li (2010); Xu and Ko (2010); Chen et al. (2014); Xu et al. (2015). Most of momentum-dependent mean-field potentials so far applied in one-body transport models are parameterized phenomenologically and they cannot be directly used in nuclear structure calculations. It is thus interesting and constructive to directly employ the momentum-dependent nuclear effective interactions, which are widely used in mean-field calculations of finite nuclei, to one-body transport models. By doing this, experimental observables from finite nuclei and heavy-ion collisions can provide crosschecks for the single nucleon potentials, and thus enhance the understanding on in-medium nuclear effective interactions and the associated nuclear matter EOS.

The Skyrme interaction Skyrme (1956, 1959) has been used very successfully in describing the ground and lowly excited state properties of finite nuclei in mean-field calculations Bender et al. (2003); Stone and Reinhard (2007), as well as heavy-ion collisions at low energies in time-dependent Hartree-Fock (TDHF) calculations Maruhn et al. (2014); Nakatsukasa et al. (2016). However, its incorrect momentum dependence at high kinetic energies (above about $300~{}\rm MeV/nucleon$), which fails to reproduce the empirical results on the nucleon optical potential obtained by Hama et al. Hama et al. (1990); Cooper et al. (1993), hinders the application of the conventional Skyrme interaction Vautherin and Brink (1972); Vautherin (1973); Chabanat et al. (1997, 1998) in transport model calculations for heavy-ion collisions at intermediate and high energies. Recent development of quasi-local energy density functional and Skyrme pseudopotential (the terminology represents effective interactions with quasi-local operators depending on spatial derivatives) makes it possible to incorporate the Skyrme interaction into the transport model calculations for intermediate to high energy heavy-ion collisions. This is achieved by introducing additional higher order derivative terms (or higher order momentum dependence) in the conventional Skyrme effective interaction. In Ref. Wang et al. (2018), based on the next-to-next-to-next-to leading order (N$3$LO) Skyrme pseudopotential, the Hamiltonian density and single nucleon potential under general non-equilibrium condition have been given, and three N$3$LO Skyrme pseudopotentials with single particle potential consistent with empirical optical potential up to $1~{}\rm GeV$ kinetic energy have been constructed. This provides the possibility to study finite nuclei and heavy-ion collisions at incident energy up to about $1~{}\rm GeV/nucleon$ (where nuclear matter with about three times of saturation density can be formed Li (2002b)) on the same footing by using the same nuclear effective interaction.

Another important aspect of developing the transport model for heavy-ion collisions is to improve the numerical stability as well as to guarantee the energy and momentum conservation during the dynamic process. The lattice Hamiltonian method Lenk and Pandharipande (1989) provides a good recipe for this motivation. The lattice Hamiltonian method conserves the total energy exactly and the total momentum to a high degree of accuracy, and has been successfully employed into the study of heavy-ion collisions Gallego et al. (1991); Xu et al. (1990, 1991); Xu (1991, 1992a, 1992b); Xu et al. (1993); Greco et al. (1999); Persram and Gale (2002); Das Gupta et al. (2013); Mallik et al. (2014, 2015). As the first step of building a transport model with the lattice Hamiltonian method by incorporating the N$3$LO Skyrme pseudopotential, in the present work we develop a one-body transport model without considering nucleon-nucleon collisions based on the lattice Hamiltonian Vlasov (LHV) method Lenk and Pandharipande (1989) with nuclear mean-field potential based on the N$3$LO Skyrme pseudopotential. In order to ensure the reliability of the LHV method, the initial phase space distribution is obtained self-consistently through varying the total energy. The ground state properties of the present LHV transport model are examined. Such a LHV transport model without nucleon-nucleon collision term can be applied to the study of dynamical evolution of quantum systems where the individual collisions are either inhibited by the Pauli exclusive principle or negligible due to the diluteness of the system. In that case, collective motions of finite nuclei are ideal sites to test the validity of the present LHV method. To this end, isoscalar giant monopole and isovector giant dipole modes of \isotope[208]Pb are studied based on the present LHV method. The obtained results are compared with that from random-phase approximation (RPA) and experimental data.

The present article is organized as follows. In Sec. II, we introduce the mean-field potential used in the LHV method, namely, N$3$LO Skyrme pseudopotential, and its Hamiltonian density under general nonequilibrium conditions. In Sec. III, the details of the LHV method are given, and the initialization of nuclear ground state and the treatment of collective motions used in the present work are also introduced. The ground state properties, isoscalar giant monopole, and isovector giant dipole modes of \isotope[208]Pb based on the LHV method with the Skyrme pseudopotentials are presented in Sec. IV. Finally, we summarize the present work and make a brief outlook in Sec. V.

II Mean-field potential

II.1 N3LO Skyrme pseudopotential

The quasilocal nuclear energy density functional (EDF) based on the density-matrix expansion provides an efficient way to investigate the universal EDF of a nuclear system. In previous literature Carlsson et al. (2008); Raimondi et al. (2011), the Skyrme interaction has been recognized as the corresponding pseudopotential of quasilocal nuclear EDF through Hartree-Fock (HF) approximation, and a mapping has been established from the N$3$LO local EDF Carlsson et al. (2008) to N$3$LO Skyrme pseudopotential. Such a mapping is worthwhile because it provides an order-by-order way to examine the validity of each term in the quasilocal effective interaction, since the precise structure of nuclear EDF can be derived from low-energy quantum chromodynamics with chiral perturbation theory Puglia et al. (2003); Kaiser et al. (2003); Finelli et al. (2006). The N$3$LO Skyrme pseudopotential Carlsson et al. (2008); Raimondi et al. (2011) is a generalization of the standard Skyrme interaction Vautherin and Brink (1972); Vautherin (1973); Chabanat et al. (1997, 1998) by adding terms that depend on derivative operators (momentum operator) up to sixth order, which corresponds to the expansion of the momentum space matrix elements of a generic interaction in powers of the relative momenta up to the sixth order. In this sense, the standard Skyrme interaction can be regarded as a N$1$LO Skyrme pseudopotential. Such kind of generalization of the Skyrme interaction has been employed to describe EOS of nuclear matter Davesne et al. (2013, 2014, 2015a, 2015b, 2016), as well as the properties of finite nuclei Carlsson et al. (2010); Becker et al. (2017).

The full Skyrme pseudopotential generally contains central, spin-orbit, and tensor components. Since in the present LHV method, and in most of the one-body transport models, only spin-averaged quantities are taken into consideration, we thus ignore the spin-orbit and tensor components, which are irrelevant to spin-averaged quantities. The corresponding extended Skyrme interaction used in the present LHV method is written as

[TABLE]

with an overall factor $\hat{\delta}(\vec{r}_{1}-\vec{r}_{2})$ omitted for the sake of clarity. The central term is expressed as

[TABLE]

where $\hat{\vec{k}}^{\prime}$ and $\hat{\vec{k}}$ are derivative operators acting on left and right, and they take the conventional form as $-(\hat{\vec{\nabla}}_{1}-\hat{\vec{\nabla}}_{2})/2i$ and $(\hat{\vec{\nabla}}_{1}-\hat{\vec{\nabla}}_{2})/2i$, respectively. $\hat{P}_{\sigma}$ represents the spin exchange operator defined as $\hat{P}_{\sigma}=\frac{1}{2}(1+\hat{\sigma}_{1}\hat{\sigma}_{2})$, where $\hat{\sigma}_{1}$ and $\hat{\sigma}_{2}$ are Pauli matrices acting on the first and second state, respectively. The density-dependent term $V^{DD}_{\text{N1LO}}$, which is introduced to mimic phenomenologically the effects of many-body interactions, is taken to be the same form as in the standard Skyrme interaction, i.e.,

[TABLE]

In the above expressions, $t^{[n]}_{i}$, $x^{[n]}_{i}$ ($n=2,4,6$ and $i=1,2$), $t_{0}$, $t_{3}$, $x_{0}$, $x_{3}$ and $\alpha$ are Skyrme parameters. In particular, the parameters $t_{0}$ and $t^{[n]}_{i}$ are measures of the mean central potential, with an spin-exchange character specified by $x_{0}$ and $x^{[n]}_{i}$. The parameters $t_{3}$ and $x_{3}$ are their analogs of the density-dependent potential, with $\alpha$ characterizing its density dependence. The parameters $t^{[n]}_{i}$ and $x^{[n]}_{i}$ are related to the derivative terms which determine the momentum-dependent parts of the Hamiltonian density and the single-nucleon potentials. With the introducing of the additional derivative terms in Eq. (II.1), the N$3$LO Skyrme pseudopotential is able to describe the empirical single nucleon potential up to a kinetic energy of $1~{}\rm GeV$Wang et al. (2018).

II.2 Hamiltonian density with N$3$LO Skyrme pseudopotential

During the heavy-ion collision process, the dinuclear system is generally far from equilibrium. In one-body transport models, such nonequilibrium conditions are described by the phase space distribution function (Wigner function) $f_{\tau}(\vec{r},\vec{p})$, with $\tau$ $=$ $1$ (or n) for neutrons and $-1$ (or p) for protons. In the LHV method, we thus need to express the Hamiltonian density ${\cal H}(\vec{r})$ of the collision system in terms of $f_{\tau}(\vec{r},\vec{p})$. The Hamiltonian density in HF approximation is obtained by calculating the expectation value of the total energy of the collision system. For N$3$LO Skyrme pseudopotential expressed in Eq. (1), it takes the following form (detailed derivation can be found in Ref. Wang et al. (2018)):

[TABLE]

with ${\cal H}^{\rm kin}(\vec{r})$, ${\cal H}^{\rm loc}(\vec{r})$, ${\cal H}^{\rm MD}(\vec{r})$, ${\cal H}^{\rm grad}(\vec{r})$, and ${\cal H}^{\rm DD}(\vec{r})$ being the kinetic, local, momentum-dependent, gradient, and density-dependent terms, respectively. The kinetic term

[TABLE]

and the local term

[TABLE]

are the same as those from the conventional Skyrme interaction. The $\rho_{\tau}(\vec{r})$ in the local term is the nucleon density, which is related to $f_{\tau}(\vec{r},\vec{p})$ through $\rho_{\tau}(\vec{r})$ $=$ $\int f_{\tau}(\vec{r},\vec{p})d\vec{p}$, while $\rho(\vec{r})$ represents the total nucleon density with $\rho(\vec{r})$ $=$ $\rho_{n}(\vec{r})$ $+$ $\rho_{p}(\vec{r})$.

The momentum-dependent term and gradient term contain the contributions from additional derivative terms in Eq.(II.1). The momentum-dependent term can be expressed as

[TABLE]

with $f(\vec{r},\vec{p})$ $=$ $f_{n}(\vec{r},\vec{p})$ $+$ $f_{p}(\vec{r},\vec{p})$. The ${\cal K}_{\rm s}(\vec{p},\vec{p}^{\prime})$ and ${\cal K}_{\rm v}(\vec{p},\vec{p}^{\prime})$ in Eq. (7) represent the isoscalar and isovector momentum-dependent kernel of mean-field potential, respectively. For N$3$LO Skyrme pseudopotential used in the present work, ${\cal K}_{\rm s}(\vec{p},\vec{p}^{\prime})$ and ${\cal K}_{\rm v}(\vec{p},\vec{p}^{\prime})$ take the following forms,

[TABLE]

In the present work, for the sake of simplicity of numerical derivatives, the gradient term is truncated at the second order, and only the isospin symmetric part is taken into account. In other words, we only keep the second order derivative of the total baryon density $\rho(\vec{r})$, i.e.,

[TABLE]

The complete gradient term in the Hamiltonian density for the N$3$LO Skyrme potential can be found in Ref. Wang et al. (2018).

In the above expressions, for convenience, we have recombined the Skyrme parameters related to the derivative terms in Eq. (II.1), namely, $t_{1}^{[n]}$, $t_{2}^{[n]}$, $x_{1}^{[n]}$, and $x_{2}^{[n]}$, into the parameters $C^{[n]}$ and $D^{[n]}$,

[TABLE]

which are related to the momentum-dependent terms, and $E^{[2]}$, which is related to the gradient terms with

[TABLE]

The density-dependent term is expressed as

[TABLE]

Based on the above expressions, one can see that the Hamiltonian density ${\cal H}(\vec{r})$ is explicitly dependent on $f_{\tau}(\vec{r},\vec{p})$, as well as the $\rho_{\tau}(\vec{r})$ and their derivatives.

The calculations in the present work are mainly based on three N$3$LO Skyrme pseudopotentials, SP$6$s, SP$6$m and SP$6$h Wang et al. (2018), which can describe the empirical single nucleon potential up to $1~{}\rm GeV$ in kinetic energy. The main difference of these three interactions is their suprasaturation behavior of the isospin dependent part of EOS, namely, the symmetry energy. Since the nucleon momenta in both the ground state properties and low energy collective excitation are not much larger than the Fermi momentum of saturated nuclear matter, conventional Skyrme interactions are still able to give the empirical single nucleon potential. In order to compare with RPA calculations, one conventional Skyrme interaction MSL$1$ Zhang and Chen (2013) is also adopted in the present work. We list the Skyrme parameters related to nuclear matter properties for the above four Skyrme interactions in Table 1, while more discussions about the gradient parameter $E^{[2]}$, which is irrelevant to nuclear matter properties, will be given in Sec. III.2. For further reference, the characteristic parameters of nuclear matter for these interactions are shown in Table 2. The definition of these quantities can be found, e.g., in Ref. Chen et al. (2009).

III Lattice Hamiltonian Vlasov method

III.1 Lattice Hamiltonian method for Vlasov equation

Quantum theory with phase-space distributions, with proper generalization or approximation, is profoundly suitable to formulate and solve many-particle dynamics Carruthers and Zachariasen (1983). It is demonstrated that in the limit involving $\hbar$ $\rightarrow$ [math], quantum theory with one-body phase-space distributions is reduced to the Vlasov equation Bertsch and Das Gupta (1988),

[TABLE]

where $f$ is the one-body phase-space distribution, or Wigner function defined as the Wigner transform of one-body density matrix $\rho(\vec{r}+\vec{s}/2,\vec{r}-\vec{s}/2)$, i.e.,

[TABLE]

In nuclear physics, Eq. (15) with an additional nucleon-nucleon collision term on the right hand side, which takes into account Fermi statistics, i.e.,

[TABLE]

is commonly referred to as the BUU equation Bertsch and Das Gupta (1988).

The Vlasov equation or BUU equation is normally solved by interpreting $f_{\tau}(\vec{r},\vec{p},t)$ as the semi-classical phase space distribution function. If we treat each volume element as nuclear matter, the (quasi)nucleons inside the volume are in momentum eigenstates obeying the Pauli principle. Under such a condition, the obtained $f_{\tau}(\vec{r},\vec{p},t)$ through Wigner transform of the density matrix turns out to be the occupation probability of the momentum eigenstates, and thus can be regarded as the classical phase space distribution function.

The test particle method Wong (1982), where the semiclassical $f_{\tau}(\vec{r},\vec{p},t)$ is mimicked by a large number of test particles, i.e.,

[TABLE]

has been introduced to solve the Vlasov equation numerically. In the conventional test particle method, the evolutions of coordinate $\vec{r}_{i}(t)$ and momentum $\vec{p}_{i}(t)$ of the $i$th test nucleon are governed by the mean-field potential, or the single particle potential, which can be obtained either by varying the Hamiltonian density or by direct parametrization. In order to obtain a smooth mean-field potential, the density of a certain cell is averaged by neighboring cells. This kind of simple smoothing technique violates the equation of motion a little bit, and fails to conserve either the total energy or the total momentum Bertsch and Das Gupta (1988).

The lattice Hamiltonian Vlasov method developed by Lenk and Pandharipande Lenk and Pandharipande (1989) overcomes this disadvantage of conventional test particle method. It conserves the total energy exactly and the total momentum to a high degree of accuracy. In the LHV method, instead of mean-field potential, the equation of motion of test nucleons is governed directly by the total Hamiltonian of the system, which is approximated by the lattice Hamiltonian, i.e.,

[TABLE]

where $\vec{r}_{\alpha}$ represents the coordinate of a certain lattice site $\alpha$ and $l_{x}$, $l_{y}$, and $l_{z}$ are lattice spacing. The above lattice Hamiltonian can be expressed in terms of the positions and momenta of test nucleons, if we write the semi-classical LHV phase space distribution at lattice site $\alpha$ as

[TABLE]

where the factor $\frac{1}{2}$ is due to spin degeneracy. $N_{\rm E}$ is the number of ensembles (or number of test particles in some literature) introduced in the calculation, and the sum runs over all test nucleons with isospin $\tau$ that contribute to the lattice site $\alpha$. This equivalently gives each test particle a form factor $S$ compared with Eq. (18). By doing this, the movement of a test particle leads to a continuous variation of the local nucleon density of the nearby lattice sites, which is useful to smooth the nucleon distribution functions in phase space. It should be noted that the form factor $S$ actually modifies the relation between test particles and the Wigner function $f$. At this point, we would like to point out that in principle, a similar form factor in momentum space can also be introduced in Eq. (20), which is expected to improve the calculations with momentum-dependent mean-field potentials. In the present work, we only adopt the form factor in coordinate space, and in the future it would be definitely interesting to perform a systematic investigation on the effects of a form factor in momentum space in heavy-ion transport model calculations. The local nucleon density at lattice sites, or LHV density, is then given by integrating $\tilde{f}_{\tau}(\vec{r}_{\alpha},\vec{p},t)$ with respect to momentum, i.e.,

[TABLE]

Note that here we distinguish the realistic phase space distribution function $f(\vec{r},\vec{p})$ and local density $\rho(\vec{r})$ from the LHV phase space distribution function and density expressed in Eqs. (20) and (21), respectively, and will explore their distinctions in Sec. III.2. Substituting Eq. (20) into Eq. (4), the lattice Hamiltonian $H_{L}$ is then expressed in terms of the coordinates and momenta of test nucleons, and subsequently they can be treated as the canonical variables of the lattice Hamiltonian. The equation of motion for the $i$th test nucleon is then governed by the Hamilton equation of total lattice Hamiltonian of all ensembles $N_{\rm E}H_{L}$, i.e.,

[TABLE]

In the above two equations, $A$ is the nucleon number of the system, while the subscripts $\alpha$ for various quantities denote their values at lattice site $\alpha$. The sums run over all lattice sites inside $V_{i}$, where the form factor of the $i$th test nucleon covers. A tilde above the Hamiltonian density, e.g., $\tilde{\cal H}^{\rm loc}(\vec{r}_{\alpha})$, denotes that in its expression in Sec. II.2, the realistic phase space distribution function and local density are replaced by the LHV phase space distribution function and density. The Coulomb interaction contributes to the Hamiltonian density through

[TABLE]

among which the minus term represents the contribution from the exchange term of Coulomb energy. Further testing shows that the Coulomb energy $\tilde{\cal H}^{\rm Cou}(\vec{r}_{\alpha})$ defined in the above equation has already converged at lattice spacing of $l_{x}$ $=$ $l_{y}$ $=$ $l_{z}$ $=$ $0.5~{}\rm fm$ used in the present work. The partial derivative of $\tilde{\rho}_{\tau,\alpha}$ in Eq. (23) can be calculated in terms of the spatial derivative of $S$, i.e.,

[TABLE]

The momentum-dependent parts of the equation of motion of test particles in Eqs. (22) and (23) are obtained by substituting the momentum-dependent part of the Hamiltonian density, i.e., Eq. (7) into Eqs. (22) and (23) after replacing $f_{\tau}(\vec{r},\vec{p})$ in Eq. (7) with the semi-classical LHV phase space distribution expressed in Eq. (20). The integrals in Eq. (7) turn out to be summations of test particles,

[TABLE]

Based on Eqs. (22) - (27), we can calculate the time evolution $\vec{r}_{i}(t)$ and $\vec{p}_{i}(t)$ of test nucleons, and then calculate physical observables based on Eq. (20).

In the present work, the factor $S(\vec{r}_{i}-\vec{r})$ is chosen to be triangle form,

[TABLE]

where $\theta$ is the Heaviside function, and $n$ is an integer which determines the range of $S$. Calculations based on lattices generally break Galilean invariance, and thus violate momentum conservation. In the present work we choose $n$ $=$ $4$, which is large enough to conserve the total momentum to a high degree of accuracy Lenk and Pandharipande (1989).

Generally speaking, the choice of $S(\vec{r}_{i}-\vec{r})$ is somewhat arbitrary. Besides the triangle form used in the present work, alternative forms used in previous literature are trapezoid Danielewicz (2000), double parabolic Persram and Gale (2002), and Gaussian Urban (2012) form. However, in order to ensure particle number conservation, $S(\vec{r}_{i}-\vec{r})$ should satisfy the following equation:

[TABLE]

It should be mentioned that in the conventional test particle method, the Hamiltonian equations of motion for the test particles are derived from the single particle Hamiltonian, which makes it difficult to exactly conserve energy of the system in the dynamic process Bertsch and Das Gupta (1988); Lenk and Pandharipande (1989). On the other hand, in the LHV method, the Hamiltonian equations of motion for the test particles, namely, Eqs. (22) and (23), are derived from the total Hamiltonian of the system of test particles, which guarantees the exact energy conservation of the system in the dynamic process Lenk and Pandharipande (1989). In addition, we would like to mention that the present LHV method is implemented based on GPU parallel computing Ruetsch and Fatica (2013), which increases the computational efficiency and makes it possible to obtain more reliable results with the use of much more ensembles.

III.2 Initialization of nuclear ground state

In the present framework, the Vlasov ground state of nuclei at zero temperature is obtained by varying the Hamiltonian with respect to the nuclear radial density. Such kind of initialization in one-body transport model is sometimes referred to as Thomas-Fermi (TF) initialization Lenk and Pandharipande (1989); Danielewicz (2000); Gaitanos et al. (2010); Lin and Danielewicz (2019). Within the one-body transport models, at zero temperature, for a nucleus in ground state, its Wigner function satisfies

[TABLE]

where $p^{F}_{\tau}(\vec{r})$ is local Fermi momentum and fulfills

[TABLE]

For simplicity in the following we assume the nucleus is spherical. The total energy of a ground state nucleus at zero temperature can be then treated as a functional of radial density $\rho(r)$ and its spatial gradients $\nabla_{r}^{n}\rho(r)$, i.e.,

[TABLE]

After varying the total energy with respect to $\rho_{\tau}(r)$ and its spatial gradients, and considering Eqs. (4) (31), we obtain the neutron/proton radial density in a ground state nucleus (note that the contribution from the Coulomb interaction in Eq. (24) for protons should also be included in the Hamiltonian density),

[TABLE]

where $\mu_{\tau}$ is the chemical potential of proton or neutron inside the nucleus and U_{\tau}\big{\{}p_{\tau}^{\rm F}\big{[}\rho_{\tau}(r)\big{]},r\big{\}} is the single nucleon potential of the nucleon at the Fermi surface. The single nucleon potential is defined as the variation of Hamiltonian density with respect to the phase space distribution function (or local density for the zero temperature case) and density gradients. For the N$3$LO Skyrme pseudopotential, the detailed expression can be found in Ref. Wang et al. (2018). Equation (34) has a very intuitive physical significance: in a classical point of view, it means within a ground state nucleus, the nucleons possessing different Fermi momentum in different radial position have the same chemical potential. The local density $\rho(\vec{r})$ for ground state spherical nucleus is obtained by solving Eq. (34) subjecting to the boundary condition

[TABLE]

where $r_{\rm B}$ is the boundary of the nucleus satisfying $\rho(r_{\rm B})$ $=$ [math], and needs to be determined when solving Eq. (34). Since the way we deal with the ground state is semiclassical, the gradient parameter $E^{[2]}$ is readjusted for each interaction in order that the solution of Eq. (34) reproduces roughly the experiment binding energy and charge radius of \isotope[208]Pb. For MSL$1$ interaction, the results in Sec. IV are based on the readjusted gradient parameters.

In the present LHV method, the initial coordinates of test nucleons are generated according to the solution of Eq. (34), while their initial momenta follow zero-temperature Fermi distribution with the Fermi momentum in Eq. (32) determined by local density.

Careful readers might realize that the obtained ground state LHV density $\tilde{\rho}(\vec{r})$ is smeared due to the form factor introduced in Eq. (20), and thus slightly different from the solution of Eq. (34). Contrary to the Gaussian wave packet in quantum molecular dynamics Aichelin (1991), we do not attach any physical meaning to the form factor $S(\vec{r}-\vec{r}^{\prime})$ and the smoothed LHV density $\tilde{\rho}$, and regard them as numerical techniques introduced in the test-particle approach so that we can obtain well-defined densities and mean fields.

In order to compensate for effects caused by the smearing of the local density due to the form factor, an additional gradient term should appear in the local density based on the following argument. The LHV density $\tilde{\rho}$ can be regarded as the convolution of the realistic local density,

[TABLE]

To express $\rho$ in terms of $\tilde{\rho}$, formally we have

[TABLE]

where $S^{-1}(\vec{r}-\vec{r}^{\prime})$ is the inverse of $S(\vec{r}-\vec{r}^{\prime})$ which satisfies

[TABLE]

The parameter $c$, defined as

[TABLE]

is a constant only depending on the form of $S$.

In the LHV method, the direct correction on $\tilde{\rho}(\vec{r})$ is not feasible since numerically the density in Eq. (37) is not always positively defined. In practice, to compensate the smearing of density due to the form factor, we introduce an additional gradient term $\tilde{E}^{[2]}\nabla^{2}\tilde{\rho}(\vec{r})$ into the Hamiltonian. To demonstrate this, one need only substitute Eq. (37) into Eq. (4), and after several necessary approximations, an additional term proportional to $c\tilde{\rho}(\vec{r})\nabla^{2}\tilde{\rho}(\vec{r})$ will show up in the Hamiltonian. Though $\tilde{\rho}$ is not a constant, in the present work, for simplicity the additional gradient term is recast effectively into $\tilde{E}^{[2]}\nabla^{2}\tilde{\rho}$, with $\tilde{E}^{[2]}$ being a constant. This is equivalent to replace the gradient term coefficient $E^{[2]}$ by $E^{[2]}$ $+$ $\tilde{E}^{[2]}$. Since in principle the rms radius in exact ground state does not change with time, $\tilde{E}^{[2]}$ here is adjusted roughly to obtain the ground state rms radius evolution with the smallest oscillation. It should be mentioned that this correction on gradient term only improves the stability of ground state evolution (rms radius and radial density profile) slightly, and does not cause much difference on the result of collective motions in Sec. IV. Needless to say, in ideal cases with $N_{\rm E}$ approaching infinity and lattice spacing approaching zero, LHV local density will approach the realistic local density, and $\tilde{E}^{[2]}$ will become zero.

In Table 3, for the interactions used in the present work, namely SP$6$s, SP$6$m, SP$6$h, and MSL$1$, we list their parameters $E^{[2]}$ and $\tilde{E}^{[2]}$, as well as, for \isotope[208]Pb, the binding energy and proton rms radius based on TF initialization, i.e., the solution of Eq. (34), and in the LHV calculations. We note from the table that after the smearing of $S$ in LHV method, the total energy decreases and the rms radius of proton increases slightly. However, this difference affects the stability of the ground state very little, as we will see in Sec. IV.1.

III.3 Nuclear giant resonance within the Vlasov equation

We consider a perturbative excitation of the Hamiltonian at the initial time, i.e.,

[TABLE]

where $\hat{Q}$ is the excitation operator we are interested in and $\lambda$ is supposed to be small. Within the linear response theory Fetter and Walecka (1971), the response of the excitation operator $\hat{Q}$ as a function of time is given by

[TABLE]

where $|0\rangle$ is the ground state for unperturbed system, $|f\rangle$ is the energy eigenstate of the excited system, $E_{0}$ and $E_{f}$ are the eigen-energy of the system before and after excitation, respectively.

We define the strength function $S(E)$ as usual through

[TABLE]

The $S(E)$ can be expressed as the Fourier integral of $\Delta\langle\hat{Q}\rangle(t)$ by taking advantage of Eq. (41), i.e.,

[TABLE]

By evaluating the time evolution of the response of $\hat{Q}$ within the LHV method, we can obtain the strength function.

We assume $\hat{Q}$ is a one-body operator, which means it can be expressed as the sum of operators $\hat{q}$ acting on each nucleon, i.e.,

[TABLE]

The expectation of $\hat{Q}$ then can be calculated as follows,

[TABLE]

where we have added two identity operators. Considering the definition of the one-body density matrix,

[TABLE]

and combining with Eq. (44), we rewrite Eq. (45) as

[TABLE]

For convenience in the following we change the variables of the integral, $\vec{r}_{1}$ $=$ $\vec{r}+\frac{\vec{s}}{2}$ and $\vec{r}_{1}^{\prime}$ $=$ $\vec{r}-\frac{\vec{s}}{2}$. Since $f(\vec{r},\vec{p})$ is the Wigner transform of density matrix, in coordinate space the density matrix can be expressed as

[TABLE]

We define the Wigner transform of $\hat{q}$ in coordinate space,

[TABLE]

where

[TABLE]

is the matrix element of $\hat{q}$ in coordinate space. Substituting Eq. (47) and the inverse transform of Eq. (48) into Eq. (46), we reduce the expectation of $\hat{Q}$ into the following form,

[TABLE]

In the LHV method, the $f(\vec{r},\vec{p})$ is replaced by $\tilde{f}(\vec{r},\vec{p})$ expressed in Eq. (20). It can be demonstrated that the effect of the external excitation $\lambda\hat{Q}\delta(t-t_{0})$ is to change the positions and momenta of the test nucleons as follows Urban (2012):

[TABLE]

For the isoscalar monopole and isovector dipole excitation, the specific form of Eq. (51) will be given in Secs. IV.2 and IV.3, respectively.

IV Result and discussion

In the following, we are going to examine the ability of the present LHV method in dealing with the (near-)equilibrium nuclear dynamics. Specifically, we will study the ground state evolution, isoscalar monopole and isovector dipole mode of \isotope[208]Pb. Under such (near-)equilibrium states, most of the nucleon-nucleon collisions are blocked according to the Pauli exclusive principle, thus in principle the BUU equation with the absence of the collision term, or the Vlasov equation, is still applicable.

IV.1 Ground state evolution of finite nuclei

As mentioned in Sec. III.2, the initial phase space information of a ground state nucleus is obtained self-consistently by varying the total energy with respect to nucleon radial density $\rho(r)$. The initial coordinates of test nucleons are generated according to the solution of Eq. (34), and the LHV density $\tilde{\rho}$ is obtain via Eq. (21). Shown in Fig. 1 is the time evolution of radial profile of LHV density for a single \isotope[208]Pb in ground state up to $200~{}{\rm fm}/c$, obtained from the present LHV method with $N_{\rm E}$ $=$ $10000$ and time step $\Delta t$ $=$ $0.4~{}{\rm fm}/c$ by using the N$3$LO Skyrme pseudopotential SP$6$m. We notice from Fig. 1 that the LHV density approximates reasonably well the realistic ground state of the Vlasov equation, or the solution of Eq. (34), since the profile of the LHV density only exhibits very small fluctuations.

To see more clearly the stability of the ground state evolution of \isotope[208]Pb within the present LHV method, the ground state evolution is continued up to $1000~{}{\rm fm}/c$, and we present in Fig. 2 the radial profiles of LHV density with a time interval of $200~{}{\rm fm}/c$. Again, only small fluctuations are observed in Fig. 2, which indicates that the present LHV method is capable of studying long-time nuclear processes, e.g., heavy-ion fusion reactions and nuclear spallation reactions.

Apart from the radial density profile, other properties concerning the stability of the LHV method are also examined. To that end, the time evolution of rms radius, fraction of bound nucleons and binding energy are presented in Fig. 3. The calculations are performed with time step $\Delta t$ $=$ $0.4~{}{\rm fm}/c$, and $N_{\rm E}$ $=$ $5000$ and $10000$, respectively. Free test nucleons are those who do not interact with other test nucleons (which means their form factors $S$ do not overlap), and they are excluded in calculating the fraction of bound nucleon and rms radius. We notice from Fig. 3(a) that both cases with different $N_{\rm E}$ give fairly stable evolution of rms radius. For the $N_{\rm E}$ $=$ $5000$ case, the rms radius starts to decrease after about $800~{}{\rm fm}/c$. This decrease is due to the evaporation of nucleons from the bound nuclei, which is demonstrated in Fig. 3(b). Such evaporation of nucleons is inevitable in transport model simulations due to the limited precision in the numerical realization, whereas it can be suppressed by increasing $N_{\rm E}$, as can be seen in Fig. 3(b), though the result with $E_{\rm E}$ $=$ $5000$ is already satisfactory. Owing to the advantage of lattice Hamiltonian framework, the binding energy of the given nucleus conserves almost exactly in both cases, as shown in Fig. 3(c).

Figure 3 indicates that the present LHV method is able to give fairly stable time evolution. Due to the high efficiency of GPU parallel computing, including more ensembles in the LHV calculation becomes possible. However, as a balance between computational resources and numerical accuracy, unless otherwise specified, the following calculations are performed with time step $\Delta t$ $=$ $0.4~{}{\rm fm}/c$ and $N_{\rm E}$ $=$ $5000$.

IV.2 Isoscalar monopole mode of \isotope[208]Pb

During the past several decades, a lot of studies on the isoscalar giant monopole resonance (ISGMR) of finite nuclei have been performed, since it provides the information about both the symmetric and asymmetric part of the nuclear matter incompressibility Youngblood et al. (1999); Shlomo and Youngblood (1993); Li et al. (2007); Patel et al. (2012, 2013); Gupta et al. (2016), which are fundamental quantities characterizing the EOS of nuclear matter. In the experimental aspect, the isoscalar monopole mode is measured through scattering off nucleus with isoscalar light particles, and recent experiments have been performed with inelastic $\alpha$-particle and deuteron scatterings Li et al. (2007); Monrozeau et al. (2008); Patel et al. (2012, 2013); Vandebrouck et al. (2014); Gupta et al. (2016).

In the one-body transport model point of view, the isoscalar monopole mode is regarded as a compressional breathing of nuclear fluid. Such mode can be generated in the LHV framework through the following procedures. For the isoscalar monopole mode, we have

[TABLE]

and thus according to Eq. (48), we obtain

[TABLE]

Note that the square root of the expectation value of $\hat{Q}_{\rm ISM}$ is the rms radius of the given nucleus. According to Eq. (51), to obtain the isoscalar monopole mode, the initial phase space information of test nucleons are changed with respect to that of the ground state by

[TABLE]

The spatial coordinators of test nucleons remain unchanged since the $q_{\rm ISM}$ is independent of momentum.

In Fig. 4 we show the time evolution of $\Delta\langle\hat{Q}\rangle_{\rm ISM}$ with three N$3$LO Skyrme pseudopotentials, namely, SP$6$s, SP$6$m, and SP$6$h, as well as one conventional Skyrme interaction MSL$1$. The perturbation parameter $\lambda$ is chosen to be $100~{}{\rm MeVfm^{-1}}/c$ in the calculation. We notice from the figure that the time evolution of $\Delta\langle\hat{Q}\rangle_{\rm ISM}$, or equivalently the rms radius, exhibits a very regular oscillation, and the rapid increase which is commonly observed in most BUU simulations using the conventional test particle method does not show up here. Besides that, only a slight damping is observed in the oscillation, which is anticipated since the only damping mechanism in the Vlasov framework is Landau damping. Landau damping is caused by one-body dissipation which is governed by a coupling of single-particle and collective motion. It should be mentioned that in the RPA framework, the damping also only comes from one-body dissipation, since the coupling to more complex states like two-particle two-hole ($2p$-$2h$) states is missing in RPA Bertsch et al. (1983). In the semiclassical framework, effects analogous to the $2p$-$2h$ excitation can be included via a nucleon-nucleon collision term Burgio and Di Toro (1988), and the width of the strength indeed increases due to the inclusion of the collision term Gaitanos et al. (2010), but this is beyond the scope of the present work, and will be pursued in the future.

The strength function is obtained based on the time evolution of the response presented in Fig. 4 and Eq. (44). When calculating the strength function, the response $\Delta\langle\hat{Q}\rangle(t)$ is multiplied by an exponential attenuation $e^{-\gamma t/2\hbar}$, which is a common practice in many giant resonance calculations within the BUU model Urban (2012); Kong et al. (2017). The reason for introducing such an attenuation is to avoid oscillations in the Fourier transforms of Eq. (43), which is caused by the finite period of the evolution. In this work the attenuation parameter $\gamma$ is set to be $2~{}\rm MeV$, as well as in Sec. IV.3 for the isovector dipole mode. However, this exponential attenuation will not affect the peak energy of the strength, on which we mainly focus.

The obtained strength functions of the isoscalar monopole mode with SP$6$s, SP$6$m, SP$6$h, and MSL$1$ are presented in Fig. 5. The grey band represents the peak energy of $13.91\pm 0.11~{}\rm MeV$ from inelastic $\alpha$-scattering off \isotope[208]Pb performed at TAMU Youngblood et al. (1999) while the cyan band of $13.7\pm 0.1~{}\rm MeV$ represents that from the Research Center for Nuclear Physics (RCNP) Patel et al. (2013). We notice from the figure that all these interactions give peak energies consistent with that of the experiment, which is a consequence of the proper nuclear incompressibility of the given interactions, as shown in Table 2 by $K_{0}$. In order to compare the results from the LHV method with that from different approaches, we calculate the strength function of MSL$1$ based on RPA. The Skyrme-RPA code by Colo et al. Colò et al. (2013) is employed. The obtained peak energy is indicated by a green arrow. The peak energies of MSL$1$ with different approaches are generally comparable, and the small discrepancy comes from the difference between semi-classical and quantum nature.

IV.3 isovector dipole mode of \isotope[208]Pb

The isovector giant dipole resonance (IVGDR) of finite nuclei is the oldest known nuclear collective motion. Systematic experimental studies on IVGDR based on photo-nuclear reactions can be traced back to more than forty years ago Berman and Fultz (1975). Recent measurements on isovector dipole response have been performed based on inelastic proton scattering at RCNP for \isotope[48]Ca Birkhan et al. (2017), \isotope[120]Sn Hashimoto et al. (2015), and \isotope[208]Pb Tamii et al. (2011), as well as by using Coulomb excitation in inverse kinematics at GSI for \isotope[68]Ni Rossi et al. (2013). It is interesting to mention that in recent years a low-lying mode called pygmy dipole resonance (PDR) has attracted a lot of attention both experimentally Ryezayeva et al. (2002); LAND-FRS Collaboration et al. (2005); Wieland et al. (2009); Endres et al. (2010) and theoretically Carbone et al. (2010); Urban (2012); Baran et al. (2013, 2014). It is well known from theoretical studies based on various models that the PDR, IVGDR, and electric dipole polarizability $\alpha_{D}$ which is dominated by these isovector dipole modes, provide sensitive probes of the density dependence of the nuclear symmetry energy Yildirim et al. (2005); Trippa et al. (2008); Piekarewicz et al. (2012); Roca-Maza et al. (2013); Zhang and Chen (2015); Roca-Maza et al. (2015); Kong et al. (2017).

Within the LHV method, the external perturbation for the isovector dipole mode can be expressed in the following form:

[TABLE]

which is defined so that the center of mass of the nucleus stays at rest. Similar to the case in the isoscalar monopole mode, the excited system is obtained by changing the initial phase space coordinates of test nucleons according to

[TABLE]

After initializing and exciting the nucleus, we observe in Fig. 6 the damped oscillations of the isovector dipole response in \isotope[208]Pb with SP$6$s, SP$6$m, SP$6$h, and MSL$1$. In the isovector dipole case, $\lambda$ is set to be $25~{}{\rm MeV}/c$. Similar to the isoscalar monopole case, in Fig. 7 the corresponding strength functions of the isovector dipole response are displayed. The vertical cyan line represents the experimental peak energy of $13.4~{}\rm MeV$ obtained from $\isotope[208]{Pb}(p,p^{\prime})$ reaction performed at RCNP Tamii et al. (2011) and the green arrow indicates the peak energy from the RPA calculation in MSL$1$. Since these interactions give a reasonable description for the empirical isospin dependent behaviors in (sub)saturation density, as indicated in Table 2, the peak energy from different interactions is consistent with experimental data. The peak energies of MSL$1$ from the LHV method and the RPA calculation are comparable as well.

In Fig. 6, the obtained responses with SP$6$s, SP$6$m, and MSL$1$ show clearly nonsingle frequency behavior. This is also evident in Fig. 7 from the additional strength or even extra bumps in the high energy region. Further tests show that the bump of the high energy part is related to the magnitude of the gradient terms in Eq. (10), which provides an additional restoring force out of phase with that from the local part. This indicates that the effect of the gradient terms, usually omitted in one-body transport models, should not be overlooked. We notice that such high-energy bumps or peaks have also been observed in previous studies based on one-body transport models with gradient terms Zheng et al. (2016a, b). One possible reason is that the larger the gradient term, the easier for the isovector dipole excitation to invoke other modes, and other modes then react upon the isovector dipole mode and thus cause the additional high-energy strength. In principle, such nonlinear modes can be avoided if the perturbation parameter $\lambda$ is small enough Simenel and Chomaz (2003). However, this is impracticable in one-body transport model simulations due to the limited precision in the numerical realization. Since the high-energy bump or peak is absent in experimental data Tamii et al. (2011), we attribute it to a numerical problem, and further investigation is needed in the future.

V summary and outlook

We have developed a one-body transport model by employing the lattice Hamiltonian method to solve the Vlasov equation with nuclear mean-field obtained from the N$3$LO Skyrme pseudopotential. The ground states of nuclei are obtained with the same interactions through varying the total energy with respect to the nucleon density distribution. Owing to the self-consistent treatment of initial nuclear ground state and the exact energy conservation of the LHV method, the present framework for solving the Vlasov equation exhibits very stable nuclear ground state evolution. As a first application of the new LHV method, we have calculated the isoscalar monopole and isovector dipole modes of finite nuclei. The obtained peak energies of ISGMR and IVGDR in \isotope[208]Pb with the N$3$LO Skyrme pseudopotentials are consistent with the experimental data. In addition, the use of the Skyrme interaction enables us to compare the LHV results with that from conventional nuclear structure method, i.e., the RPA calculation, and the obtained peak energies of ISGMR and IVGDR in \isotope[208]Pb have been shown to be comparable in the LHV method and the RPA calculation.

Our results have demonstrated the capability of the present LHV method in dealing with the ground state of finite nuclei and the near-equilibrium nuclear dynamics. The present work provides a solid foundation not only for long-time Vlasov calculation of low energy heavy-ion reactions, but also for the BUU calculation with a nucleon-nucleon collision term. Based on the future lattice Hamiltonian BUU method, one can use the Skyrme pseudopotentials to simultaneously explore both the structure properties of finite nuclei and heavy-ion collisions at intermediate to high energies. Crosschecks for nuclear effective interactions and thus for the EOS of asymmetric nuclear matter from the nuclear structure and heav-ion collisions thus become possible. Such studies are in progress and will be reported elsewhere.

Acknowledgments

We thank P. Danielewicz and B.-A. Li for helpful discussions, and M. Gao and X. Zhang for the maintenance of the GPU severs. This work was supported in part by the National Natural Science Foundation of China under Grant No. 11625521, the Major State Basic Research Development Program (973 Program) in China under Contract No. 2015CB856904, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, and the Science and Technology Commission of Shanghai Municipality (11DZ2260700).