Nuclear collective dynamics in transport model with the lattice Hamiltonian method

We review the recent progress on studying the nuclear collective dynamics by solving the Boltzmann-Uehling-Uhlenbeck (BUU) equation with the lattice Hamiltonian method treating the collision term by the full-ensemble stochastic collision approach. This lattice BUU (LBUU) method has recently been developed and implemented in a GPU parallel computing technique, and achieves a rather stable nuclear ground-state evolution and high accuracy in evaluating the nucleon-nucleon (NN) collision term. This new LBUU method has been applied to investigate the nuclear isoscalar giant monopole resonances and isovector giant dipole resonances. While the calculations with the LBUU method without the NN collision term (i.e., the lattice Hamiltonian Vlasov method) describe reasonably the excitation energies of nuclear giant resonances, the full LBUU calculations can well reproduce the width of the giant dipole resonance of $^{208}$Pb by including a collisional damping from NN scattering. The observed strong correlation between the width of nuclear giant dipole resonance and the NN elastic cross section suggests that the NN elastic scattering plays an important role in nuclear collective dynamics, and the width of nuclear giant dipole resonance provides a good probe of the in-medium NN elastic cross section.


INTRODUCTION
Transport models deal with the time-evolution of the Wigner function or phase-space distribution function f ( r, p, t) which stems from the Wigner representation of the Schrödinger equation [1,2], and provide a successful semi-classical time-dependent approach to nuclear dynamics, especially to heavy-ion collisions (HICs). One of the main ingredients of transport models is the mean field potential which embodies information on the nuclear equation of state (EOS) or the in-medium effective nuclear interaction. Therefore, transport models serve as an important theoretical tool to explore the EOS of asymmetric nuclear matter from observables in HICs. A lot of information on the nuclear EOS from sub-saturation [3,4,5] to supra-saturation densities up to about 3 ∼ 5 times saturation density [6,7,8,9,10,11,12,13,14,15,16,17,18], has been obtained from transport model analyses of various obseverables, e.g., collective flows and particle production, in intermediate and high energy HICs. The exact information about the nuclear EOS is crucial in describing reaction dynamics of exotic nuclei [19,20], various properties of both finite nuclei (e.g., neutron skin thickness [21,22,23] and drip lines [24,25]) and neutron stars (e.g., masses and cooling mechanisms [26,27,28,29,30]), and astrophysical processes such as supernova explosion scenarios [31,32,33]. In particular, it should be mentioned that the first gravitational wave signal GW170817 [34] of binary neutron star merger has been recently observed and localized by the LIGO and Virgo observatories, and it inaugurates a new era of multimessenger astronomy and gives important constraints on the dense nuclear matter EOS [35,36,37,38,39]. Moreover, very recently, using the X-ray data from NASA's Neutron Star Interior Composition Explorer (NICER), the mass and radius of the millisecond pulsar PSR J0030+0451 have been simultaneously estimated [40,41] and its implications on the dense nuclear matter EOS has been analyzed [42]. In addition, a new record for the maximum mass of neutron stars, namely, a millisecond pulsar J0740+6620 with mass 2.14 +0. 10 −0.09 M ⊙ (68.3% credibility interval), has been recently reported [43], and this heaviest neutron star observed so far can rule out many soft nuclear matter EOS's and especially the supersoft high-density symmetry energy [44].
The time-dependent-Hartree-Fock (TDHF) theory provides a very successful quantum many-body approach at the mean-field level to describe low-energy nuclear reaction dynamics including the nuclear collective dynamics (see, e.g., Refs. [45,46] for recent review). Given that the Vlasov equation, i.e., Boltzmann-Uehling-Uhlenbeck (BUU) equation without the nucleon-nucleon (NN) collision term, corresponds to the semi-classical limit of the TDHF equation, transport models can thus be seen as an efficient semi-classical approach to study nuclear collective dynamics. In particular, the two-particletwo-hole (2p-2h) correlation beyond the mean-field approximation, which dominates the collisional damping of nuclear giant resonances, can be effectively taken into account in transport models via binary collisions. In literature, there have been a lot of works studying nuclear giant resonances based on the pure Vlasov equation [47,48,49], the Vlasov equation with a collision relaxation time [50], and the full transport model with both the mean-field and the NN scatterings [51,52,53]. For example, based on simulations of transport models, the excitation energies of nuclear giant resonances have been used to extract information on the nuclear EOS and neutron-proton effective mass splitting [54], while the width of nuclear giant dipole resonances (GDR) has been proposed as an effective probe of the in-medium NN elastic cross section [55]. The width of nuclear GDR can also serve as a fingerprint of α-particle clustering configurations in nuclei [56].
Although transport models have been extensively used in the study of nuclear giant resonances, the accurate description of the giant resonances within transport models is still a challenge. In transport model calculations, unlike the simulations of HICs at intermediate and high energies, the calculation of nuclear giant resonances, which are the collective excitation states with an excitation energy of about 20 MeV, requires a more proper description of nuclear ground-state and a rather accurate implementation of Pauli blocking. In particular, the Pauli blocking is intimately related to the collisional damping and thus the width of nuclear giant resonances in the transport model calculations. In this sense, the nuclear collective motion provides an ideal topic to examine and improve transport models, since the effects of several deficiencies, such as the inaccurate treatment of Pauli blocking, are more pronounced in nuclear collective dynamics with small amplitude oscillations. Transport models for HICs can be roughly divided into two categories, the BUU equation (see, e.g., Ref. [2]) and the quantum molecular dynamics (QMD) model (see, e.g., Ref. [57]). From the viewpoint of transport models, the essential difference between these two types is that the BUU-type transport models mimic f ( r, p, t) by a large number of ensembles or test particles for each nucleon while the QMD-type by a Gaussian wave packet for each nucleon. Recently, the transport model community has started the code comparison project [58,59,60] to try to understand the source of the discrepancies in various transport model codes, and eventually reduce the uncertainties of transport models. For the issue of Pauli blocking, the QMD-type transport models seem not to be as good as BUU-type models [59], therefore the BUU-type transport models are more suitable for the study of nuclear collective motions, especially for the calculation of the spreading width, in which the accurate treatment of Pauli blocking is essential.
In order to study (near-)equilibrium nuclear dynamics within the framework of transport models, a BUUtype transport model, namely, the lattice BUU (LBUU) method [61,55], has been recently developed, and it can achieve good stability for the ground state evolution [61] and treat Pauli blocking with very high accuracy [55]. The resulting LBUU framework has the following features: 1) a smearing of the local density, which is common practice in transport models to obtain a smooth mean-field, is included selfconsistently in the equations of motion through the lattice Hamiltonian (LH) method; 2) the ground state of a nucleus is obtained through varying the total energy with respect to the nucleon density distribution based on the same Hamiltonian that governs the system evolution; 3) the NN collision term in BUU equation is implemented through a full-ensemble stochastic collision approach. The above techniques as well as the sufficiently large number of ensembles make it possible to solve the BUU equation almost exactly, and thus one can obtain very accurate results for the nuclear collective motions within the BUU equation. We note that the high accuracy of the LBUU method relies on increasing computational resources, therefore the high-performance GPU parallel computing [62] has been employed in the LBUU framework to improve the computing efficiency. This paper is organized as follows. In Sec. 2, we first introduce the LBUU method for solving the BUU equation, including the mean-field, the collision integral and the initialization for the nuclear ground state, and then describe how to deal with the nuclear giant resonances within transport models. In Sec. 3, we present results on the peak energies of the nuclear giant resonances from lattice Hamiltonian Vlasov (LHV) calculations, i.e., the LBUU calculations without the NN collision term, and then compare these with the results from the random-phase approximation (RPA). In Sec. 4, we show the results of the strength function and the width of GDR from the full LBUU calculations, and compare these to the experimental data from 208 Pb( p, p ′ ) reaction carried out at the Research Center for Nuclear Physics in Osaka, Japan (RCNP) [63]. Finally we give a brief summary and outlook in Sec. 5.

