Relativistic Coulomb excitation within Time Dependent Superfluid Local Density Approximation

Relativistic Coulomb excitation within Time Dependent Superfluid Local Density Approximation

I. Stetcu Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA    C. A. Bertulani Department of Physics and Astronomy, Texas A & M University - Commerce, Commerce, TX 75429, USA    A. Bulgac Department of Physics, University of Washington, Seattle, WA 98195–1560, USA    P. Magierski Department of Physics, University of Washington, Seattle, WA 98195–1560, USA Faculty of Physics, Warsaw University of Technology, ulica Koszykowa 75, 00-662 Warsaw, POLAND    K.J. Roche Department of Physics, University of Washington, Seattle, WA 98195–1560, USA Pacific Northwest National Laboratory, Richland, WA 99352, USA
August 17, 2019

Within the framework of the unrestricted time-dependent density functional theory, we present for the first time an analysis of the relativistic Coulomb excitation of the heavy deformed open shell nucleus U. The approach is based on Superfluid Local Density Approximation (SLDA) formulated on a spatial lattice that can take into account coupling to the continuum, enabling self-consistent studies of superfluid dynamics of any nuclear shape. We have computed the energy deposited in the target nucleus as a function of the impact parameter, finding it to be significantly larger than the estimate using the Goldhaber-Teller model. The isovector giant dipole resonance, the dipole pygmy resonance and giant quadrupole modes were excited during the process. The one body dissipation of collective dipole modes is shown to lead a damping width MeV and the number of pre-equilibrium neutrons emitted has been quantified.

25.70.De, 21.60.Jz, 25.70.-z, 24.30.Cz
preprint: LA-UR-14-21210, NT@UW-14-06

Coulomb excitation represents an ideal method to probe the properties of large amplitude nuclear motion, because the excitation process is not obscured by uncertainties related to nuclear forces. The excitation probabilities are governed by the strength of the Coulomb field only and they can be fully expressed in terms of the electromagnetic multipole matrix elements Winther and Alder (1979); Baur and Bertulani (1986); Emling (1994); *abe; Bertulani and Ponomarev (1999); Lanza et al. (1999); Hussein et al. (2004). From the theoretical point of view, Coulomb excitation can be treated as a textbook example of a nuclear system being subjected to an external, time-dependent perturbation. However, in order to be able to probe nuclear collective modes involving multi-phonon states for example Boretzky et al. (2003); Ilievski et al. (2004), a large amount of energy has to be transferred to the nuclear system. Thus the interaction time should be relatively short and the velocity of the projectile has to be sufficiently large for an efficient excitation of nuclear modes of frequency , the collision time has to fulfill the condition that . Here is the impact parameter, is the projectile velocity, and is the Lorentz factor. One of the best known examples of collective nuclear motion is the isovector giant dipole resonance (IVGDR). A reasonably good estimate of the IVGDR vibrational frequency is for spherical nuclei. It implies that the excitation of a heavy nucleus to such energies requires a relativistic projectile.

We report on a new and powerful method to study relativistic Coulomb excitation and nuclear large amplitude collective motion in the framework of Time Dependent Superfluid Local Density Approximation (TDSLDA). This is a fully microscopic approach to the problem based on an extension of the Density Functional Theory (DFT) to superfluid nuclei and time-dependent external probes, where all the nuclear degrees of freedom are taken into account on the same footing, without any restrictions and where all symmetries (translation, rotation, parity, local Galilean covariance, local gauge symmetry, isospin symmetry, minimal gauge coupling to electromagnetic (EM) fields) are correctly implemented Bulgac (2013); Stetcu et al. (2011). The interaction between the impinging U projectile and the U target is very strong (, where is the fine structure constant), which thus require a non-perturbative treatment, and the excitation process is highly non-adiabatic. We assume a completely classical projectile straight-line motion since its de Broglie wavelength is of the order of fm for . In evaluating the EM-field created by the uranium projectile with a constant velocity along the z-axis, we neglect its deformation. The projectile produces an EM-field described by scalar and vector Lienard-Wiechert potentials. These fields couple to a deformed U target nucleus residing on a spatial lattice, see Ref. sup . The interaction leads to a CM motion of the target as well as to its internal excitation and full 3D dynamical deformation of the target. In order to follow the internal motion for a long enough trajectory that allows the extraction of useful information, we perform a transformation to a system in which the lattice moves with the CM. The required transformation for each single particle wave function reads , with describing the CM motion and the momentum operator. The equation of motion (simplified form here) for becomes