MODEL DESCRIPTION
The starting point for calculating the nuclear collective motion is the BUU equation with a momentumdependent mean-field potential U( r, p), i.e., where f (i.e., Wigner function) is the Fourier transform of one-body density matrix ρ( r + s/2, r − s/2), i.e., In the local density approximation, f ( r, p) is reduced to the classical one-body phase-space distribution function. The collision term I c , which takes into account the Pauli principle due to nucleons' Fermi statistics, reads where g = 2 is the spin degeneracy factor, and M 12→34 is the in-medium transition matrix element. Note that we have ignored the isospin index in the above three equations, but it can be restored easily. The BUU equation without the collision term I c is referred to as the Vlasov equation, which is the semi-classical limit of the quantum transport theory with the system described by the one-body phase-space distribution function [1,2], whereas the quantum corrections can be included perturbatively [64,65].
We use the LH method, originally proposed by Lenk and Pandharipande [66] in 1989, to solve the BUU equation. The LH method has been successfully employed in the study of HICs [67,68]. It improves the sample smoothing technique of the usual test particle approach [69], and conserves the total energy almost exactly. In the LH method, the phase-space distribution function f τ ( r, p, t) is mimicked by A × N E test nucleons with a form factor S in the coordinate space to modify the relation between the test nucleons and the Wigner function, i.e., where A is the mass number of the system, and N E is the number of ensembles or number of test particles, usually a very large number, used in the calculation. The sum in the above expression runs over all test nucleons with isospin τ . The form factor S can take a Gaussian form, or a certain form with a finite range that ensures the particle number conservation. By giving each test nucleon a form factor, the movement of a test nucleon leads to a continuous variation of the local nucleon density of the nearby lattice sites, which is useful to smoothen the nucleon distribution functions in phase space. A similar form factor in momentum space [here the δ-function is used in Eq. (4)] could be introduced and this might be helpful to reduce fluctuations if momentum-dependent mean-field potential is employed, and it would be interesting to carry out a systematic investigation of the effects of a form factor in momentum space in the future. The equations of motion of the test nucleons are governed by the total Hamiltonian, and we approximate the latter by the lattice Hamiltonian, i.e., where r α denotes the coordinate of lattice site α, and l x , l y and l z are the lattice spacings. Therefore in the LH method, only the values of the phase-space distribution function at lattice sites f τ ( r α , p, t) need to be calculated.
By solving the BUU equation or Vlasov equation based on the LH method, one obtains the time evolution of f ( r, p, t), or the test nucleons' coordinates r i and momenta p i , and then the time evolution of other physical quantities can be calculated accordingly.

Mean fields
We employ the Skyrme pseudopotential to calculate the lattice Hamiltonian in Eq. (5). The next-tonext-to-next leading order (N3LO) Skyrme pseudopotential [70], which is a mapping of N3LO local energy density functional [71], generalizes the standard Skyrme interaction [72] and can reproduce the empirical nuclear optical potential up to about 1 GeV in kinetic energy [73], for which the standard Skyrme interactions fail to describe. The Hamiltonian density from the N3LO Skyrme pseudopotential contains the kinetic term H kin ( r), the local term H loc ( r), the momentum-dependent term H MD ( r), the density-dependent term H DD ( r) and the gradient term H grad ( r). The kinetic term and the local term are the same as those from the standard Skyrme interaction. The momentum-dependent term is written in the following form, with f ( r, p) = f n ( r, p) + f p ( r, p). The quantities K s ( p, p ′ ) and K v ( p, p ′ ) in Eq. (8) represent the isoscalar and isovector kernels of the momentum-dependent part of the mean-field potential, respectively. K s ( p, p ′ ) and K v ( p, p ′ ) for the N3LO Skyrme pseudopotential are expressed as If we keep only the C [2] and D [2] terms, the N3LO Skyrme pseudopotential reduces to the standard Skyrme effective interaction. For the sake of simplicity of performing numerical derivatives, we truncate at the second order of the spatial gradient of ρ( r), In the second line we have introduced, g [2] = E [2] + 1 2 F [2] , g [2] iso = 1 2 F [2] , and ρ δ = ρ n − ρ p . We neglect the second term in Eq. (11) since it is much smaller than the first term, in other words, we keep only the second order spatial derivative of the total nucleon density ρ( r). The density-dependent term for the N3LO Skyrme pseudopotential takes its form in the standard Skyrme interaction One can see that the Hamiltonian density H( r), expressed as the sum of Eqs. (6)-(8), Eqs. (11) and (12), is explicitly dependent on f τ ( r, p), as well as the densities ρ τ ( r) and their derivatives.
In the above expressions, the parameters C  2 , which are related to derivative terms in the Skyrme two-body potential v Sk ( r 1 , r 2 ) [72], i.e., Specifically, we obtain the coefficient of the gradient term Substituting f ( r, p, t) as expressed in Eq. (4) into Eqs. (6)- (12), and noting that the local nucleon density ρ τ ( r) is given by integrating f τ ( r, p, t) with respect to momentum, we can express the lattice Hamiltonian H L in Eq. (5) in terms of the coordinates and momenta of the test nucleons. Since the coordinate and momentum of the test nucleons r i and p i can be regarded as the canonical variables of the lattice Hamiltonian, their time evolution is then governed by the Hamilton equation for all ensembles, In the above two equations, the quantities with a subscript α refers to their values at lattice site α. The V i under the summation represents the volume, which the form factor of the i-th test nucleon covers, and the sums run over all lattice sites inside V i . The Coulomb interaction contributes to the Hamiltonian density through the term where the second term represents the contribution from the Coulomb exchange energy. Further tests show that the Coulomb energy H Cou ( r α ) converges at a lattice spacing of l x = l y = l z = 0.5 fm used in the present LBUU simulations. The gradient term H grad where have integrated by parts, in order to obtain the second line. The spatial derivative of ρ τ,α in Eq. (20) is related to the spatial derivative of S through Substituting f τ ( r, p) expressed in Eq. (4) into Eq. (8), we obtain the momentum-dependent parts of the equation of motion for the test nucleons, and these are expressed in terms of the sums over the test nucleons by Based on Eqs. (19) - (25), one can evaluate the time evolution of coordinates and momenta r i (t) and p i (t) of the test nucleons, then obtain f ( r, p, t) through Eq. (4), based on which physical observables can be calculated.
The choice of the form factor S( r i − r) should ensure particle number conservation In the present LBUU framework, we use a triangular form where θ is the Heaviside function, and n is an integer which determines the range of S. Generally speaking, calculations on lattices violate momentum conservation since they break Galilean invariance. Early studies have shown that the total momentum can be conserved to a high degree of accuracy if n 4 [66].
It should be mentioned that compared with the conventional test particle method, in which the equations of motion for the test nucleons are derived from single-particle Hamiltonian, the equations of motion for the test nucleons in the LH method, Eqs. (19) and (20), are derived from the total Hamiltonian of the system. For the former way it is difficult to conserve energy exactly [2,66], while the latter can ensure the exact energy conservation , in the dynamic process [66].