where is the CM velocity and the total current density.

The target nucleus is described within the SLDA and its time evolution is governed by the TD meanfield-like equations (spin degrees of freedom are not shown):


The single-particle Hamiltonian and the pairing field are obtained self-consistently from an energy functional that is in general a function of various normal, anomalous, and current densities. The external electromagnetic (EM) field has the minimal gauge coupling (and similarly for the time-component) in all terms with currents, as well as in the definition of the momentum operator in Eq. (1), details in sup . In the current calculation, the Skyrme SLy4 energy functional Chabanat et al. (1998a); *Sly4E was adopted, with nuclear pairing as introduced in Ref. Yu and Bulgac (2003), which provides a very decent description of the IVGDR in U Stetcu et al. (2011). The coupling between the spin and the magnetic field was neglected. The Coulomb self-interaction between protons of the target nucleus is taken into account using the modification of the method described in Ref. Castro et al. (2003), so as not to include contributions from images in neighboring cells. For the description of the numerical methods see Refs. Bulgac and Roche (2008); Bulgac et al. (2011) and many other technical details can be found in sup .

The DFT approach to quantum dynamics has some peculiar characteristics. Unlike a regular quantum mechanics (QM) treatment one does not have access to wave functions, but instead to various one-body densities and currents. Within a DFT approach quantities trivial to evaluate in QM become basically impossible to calculate. For example, by solving the Schrödinger equation one can evaluate at any time the probability that a system remained in its initial state from , where is the solution of the Schrödinger equation (or some of its approximations). Within DFT one has access to the one-body (spin-)density and one-body (spin-)current with no means to compute the probability . One can calculate for example a quantity such as , but there is no obvious way to relate it to the probability . One might try to define instead through the overlap of the initial and current “Slater determinants” constructed through the fictitious single-particle wave functions entering the DFT formalism, which is a rather arbitrary postulate. One can find quite often in literature various formulas used within DFT treatment of nuclei, which are simply “copied” from various mean field approaches, without any solid justification provided. Restrictions inherent to a DFT approach, prevent us from being able to calculate various quantities, which within a QM approach are easy to evaluate. Even though we evaluate accurate densities and currents well beyond the linear regime, within a DFT approach we cannot separate for example the emission of one and two photons from an excited nucleus, which however could be estimated relative easily within a perturbative linear response approach such a (Q)RPA. On the other hand a DFT approach has unquestionable advantages, allowing us to go far into the non-linear regime and describe large amplitude collective motion.

Figure 1: (color online) The emitted energy rate via EM radiation for a collision with impact parameter fm, for three orientations. In two cases the nuclear symmetry axis is parallel to the reaction plane and perpendicular (dot-dashed line) or parallel (dashed line) with respect to the incoming projectile, while in the third it is both perpendicular to the reaction plane and the incoming projectile (dotted line). These configurations are denoted by , and , respectively. We show time-averaged quantities, while in the insert, for one configuration, we also show the raw, strongly oscillating, data. The rate at which this quantity changes is directly related to the characteristic damping time, which we estimate at , leading to a width MeV.

The incoming projectile excites various modes in the target nucleus and the axial symmetry of the initial ground state is lost. Because U is highly deformed the energy of the first is 45 keV, corresponding to a very long rotational period, and thus during simulation time considered here ( sec.) it can be considered fixed. The identification of these modes requires certain care, since during the collision the system is beyond the linear regime and the analysis using the response function is not applicable in general. However, the information about the excited nuclear modes is carried in the subsequent EM radiation leading to nuclear de-excitation. De-excitation to the ground state via photon emission requires times of about sec., which is four orders of magnitude longer than in the current calculations. However, it is possible to compute the spectrum of the pre-equilibrium neutrons and gamma radiation, which allows the identification of the excited nuclear modes. We can accurately treat the system as a classical source of electromagnetic radiation and the time dependence of the proton current governs the rate of emission, see Refs. sup ; Baran et al. (1996); Oberacker et al. (2012):