Collision integral
In the present LBUU method, instead of the commonly used geometric method, the stochastic collision method [74] is implemented for the NN collision term in the BUU equation. In the stochastic collision approach, the collision probability of two test nucleons can be derived directly from the NN collision term, I c in Eq. (3), as follows. Considering nucleons around lattice site r α from two momentum space volume elements V p 1 = p 1 ± 1 2 ∆ 3 p 1 and V p 2 = p 2 ± 1 2 ∆ 3 p 2 , respectively, one can average over momentum space volume V p i to obtain the distribution function f ( r α , p i ) according to Eq. (4), The number of collisions between nucleons from these two momentum space volumes that happen in a time interval ∆t is The quantities are the changing rates of f ( r α , p 1 ) and f ( r α , p 2 ), respectively, caused by two-body scatterings between the nucleons in V p 1 and V p 2 . These terms can be obtained directly from Eq. (3), i.e., the NN collision term in the BUU equation as where we have substituted by the definition of the cross section with v rel being the relative velocity of the test nucleons in the two momentum space volumes and σ * NN is the scattering cross section in the two-nuclei center-of-mass frame. Here, we obtain the in-medium NN cross section σ * NN through multiplying the free NN cross-section σ free NN by a medium-correction factor. The NN elastic scattering cross section in free space σ free NN is taken from the parameterization in Ref. [75], with cutoff of σ free NN (p lab ≤ 0.1 GeV/c) = σ free NN (p lab = 0.1 GeV/c) for neutron-neutron (nn) or proton-proton (pp) scatterings and σ free NN (p lab ≤ 0.05 GeV/c) = σ free NN (p lab = 0.05 GeV/c) for neutron-proton (np) scatterings, respectively, since the parameterization is shown to be valid for nucleon momentum p lab down to the corresponding cutoff [75]. We note that the p lab cutoff actually corresponds to only a few MeV of incident kinetic energy (i.e., 1.3 MeV for p lab = 0.05 GeV/c and 5.3 MeV for p lab = 0.1 GeV/c) and these very low energy scatterings are not important in the present transport model calculations. Since this parameterization of σ * NN is given in the two-nucleon center-of-mass frame, its value in the two-nuclei center-of-mass frame can be obtained through the Lorentz invariant quantity where the ∆N coll ij denotes the number of physical collisions from the scattering of the i-th and j-th test nucleon. Given that every test nucleon is 1/N E of a physical nucleon, one obtains the collision probability of the i-th and j-th test nucleons as One can reduce statistical fluctuations of the collision events by allowing the collisions of test nucleons that come from different ensembles. Under such a circumstance, the collision probability is reduced, In our case, the NN scattering probabilities are very small within one time step, therefore instead of evaluating probabilities of all possible collisions of test nucleons, we divide randomly the test nucleons that are available for scattering around the lattice site α into many pairs for scattering, and accordingly amplify the corresponding scattering probabilities, which is a common scheme when allowing the scattering of test nucleons from different ensembles [74,76]. The amplified scattering probabilities are given by where N α is the number of test nucleons that contribute to lattice site r α , and N ′ α is the number of test nucleons available for scattering. Since we choose a finite range form factor for coordinate in the LBUU framework, one test nucleon can be involved in different collision events from different lattice sites. Those test nucleons that have already collided on other lattice site are excluded in the scattering of the present lattice site, therefore N ′ α is not necessarily equal to N α . The time step ∆t needs to be sufficiently small to pin down the effect of such an exclusion by suppressing the chance of multi-scattering attempts, as well as to keep P ′ ij less than unity. In the present LBUU framework, we choose ∆t = 0.2 fm/c for the full LBUU calculations, and 0.4 fm/c for the Vlasov calculations (i.e., the LBUU calculations without the NN scatterings).
To verify the accuracy of the stochastic collision treatment within the present LBUU framework, we simulate collisions of nucleons confined in a cubic box of volume V = 10 × 10 × 10 fm 3 with periodic boundary condition. In this simulation, we ignore the nuclear mean-field potential and the quantum nature of nucleons. Initially, 80 neutrons and 80 protons are uniformly distributed over the box, corresponding to the nucleon density of ρ = 0.16 fm 3 . Their momenta are generated according to the relativistic Boltzmann distribution, with m = 939 MeV the free nucleon mass. Here, the temperature T is taken to be 14.24 MeV so that the system has the same kinetic energy density as the zero-temperature isospin-symmetric Fermi gas of nucleons.
Using the NN elastic scattering cross section in free space [75], we simulate the time evolution of such a system up to 1 fm/c with the time-step of 0.2 fm/c and N E = 1000. It is constructive to see the collision rate as a function of the c.m. energy √ s of the colliding nucleon pair. The √ s distributions of the collision rates for np and nn plus pp are shown as red circles in the left and right panels of Fig. 1, respectively. Theoretically, considering two species of particles following relativistic Boltzmann distributions, the √ s distributions of their collision rate can be derived as where K n is the n-th order modified Bessel function, N i (N j ) is the number of particle i(j) in the volume V , and σ is their scattering cross section. The expected distributions are shown as the black solid lines in Fig. 1   Given the quantum nature of nucleons, we treat Pauli blocking in the LBUU method as follows. If the NN scattering between the i-th and j-th test nucleons happens on the lattice site r α according to P ij or P ′ ij , the directions of their final momenta p 3 and p 4 are determined according to the differential cross-section given in Ref. [75], and then the Pauli blocking factor [1 − f ( r α , p 3 )] × [1 − f ( r α , p 4 )] is used to determine whether the collision is blocked by the Pauli principle. The distribution function f τ ( r α , p) is calculated according to Eq. (28). For the momentum-space volume ∆ 3 p i , we adopt a sphere with radius R p τ ( r α , p) centered at p i . In typical BUU transport models, R p τ ( r α , p) is a constant of about a hundred MeV. For the calculation of small-amplitude nuclear collective dynamics near the ground state, a specifically proposed R p τ ( r α , p) is more suitable [51], i.e., where ∆p is a constant that should be sufficiently small, and p F τ = (3π 2 ρ τ ) 1/3 is the nucleon Fermi momentum.

Ground state initialization and evolution stability
In the present LBUU method, we obtain the ground state of nuclei at zero temperature through varying the Hamiltonian with respect to the nuclear radial density, which is sometimes called as Thomas-Fermi (TF) initialization [66,77,51,78] in one-body transport model. We assume that for a ground state nucleus at zero temperature, its Wigner function satisfies where p F τ ( r) is local Fermi momentum given by It should be noted that in principle, with the inclusion of NN scatterings, which goes beyond mean-field correlations, the nucleon momentum distribution in the ground state may differ slightly from the zerotemperature Fermi distribution. If we assume that the nucleus is spherical for simplicity, the total energy of a ground state nucleus at zero temperature can be regarded as a functional of radial density ρ τ (r) and its spatial gradients, We can obtain the neutron (proton) radial density in a ground state nucleus by varying the total energy with respect to ρ τ (r) [note that for protons the contribution from the Coulomb interaction in Eq. (21) should also be included in the Hamiltonian density] as, where µ τ is the chemical potential of proton or neutron inside the nucleus and its value is determined by the given proton number Z or neutron number N. The quantity U τ p F τ ρ τ (r) , r refers to the single nucleon potential of the nucleon with local Fermi momentum. The single nucleon potential is derived by varying the Hamiltonian density in Eqs. (6)- (12) with respect to the phase space distribution function and density gradients, and its detailed expression for the N3LO Skyrme pseudopotential is shown in Ref. [73]. The physical significance of Eq. (41) is very intuitive: in a classical picture, in a ground state nucleus at zero temperature, the nucleons in the Fermi surface of different radial position have the same chemical potential. The local density ρ τ ( r) for a ground state spherical nucleus is obtained by solving Eq. (41) subject to the following boundary condition on the total local density ρ(r) = ρ n (r) + ρ p (r), Here, r B is the boundary of the nucleus and it satisfies ρ(r B ) = 0.
In the present LBUU framework, the initial coordinates of test nucleons are generated according to the obtained ρ τ ( r), while their initial momenta to a zero-temperature Fermi distribution with the Fermi momentum given in Eq. (39). Due to the presence of the form factor S( r − r ′ ) introduced in Eq. (4), the density is smeared slightly in the LBUU calculations compared with the realistic local density. Thus the initial ground state radial density distribution is slightly different from the solution of Eq. (41). Unlike the Gaussian wave packet which is used to mimic the Wigner function in QMD model [57], the form factor S( r − r ′ ) does not have any physical meaning, and it can be regarded as a numerical technique introduced in the test-particle approach so that one can obtain well-defined densities and mean fields. As shown in the following, an additional gradient term in the local density can compensate for the effects caused by the smearing of the local density due to the form factor. In the following of this subsection, we will usẽ ρ( r) to represent the local density in the LBUU calculation while ρ( r) to the realistic local density. The local densityρ( r) can be regarded as a convolution of the realistic local density with the form factor, To express ρ( r) in terms ofρ( r), we have formally where we have truncated at next-to-leading order [the ∇ρ( r) term vanishes because of the symmetry of the integral], and S −1 ( r − r ′ ) is the inverse of S( r − r ′ ) which satisfies The parameter c defined by is a small constant that only depends on the form of S. In the LBUU framework, to obtain ρ( r) by the direct correction onρ( r), is not feasible since numerically the density in Eq. (44) is not always positive. If we substitute Eq. (44) into the total Hamiltonian, with several necessary approximations, we obtain an additional term that is proportional to cρ( r)∇ 2ρ ( r). This term leads to an additional gradient term E [2] ∇ 2ρ in the equations of motion Eq. (20). Therefore, in practice we can add the additional gradient termsẼ [2] ∇ 2ρ to the equations of motion, to compensate for the smearing of density due to the form factor. In principle, the parameterẼ [2] should contain higher order effects, therefore we adjust it to roughly obtain the ground state rms radius evolution with the smallest oscillation, since the rms radius in the exact ground state should not change with time. NormallyẼ [2] is a small parameter around 15 MeV for various (N3LO) Skyrme parameter sets. It should be mentioned that this correction on the density gradient term only improves the stability of the ground state evolution (rms radius and radial density profile) slightly, and does not cause much difference on the results for collective motions. In ideal cases with N E → ∞ and l x , l y , l z → 0, the local density in LBUU calculation will approach the physical local density, andẼ [2] will become zero. Since all the LBUU calculations are based onρ( r), we do not distinguishρ( r) and ρ( r), and ρ( r) should be interpreted asρ( r) in the rest of the article.
We first examine the ground-state evolution stability of the LHV calculation, i.e., the LBUU calculation without the collision term, since in principle all NN scatterings should be blocked in the ground state. We show in Fig. 2 the time evolution of the radial density profile from the LHV calculation for the nucleus 208 Pb in ground state up to 1000 fm/c, obtained with N E = 10000 and a time step ∆t = 0.4 fm/c by using the N3LO Skyrme pseudopotential SP6m. We notice from Fig. 2 that the profile of the radial density exhibits only very small variations with time, which indicates the success of the above initialization method. It also shows that the smearing of the local density caused by the inclusion of the form factor S does not affect the dynamic evolution significantly. Such features indicate that the present LBUU method of solving the BUU equation can be used to study long-time nuclear processes such as nuclear spallation and heavy-ion fusion reactions.  Figure 2. Time evolution of the radial density profile of the ground state of 208 Pb based on the LHV calculation (i.e., the LBUU calculation without NN scatterings) with N3LO Skyrme pseudopotential SP6m up to 1000 fm/c. Taken from Ref. [61] with permission from the American Physical Society.
Apart from the radial density profile, other properties concerning the ground-state evolution stability are also examined. In Fig. 3, we present the time evolution of the rms radius, the fraction of bound nucleons, and the binding energy of the the LHV calculation (i.e., the LBUU calculation without NN scatterings). The calculations are performed with time step ∆t = 0.4 fm/c, and N E = 5000 and 10000, respectively. The test nucleons whose form factor do not overlap with that of others are considered as free test nucleons, and they are excluded in calculating the fraction of bound nucleons and the rms radius. We notice from Fig. 3(a) that although in N E = 5000 case, the rms radius starts to decrease after about 800 fm/c, the LHV calculation gives a fairly stable time evolution of rms radius. This decrease is due to the evaporation of test nucleons from the bound nuclei, which is illustrated in Fig. 3(b). Such an evaporation of test nucleons is inevitable in transport model calculations due to the limited precision in the numerical realization, whereas it can be suppressed by increasing N E , as can be seen in Fig. 3(b), though the result with E E = 5000 is already satisfactory [61]. As shown in Fig. 3(c), the LH method ensures the energy conservation to a very high degree. The difference between the cases N E = 5000 and N E = 10000 is mainly caused by the numerical precision of the gradient term in the Hamiltonian. It is seen from Fig. 3 that the present LBUU framework can give a fairly stable ground-state time evolution. Due to the high efficiency of the GPU parallel computing, it becomes possible to include more ensembles or test particles in the LBUU calculation. As one will see in the following, to obtain the correct GDR width, as much as 30000 ensembles are needed in the full LBUU calculation with NN scatterings.  For the stability of the ground-state evolution in the full LBUU calculation, we note that for σ free NN with N E = 30000, the rms radius and the ground-state energy of 208 Pb vary less than 3.6% (0.2 fm) and 3.2% (50 MeV), respectively, during the time evolution of 0 − 500 fm/c [55]. The stability of the rms radius in the full LBUU calculation is not as good as in the LHV case, and this may be caused by the fact that with the inclusion of NN scatterings, i.e., beyond mean-field correlations, the nucleon momentum distribution in the ground state may differ slightly from the zero-temperature Fermi distribution. Apart from this, although the LH method can conserve the energy almost exactly for the mean-field evolution without NN collisions, the non-perfect energy conservation in the LBUU calculation might be caused by the NN scattering processes, which usually violate the energy conservation when the momentum-dependent meanfield potentials are used. Both problems require further investigation of the transport model calculations in the future.