with . Here , is the spherical Bessel function of order , and is the proton current. The emission rate is plotted in Fig. 1. The magnitude of this quantity indicates that the total amount of radiated energy during the evolution time (about fm/c) is rather small compared to the total absorbed energy and does not exceed MeV, which is about of the deposited energy reported in Table 1 below. This implies that the effect of damping of nuclear motion due to the emitted radiation can be neglected for such short time intervals. Consequently, the decreasing intensity of the radiation, see Fig. 1, is merely related to the rearrangements of the intrinsic structure of the excited nucleus caused by damping of collective modes due to the one-body dissipation mechanism. It has to be emphasized that within the framework of the presented approach one is able to extract only a small fraction of the spreading width , which is due to the one-body dissipation mechanism. The two-body effects require e.g. stochastic extension of TDSLDA which would allow for a dynamic hopping between various mean-fields, and thus could account for collisional damping as well.

TDSLDA provides the EM power spectrum sup ; Baran et al. (1996); Oberacker et al. (2012), arising from the multipole expansion in Eq. (4). This quantity is different from what one would compute within a linear response approach or first order perturbation theory, see, e.g., Refs. Winther and Alder (1979); Baur and Bertulani (1986); Emling (1994); *abe; Bertulani and Ponomarev (1999); Lanza et al. (1999); Hussein et al. (2004), which provides the excitation probability only , where is the transition density and - the external field. is proportional to the excitation probability, here in the non-linear regime, and the subsequent photon emission probability as well. A typical example of the emitted EM radiation for a given impact parameter is shown in Fig. 2, here due to the internal excitation of the system alone. The EM radiation due the CM motion has been calculated separately (see Table 1 and sup ).

Figure 2: (color online) (a) The total energy spectrum (solid line) of emitted EM radiation, averaged over the target-projectile configurations, at the impact parameter fm. We show the total quadrupole contribution (double-dotted line), as well as the contributions from the three target-projectile orientations using the same symbols as in Fig. 1. (b) The radiation emitted from the target nucleus when only the dipole component of the projectile electromagnetic field is used. The insert shows the pygmy resonance contribution to the emitted spectrum, visible in the main figure as the slope change at low energies.

In Fig. 2(a) the emitted radiation shows a well defined maximum at energy MeV which corresponds to the excitation of IVGDR. We have applied a smoothing of the original calculations using the half-width of MeV. Therefore, the original separate peaks split due to the deformation of U merge into a single wider peak. However, at larger frequencies another local maximum exists which we associate with the isovector giant quadrupole resonance (IVGQR). In order to rule out other possibilities we have repeated the calculation by retaining only the dipole component of the electromagnetic field produced by the projectile sup . The results are shown in Fig. 2(b). In this case, the high-energy structure above 20 MeV disappears, evidence that the high energy peak is related to the IVGQR. Noticeable contribution to the total radiation is coming from the quadrupole component of radiated field.

At low energies a change of slope occurs at about MeV, present at the same energy for all impact parameters and orientations, see Ref. sup . It indicates a considerable amount of strength at low energies, giving rise to an additional contribution to the EM radiation. We attribute this additional structure to the excitation of the pygmy dipole resonance (PDR). The inset of the figure 2 shows the spectrum of emitted radiation due to this mode. The contribution to the total radiated energy coming from the PDR is rather small and reads: , , and keV for impact parameters , , and fm, respectively. It corresponds to about %, %, %, % of the emitted radiation (due to internal motion) respectively. The relatively small amount of E1 strength obtained in our calculations, in the region where the PDR is expected, agrees with recent measurements Hammond et al. (2012).

The comparison between the average energy transferred to the internal motion of the target nucleus for three values of the impact parameter obtained within TDSLDA and within a simplified Goldhaber-Teller model Goldhaber and Teller (1948) presented in Table 1 shows that significantly more energy is deposited by the projectile within the TDSLDA. The Goldhaber-Teller model is equivalent to a linear response result, assuming that all isovector transition strength is concentrated in two sharp lines, corresponding to an axially deformed target. An exact QRPA approach would therefore severely underestimate the amount of internal energy deposited, one reason being the non-linearity of the response, naturally incorporated in TDSLDA. Another reason is the fact that the present microscopic framework describing the target allows for many degrees of freedom to be excited, apart from pure dipole oscillations. At the same time, the CM target energy alone is approximately the same as obtained in a simplified point particles Coulomb recoil model of both the target and projectile.