Nuclear giant resonances within transport models
We consider a small excitation of the Hamiltonian whereQ is the excitation operator for a given mode, and λ is the initial excitation parameter supposed to be small. Within linear response theory [79], the response of the excitation operatorQ as a function of time is given by where |0 is the unperturbed nuclear ground state with energy E 0 , |0 ′ is the nuclear state after the perturbation, and |F is the energy eigenstate of the excited nucleus with eigen-energy E F . The strength function, which is defined as can be expressed as a Fourier integral of ∆ Q (t) in Eq. (48) by By evaluating the time evolution of ∆ Q (t) within the transport model, we can obtain the strength function, and subsequently other quantities like the peak energy, width and energy-weighted sum rules. The time evolution of ∆ Q (t) can be expressed in terms of the Wigner function f ( r, p) as follows.
If we assumeQ is a one-body operator, then it can be written as the sum of single particle operatorŝ q acting on each nucleon,Q = A iq , and the expectation value ofQ for a given state is evaluated as follows, where we have added two identity operators. Considering the definition of the one-body density matrix, ρ( r 1 , r ′ 1 ) = A r 1 r 2 · · · r N |Φ Φ| r ′ 1 r 2 · · · r N d 3 r 2 · · · d 3 r N , and combining it with the one-body operator conditionQ = A iq , we can rewrite Eq. (51) as The density matrix can be expressed in coordinate space as the inverse Fourier transform of f ( r, p) by In the above equation we have changed the integration variables, r 1 = r + s 2 and r ′ 1 = r − s 2 . We define the Wigner transform ofq in coordinate space, where q r+ s 2 , r− s 2 = r+ s 2 |q| r− s 2 represents the matrix element ofq in coordinate space. Substituting Eq. (53) and the inverse transform of Eq. (54) into Eq. (52), the expectation ofQ can be written into the following form, which means the time evolution of Q can be calculated through the time evolution of f ( r, p).
In the transport model, different external excitation λQδ(t − t 0 ) can be generated by changing the positions and momenta of the test nucleons as follows [47]: The detailed forms of q( r i , p i ) for different collective modes and their corresponding initialization in the transport model will be given later.

LATTICE HAMILTONIAN VLASOV CALCULATIONS
In this section we compare the peak energy of nuclear giant resonances obtained from LBUU calculations without the NN scatterings, i.e., the LHV calculations, with that from the RPA, since the 2p-2h correlation is absent in both cases. Both the isoscalar monopole and the isovector dipole modes of 208 Pb are examined.

Isoscalar monopole mode
Since the isoscalar giant monopole resonance (ISGMR) provides information about the nuclear matter incompressibility [80,81,82,83,84,85], which is a fundamental quantity that characterizing the EOS of symmetric nuclear matter, it is interesting to study the ISGMR within the transport model to make a cross check with the incompressibility extracted from the HICs.
From point of view of the one-body transport model, the isoscalar monopole mode is regarded as a compressional breathing of the nuclear fluid. The excitation operator for the isoscalar monopole modê Q ISM and its one-body operatorq ISM take the form Through Eq. (54), we obtain the Wigner transform ofq ISM as According to Eq. (56), we can generate in the transport model the initial isoscalar monopole excitation by changing the initial phase space information of test nucleons with respect to that of the ground state, The spatial coordinates of the test nucleons remain unchanged since q ISM in Eq. (58) is independent of momentum. Note that the rms radius of a nucleus, shown in Fig. 3, is given by the square root of the expectation value ofQ ISM . The results correspond to three N3LO Skyrme pseudo potentials, SP6s, SP6m, SP6h, and one conventional Skyrme interaction MSL1, respectively. Taken from Ref. [61] with permission from the American Physical Society.
We show in Fig. 4 the time evolution of ∆ Q ISM , i.e., the difference of the expectation value of Q ISM between the excited and the ground state from the LHV calculations. The results are from one conventional Skyrme interaction MSL1, and three N3LO Skyrme pseudopotentials, SP6s, SP6m, and SP6h. In the calculation, we set the number of ensembles N E to be 5000, and the initial excitation parameter λ to be 100 MeV · fm −1 /c. One sees from the figure that the time evolution of ∆ Q ISM , or equivalently the rms radius, displays a very regular oscillation, and the quick increase of the radius with time, which is generally seen in most BUU calculations using the conventional test particle method, does not show up here. Besides that, since the only damping mechanism in the LHV calculation is Landau damping, the amplitude of the oscillation only decreases slightly. 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 comes only from one-body dissipation, since the coupling to more complex states, like 2p-2h states, is missing in RPA [86]. We obtain the peak energy of the GMR through Fourier transform of the time evolution of ∆ Q ISM shown in Fig. 4. The obtained peak energy is 13.8 MeV for SP6s, 13.6 MeV for SP6m, 13.9 MeV for SP6h, and 13.5 MeV for MSL1. In order to compare the result from the LHV calculation with that from RPA, we calculate the strength function of GMR using the Skyrme-RPA code by Colo et al. [87] with the MSL1 interaction. The obtained peak energy 14.1 MeV is comparable to that from the LHV calculation with MSL1, and the small discrepancy may reflect the difference between the semi-classical and quantum nature.

Isovector dipole mode
The isovector giant dipole resonance (IVGDR) of finite nuclei is the earliest observed nuclear collective excitation. Systematic experimental investigation of the IVGDR with photon-nuclear reactions has been done decades ago [88]. Recent precise measurements of the isovector dipole response have been performed at RCNP for 48 Ca [89], 120 Sn [90], and 208 Pb [63] with inelastic proton scattering, as well as at GSI for 68 Ni [91] by using Coulomb excitation in inverse kinematics. Recently, a low-lying mode called pygmy dipole resonance (PDR) has been observed experimentally [92,93,94,95], and this effect has already been studied based on the Vlasov equation [47]. The IVGDR [96,97,54], PDR [98,99], and electric dipole polarizability α D [100,101,102,103], which are dominated by these isovector dipole modes, provide sensitive probes to constrain the density dependence of the nuclear symmetry energy. Isovector dipole response of 208 Pb Figure 5. Same as Fig. 4 but for isovector dipole mode with λ = 25 MeV/c. Taken from Ref. [61] with permission from the American Physical Society.
For the isovector dipole mode, the external perturbation can be written in the following form, The coefficients in front of the single particle position operator are defined to keep the center of mass of the nucleus stays at rest. According to Eq. (56), the excited nucleus can be obtained in transport models by changing the initial phase space coordinates of test nucleons, We show in Fig