12.2 39.29 0.668 0.911 0.960 17.58 24.68
14.6 19.2 0.608 0.567 0.963 10.32 14.51
16.2 12.87 0.547 0.411 0.963 7.41 10.43
20.2 5.41 0.444 0.199 0.961 3.43 4.84
12.2 25.11 0.588 0.5 0.941 12.94 18.17
14.6 13.16 0.498 0.306 0.942 7.22 10.16
16.2 8.97 0.470 0.217 0.939 5.02 7.07
20.2 3.8 0.367 0.106 0.934 2.16 3.05
12.2 24.21 0.591 0.407 0.930 12.36 17.33
14.6 12.58 0.513 0.245 0.929 6.65 9.34
16.2 8.5 0.464 0.175 0.926 4.49 6.32
20.2 3.5 0.353 0.085 0.919 1.78 2.51
Table 1: Internal excitation energy in TD-SLDA () and in the Goldhaber-Teller model (), as well as the EM energy radiated () from the excited nucleus during time interval fm/c after collision, for four values of impact parameters and three orientations of the nucleus with respect to the beam. We also list their respective ratios to the total transferred energy. Finally, the Goldhaber-Teller model results () for effective mass are presented in the last column. All energies are in MeV.

The average energy radiated due to the internal excitation represents only part of the total radiated energy. (One should remember that a straightforward DFT approach provides no measure for the variance.) Also, because of the spreading of the strength due to one-body dissipation only a fraction of the energy (where is the EM-width alone and the total width of the IVGDR) is emitted as a pulse, as shown in Fig. 1. A subsequent pulse of reduced amplitude is to be expected after a delay fm/c. Since our simulations times are much shorter we are not able to see emission of the second photon, as reported in experiment Boretzky et al. (2003); Ilievski et al. (2004), where two photons were measured in coincidence. In our calculations we have followed the nuclear evolution during approximately fm/c after collision. The other component of the EM radiation arises from the CM acceleration as a result of collision (Bremsstrahlung), during the relatively short time interval . The radiation emitted from the internal motion has a much longer time scale.

We can estimate the cross section for the emission of radiation by assuming that the asymptotic transition probability for a given impact parameter is given by . Here , and are the total energies transferred to the target nucleus during the collision at the impact parameter and for the three independent orientations. Our simulation yields mb. A detailed comparison of intensities of radiation for various impact parameters and orientations is shown in Table 1. It is evident that although the intensity of radiation decreases with increasing impact parameter, the ratio between the intensities due to the internal modes with that of the CM motion remains fairly constant and depends slightly on orientation For the orientation perpendicular to the beam and parallel to the reaction plane the target nucleus is the most efficiently excited which results in a larger ratio.

Figure 3: (color online) The number of neutrons inside the sphere of radius fm around the target nucleus as a function of time for the four impact parameters. We use the same convention as in Fig. 1 for the possible orientations. The emission rate is inverse proportional with the value of the impact parameter.

The average energy deposited in the target nucleus is of the order of the neutron separation energy. In Fig. 3 we plot the total number of neutrons inside a (smoothed) sphere of radius fm which is slightly larger than the nuclear diameter (see Ref. sup for details). For all these impact parameters neutrons can leak from the excited system. Since more energy is deposited in the nucleus with perpendicular orientation with respect to the beam, the rate of emitted neutrons is larger in that case. However, the simulation box is large enough (about 40 times bigger than the nucleus) so that during the evolution the calculations are not affected by the emitted neutrons.