SPREADING WIDTH OF GIANT DIPOLE RESONANCE AND COLLISIONAL DAMPING
It is generally thought that in low-energy HICs with incident energy of only a few MeV/nucleon, the NN scattering can be safely neglected since they are mostly blocked by the Pauli principle. However, when it comes to the width of the GDR , the collisional damping caused by NN scatterings is an essential mechanism to enhance the insufficient GDR width obtained through the pure Vlasov calculation [55]. Nevertheless, to properly implement the damping mechanism caused by NN scatterings in transport models requires a rather accurate treatment of the Pauli blocking, which is a challenge to transport model calculations. The difficulty mainly comes from calculating local momentum distributions f τ ( r α , p) accurately in transport models. The inaccuracy of f τ ( r α , p) affects the accuracy of the Pauli blocking and leads to spurious collisions, which enhance the collisional damping and thus overestimate the width of nuclear giant resonances. There are three main origins of the inaccuracy of the calculated f τ ( r α , p) and thus the spurious collisions in transport models: 1) Fluctuations in calculating f τ ( r α , p) through Eq. (28) caused by too small N E ; 2) Spurious temperature caused by a finite ∆p in calculating f τ ( r α , p) (also see Ref. [51]); 3) The finite lattice spacing l causes the diffusion in local momentum space due to the average of different local lattice densities in the nuclear surface region.
In order to obtain the spreading width with high accuracy with the BUU equation, one should choose a large N E together with the sufficiently small l and ∆p. After a careful test, it is found [55] that to get a convergent GDR width, l should be smaller than 0.5 fm, ∆p smaller than 0.05 GeV, and N E larger than 30000. Further reducing ∆p and l or increasing N E leads only to a negligible decrease of the calculated GDR width. Therefore, in the following full LBUU calculations of the GDR width, we take l = 0.5 fm, ∆p = 0.05 GeV, and N E = 30000.
The collisional damping or NN scattering can have a significant effect on the width of nuclear giant resonances. Shown in Fig. 6 is the time evolution of isovector dipole response ∆ Q IVD of 208 Pb and its strength function from the LHV calculation and the full LBUU calculation with the free NN elastic scattering cross section [55]. In both cases, the N3LO Skyrme pseudopotential SP6h is adopted, and the same initial excitation with λ = 15 MeV/c is employed (we note varying λ by 2/3 leads to almost the same value of the GDR width). Shown in the left window of Fig. 6 by the dotted line is the time evolution of the expectation value 0|Q IVD |0 in the ground state of 208 Pb as obtained from the LBUU calculation with the free-space NN cross section. The expectation value 0|Q IVD |0 in the ground state of 208 Pb is negligible compared with that in the GDR cases with and without NN scatterings. It is seen that including NN scatterings enhances the damping of the oscillations significantly , and leads to a much larger width. From the strength functions, the obtained GDR width of 208 Pb is 1.5 MeV in the Vlasov calculation and 6.5 MeV in LBUU calculation with NN scatterings. We also notice from the right panel of Fig. 6 that the peak shifts to higher energy when we include the NN scatterings. The impact of NN scatterings on the width indicates that they might also affect some particular observables in low-energy HICs, e.g., the nuclear stopping of HICs in Fermi energy region, which can be studied further. Recent experiments of the 208 Pb( p, p ′ ) reaction performed at RCNP [63] have measured the GDR width of 208 Pb accurately, and a value of 4.0 MeV has been extracted. Therefore, the LBUU calculation with the free NN elastic cross section (which predicts a GDR width of 6.5 MeV) significantly overestimates the GDR width of 208 Pb. It is well known that the NN elastic cross section is suppressed in the nuclear medium, therefore the overestimation of the GDR width with σ free NN is understandable since the medium effects on the NN elastic cross section will weaken the collisional damping and consequently result in smaller GDR width. As shown in Ref. [55], in order to reproduce the experimental GDR width of 208 Pb of RCNP, a strong medium reduction of the NN cross section is needed. There are many parameterizations for the medium reduction of the NN cross section [104,105,106,107], which could be density-, collision energy-, isospin-dependent. As an example, we choose the FU4FP6 parameterization [107] for the medium reduction of the cross section to calculate the strength function and width of the GDR in 208 Pb.The FU4FP6 parameterization of the medium reduction is density-, momentum-and isospindependent, and it is preferred by the nucleon induced nuclear reaction cross section data [107] and predicts a very strong in-medium reduction of NN scattering cross sections. The strength function of the GDR in 208 Pb from the LBUU calculation is shown in Fig. 7 and compared with the RCNP data [63]. The obtained GDR width from the LBUU calculation through the full width at half maximum (FWHM) of the strength function is 4.32 MeV, which is consistent with 4.0 MeV from the RCNP experiment. However, the Skyrme pseudopotential SP6h adopted in the calculation overestimates the peak energy by about 1.5 MeV. Using effective interactions with a different symmetry energy slope parameter L or nuclear effective masses can easily reproduce the correct peak energy [54]. In order to compare the shape of the strength function and the value of the width of the GDR in 208 Pb, we shift the strength function from LBUU calculation to meet the experimental peak energy. We conclude from Fig. 7 that the present LBUU method with the FU4FP6 parameterization [107] for the medium reduction of the NN scattering cross section can nicely reproduce the measured shape of the strength function and the width of GDR in 208 Pb. It should be stressed that the FU4FP6 parameterization suggests a very strong in-medium reduction of NN scattering cross sections, consistent with the conclusions obtained in Ref. [55].  Figure 7. Strength function of the GDR in 208 Pb after a perturbationĤ ex = λQ IVD δ(t − t 0 ) with λ = 15MeV/c from the LBUU calculation using the FU4FP6 parameterization [107] for the in-medium NN scattering cross section. The strength function measured in the RCNP experiment [63] is also shown for comparison.

SUMMARY AND OUTLOOK
We have reviewed the recent progress in calculating the nuclear collective motions by solving the BUU equation with the LH method. In order to calculate the nuclear collective motions accurately with the BUU equation, the present LBUU framework contains the following features: 1) the smearing of the local density is considered in the equations of motion self-consistently through the lattice Hamiltonian method; 2) the initialization of a ground state nucleus is carried out according to a nucleon radial density distribution obtained by varying the same Hamiltonian that governs the evolution; 3) the NN collision term in the BUU equation is implemented through a full-ensemble stochastic collision approach; 4) the high-performance GPU parallel computing is employed to increase the computing efficiency. The present LBUU framework with these features brings us to a new level of precision in solving the BUU equation.
Within the LBUU framework, it has been shown that the peak energies of ISGMR and IVGDR from the pure Vlasov calculation are consistent with the RPA calculation, and the full LBUU calculation is able to give reasonable GDR strength function compared with the experimental data. The peak energies can be used to extract information about the nuclear EOS, while the width of the GDR can constrain the medium reduction of the elastic NN scattering cross section. The success of the present LBUU framework in describing the nuclear collective motions has demonstrated its capability in treating the stability of ground-state nuclei and the nuclear dynamics near equilibrium. Therefore the present LBUU framework provides a solid foundation for studying the long-time process of heavy-ion reactions at low energies, e.g., heavy-ion fusion and multi-nucleon transfer reactions at near-barrier energies, based on solving the BUU equation. The significant effects of the collisional damping on the width of the nuclear GDR indicate that NN scatterings should play a crucial role in nuclear collective dynamics with small amplitude oscillations.
The present LBUU framework has been shown to significantly reduce the uncertainties of the transport model simulations for HICs in various aspects, especially for the stability of the nuclear ground state evolution and the very accurate treatment of NN scatterings as well as the Pauli blocking. This is rather important for various studies of HICs based on transport model calculations, e.g., the extraction of the nuclear EOS and the in-medium NN scattering cross sections. Further studies of HICs from low to intermediate energies within the present LBUU framework are in progress and it is expected that more reliable information on the nuclear EOS, the in-medium NN scattering cross sections, and the effective nuclear interactions could be extracted in near future.