We thank G.F. Bertsch for a number of discussions and reading the manuscript. We acknowledge support under U.S. DOE Grants DE-FC02-07ER41457, DE- FG02-08ER41533, NSF grant PHY-1415656, Polish National Science Centre (NCN) Grant, decision no. DEC-2013/08/A/ST3/00708, and in part by the ERANET-NuPNET grant SARFEN of the Polish National Centre for Research and Development (NCBiR). Part of this work was performed under the auspices of the National Nuclear Security Administration of the US Department of Energy at Los Alamos National Laboratory under contract No. DE-AC52-06NA25396. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. Some of the calculations reported here have been performed at the University of Washington Hyak cluster funded by the NSF MRI Grant No. PHY-0922770. KJR was supported by the DOE Office of Science, Advanced Scientific Computing Research, under award number 58202 ”Software Effectiveness Metrics” (Lucille T. Nowell).


Supplemental material to:

Relativistic Coulomb excitation within Time Dependent Superfluid Local Density Approximation

I. Stetcu, C. Bertulani, A. Bulgac, P. Magierski and K.J. Roche

August 17, 2019

Density functional in TDSLDA and coupling to electromagnetic field

Here we present various definitions and conventions which we have used in the manuscript. The density functional is constructed from the following local densities and currents:

  • density:

  • spin density:

  • current:

  • spin current (2nd rank tensor):

  • kinetic energy density:

  • spin kinetic energy density:

  • anomalous density:



The coupling of the nuclear system to the electromagnetic field:


requires the following transformation of proton densities and currents:

  • density:

  • spin density:

  • current:

  • spin current (2nd rank tensor):

  • spin current (vector):

  • kinetic energy density:

  • spin kinetic energy density:

As a result the proton single particle hamiltonian has the form:




Numerical Implementation

We build a spatial three-dimensional Cartesian grid in coordinate space with periodic boundary conditions, and derivatives evaluated in momentum (Fourier-transformed) space. This method represents a flexible tool to describe large amplitude nuclear motion as it contains the coupling to the continuum and between single-particle and collective degrees of freedom. For the present problem, we have considered a box size of with the lattice constant 1 fm. The time step has been set to with a total time interval of about . The projectile is initially placed at such a distance from the target nucleus that the collision occurs after . Even though inially the projectle is far enough from the target and hence the EM fields are weak, spurious excitations produced by a sudden switch of the EM interaction at are possible. They were avoided by multiplying the the EM potentials in Eq. (LABEL:LW) by the smoothing function , where fm, fm. This ensures that the EM field varies smoothly within the distance , but stay approximately equal to its physical value within the distance .

Coulomb potential on the lattice

Here we describe the method used to describe the Coulomb self-interaction of the target nucleus.

Consider the charge distribution :


Note that above we have defined as (note in the formula). After the Fourier transform one gets:


The above prescription generates however the spurious interaction between neighbouring cells.

Therefore we define the modified potential ( denote number of equidistant lattice points in each direction, , , is lattice constant):


Clearly the Fourier transform is:


and moreover


where in the last term is the Fourier transformed density on the lattice . In practice it means that one has to perform forward and backward Fourier transforms on the lattice three times bigger in each direction.

This may however be avoided if one realizes that the Fourier transform of the density in the larger lattice can be expressed through the Fourier transforms in the smaller lattices:


and we need to perform 27 FFTs on the smaller lattice for of the following quantities:

Subsequently we obtain the potential through the relation:


Dipole component of the electromagnetic field produced by the projectile

Figure 4: Coordinate system used to describe the reaction

Coordinates (see fig. 7):


Interaction potential:




and . is the charge density of the nucleus at location and are proton wavefunctions . The vector potential is given by


In order to extract the dipole component we used the interaction Hamiltonian:


where one subtracts the second term which is responsible for the c.m. scattering (i.e. monopole field).

Consequently the dipole term reads:


where r is the coordinate of one of the protons in the target. A sum over has to be performed, i.e.


where is the part of the interaction which does not involve the intrinsic structure of the nucleus:


Electromagnetic radiation from a nucleus described within TDSLDA

Let us consider the proton density and current (we use Gauss units):




Maxwell equations:


and spatial Fourier transforms:


Hence clearly:








Therefore one has a freedom to choose (gauge transformation) whereas is the gauge invariant part of the vector potential.

We choose the Coulomb gauge:




and only perpendicular components of electric and magnetic fields are responsible for emission of radiation. The important equation in this case is the fourth Maxwell equation:


Since the lhs represents the vector of type therefore:




Substituting the potential :




Therefore in the Coulomb gauge


and in the far zone