Selfconsistent collective coordinate for reaction path and inertial mass
Abstract
We propose a numerical method to determine the optimal collective reaction path for the nucleusnucleus collision, based on the adiabatic selfconsistent collective coordinate (ASCC) method. We use an iterative method combining the imaginarytime evolution and the finite amplitude method, for the solution of the ASCC coupled equations. It is applied to the simplest case, the scattering. We determine the collective path, the potential, and the inertial mass. The results are compared with other methods, such as the constrained HartreeFock method, the Inglis’s cranking formula, and the adiabatic timedependent HartreeFock (ATDHF) method.
pacs:
24.60.k, 24.10.Lx, 25.60.Pj, 25.70.LmI Introduction
The timedependent HartreeFock (TDHF) method has been extensively applied to studies of heavyion reaction Negele (1982); Simenel (2012); Nakatsukasa (2012); Maruhn et al. (2014); Nakatsukasa et al. (2016). The TDHF provides a successful description for the time evolution of onebody observables. Its small amplitude limit corresponds to the randomphase approximation Ring and Schuck (1980); Blaizot and Ripka (1986), which is a current leading theory for the nuclear response calculations. However, beyond the linear regime, it is not trivial to extract the quantum mechanical information from the TDHF trajectories of given initial values. It is also wellknown that the TDHF has some drawbacks due to its semiclassical nature Negele (1982); Ring and Schuck (1980). For instance, the realtime description of subbarrier fusion and spontaneous fission processes is practically impossible, because a single Slater determinant with a single average meanfield potential is not capable of describing quantum mechanical processes in rare channels.
The “requantization” of TDHF is a possible solution to these problems, that was proposed from a view point of the path integral Reinhardt (1980); Negele (1982). However, the original quantization prescription requires the identification of periodic TDHF trajectories, which is a very difficult task. As far as we know, there have been no application of the theory to realistic nuclear problems Baranger et al. (2003). A family of the periodic TDHF trajectories is associated with a collective subspace decoupled from the other intrinsic degrees of freedom. If we identify the collective subspace spanned by a small number of canonical variables, the requantization becomes much easier than finding the periodic orbits Nakatsukasa et al. (2016). In fact, the theory of the adiabatic TDHF (ATDHF) was aiming at determining such an optimum collective subspace Brink et al. (1976); Villars (1977); Baranger and Vénéroni (1978); Goeke and Reinhard (1978). The ATDHF, however, encounters a “nonuniqueness” problem, namely cannot provide a unique solution for the collective subspace. In order to uniquely fix the solution, a prescription, socalled validity condition, was proposed Reinhard and Goeke (1978). Goeke, Reinhard, and coworkers have developed a numerical recipe for the reaction path and inertial mass solving the ATDHF equations of the initialvalue problem Goeke et al. (1983); Reinhard and Goeke (1987). Their procedure requires us to calculate a large number of trajectories with different initial states, then, to obtain the optimal collective path as an envelope curve of those Goeke et al. (1983).
The selfconsistent collective coordinate (SCC) method, originally proposed by Marumori and coworkers Marumori et al. (1980), is solely based on the invariance principle of the TDHF equation in the collective subspace, which treats the collective coordinate and the momentum on an equal footing. The SCC method is able to determine the unique collective path. In addition, the AndersonNambuGoldstone (ANG) modes are properly decoupled in the SCC method Matsuo (1986); Nakatsukasa et al. (2016). Its weak point is that practical solutions to the basic equation was restricted to a perturbative expansion around the HF state. To overcome this perturbative nature of the SCC, a method treating the coordinate in a nonperturbative way but expanding with respect to momenta has been later proposed. It is named “adiabatic selfconsistent collective coordinate (ASCC) method” Matsuo et al. (2000). The ASCC provides an alternative practical solution to the SCC Matsuo et al. (2000): The state is determined at each value of by solving the equation expanded up to the second order in . The ASCC method has been successfully applied to nuclear structure problems with largeamplitude shape fluctuations/oscillations for the Hamiltonian of the separable interactions Hinohara et al. (2008, 2009, 2010, 2011, 2012); Sato et al. (2012); Matsuyanagi et al. (2016); Nakatsukasa et al. (2016). It should be noted that a solution to the nonuniqueness problem of the ATDHF was given by higherorder equations with respect to momenta Mukherjee and Pal (1982); Klein et al. (1991), which are similar to the ASCC equations.
In this paper, we apply the ASCC method to nuclear reaction studies, then, selfconsistently determine the optimal reaction path, the internuclear potential, and the inertial mass. Since the separable interactions, such as the pairingplusquadrupole interaction, are not applicable to a system with two colliding nuclei, we need to treat the Hamiltonian of a nonseparable type. For this purpose, we develop a computer code of a novel numerical technique. We use a combining procedure of the imaginarytime method Davies et al. (1980) and the finiteamplitude method Nakatsukasa et al. (2007); Avogadro and Nakatsukasa (2011, 2013) for the solution of the ASCC equations. We show that this method nicely works for the threedimensional (3D) coordinatespace representation, taking a reaction of Be as an example.
The paper is organized as follows. In Sec. II, we give the formulation of the basic ASCC equations to determine the onedimensional (1D) collective path and the canonical variables . In Sec. III, we show the numerical results and compare with those of conventional methods. Summary and concluding remarks are given in Sec. IV.
Ii Theoretical framework
ii.1 The adiabatic selfconsistent collective coordinate (ASCC) method
The SCC method aims at determining a collective submanifold embedded in the large dimensional TDHF space of Slater determinants, which is maximally decoupled from the remaining intrinsic degrees of freedom. For the 1D collective path, a pair of canonical variables along the collective path are introduced by labeling the Slater determinants as , where and respectively represent the coordinate and the conjugate momentum. Once the states are determined, the (classical) collective Hamiltonian is given by , where is the total Hamiltonian of the system. Therefore, the main task is to determine on a decoupled collective submanifold.
In the ASCC Matsuo et al. (2000), the wave function is written in a form
(1) 
using a local generator which is defined as
(2) 
The collective coordinate operator is an infinitesimal generator of “accelerating” the system. The momentum operator is introduced in the similar way as an infinitesimal generator for “translating” the system, .
Since the Thouless theorem guarantees that small variation of a Slater determinant can be generated by the particlehole (ph) excitations Ring and Schuck (1980), the local generators, and , can be written in terms of ph and hp operators as
(3)  
(4) 
In this paper, the indexes and refer to the hole and particle states with respect to , respectively, Hereafter, the creation and annihilation operators are denoted as instead of for simplicity. These generators are required to follow the weak canonicity condition
(5) 
In the ASCC, the collective momentum is assumed to be small. Keeping the expansion with respect to up to the second order, the invariance principle of TDHF equation leads to a set of ASCC equations Matsuo et al. (2000); Hinohara et al. (2007, 2008, 2010); Nakatsukasa (2012); Nakatsukasa et al. (2016) to determine the wave function and the local generators selfconsistently along the collective path. In this paper, we consider only the 1D collective motion, without taking the pairing correlations into account. The equations in the zeroth, first, and second order in momentum read, respectively,
(6)  
(7)  
(8) 
where is the “moving” Hamiltonian. The collective potential is defined as
(9) 
and is the inertial mass of the collective motion. Equation (6) is called “moving meanfield equation” (“moving HartreeFock (HF) equation”), and Eqs. (7) and (8) are “moving randomphase approximation (RPA)”.
In fact, to derive the secondorder equation (8), an additional term called “curvature term” Matsuo et al. (2000) is neglected. Although the exact treatment of the curvature term is possible, it numerically involves iterative tasks and has only minor effect on the final result Hinohara et al. (2007). Here, we neglect this curvature term, which leads to Eq. (8). Equation (6) looks similar to a constrained HartreeFock (CHF) equation. However, the constraint operator changes along the collective path , which is selfconsistently determined by the moving RPA equations (7) and (8).
Substituting and of Eqs. (3) and (4) into Eqs. (7) and (8) leads to
(10)  
(11) 
where the and matrix elements are defined as
(12) 
When all of these matrix elements are real, Eqs. (10) and (11) can be recast into an eigenvalue equation
(13) 
with
(14) 
where is the movingRPA eigenfrequency. can be pure imaginary (). The generator can be obtained from a matrix equation for ,
(15) 
Equation (13) has many solutions, among which we choose the collective mode of our interest. For instance, in numerical calculation for the scattering Be in Sec. III, the lowest quadrupole mode of excitation is selected.
Since the scale of the coordinate is arbitrary, the ASCC equations (6), (7), and (8) and the weak canonicity condition (5) are invariant with respect to the scale transformation of the collective coordinate, (). The generators and the collective inertial mass are transformed as , , and , respectively. Therefore, when we perform the numerical calculation to determine the collective coordinate , we need a condition to fix the scale of the coordinate . A convenient choice could be the condition that the mass is unity, which we adopt in the present paper. Then, the eigenvalue of Eq. (13) gives the second derivative of with respect to .
In this way, we obtain a series of states on the collective path, the collective potential of Eq. (9), and the collective inertial mass equal to unity by tuning the scale of . Thus, the collective Hamiltonian is constructed as
(16) 
The canonical quantization of this Hamiltonian immediately leads to
(17) 
ii.2 Mapping to different variables
In order to obtain a physical picture of the collective dynamics, it is often convenient to adopt an “intuitive” variable, such as the distance between two nuclei, . Of course, the optimal collective coordinate , determined by the ASCC solutions, is different from , in general. Nevertheless, as far as the onetoone correspondence between and is guaranteed, we may use the variable to modify the scale of the coordinate. Without affecting the collective dynamics, the collective Hamiltonian, Eqs. (16) and (17), is rewritten in terms of .
Let us denote a new variable as , defined by the expectation value of the corresponding onebody Hermitian operator . For instance, the operator of the relative distance between two symmetric nuclei () is given by
(18) 
where is the step function. We also assume the onetoone mapping between and , and its inverse function . The transformation of the collective potential, , is trivial: . In contrast, the inertial mass is transformed as
(19) 
where we use . The inertial mass is not constant but depends on . The collective Hamiltonian is rewritten as
(20) 
The quantization identical to Eq. (17) is given by the Pauli’s prescription Pauli (1933)
(21) 
The mass requires the calculation of the derivative, or , in Eq. (19). These quantities can be obtained by use of the local generator ,
(22)  
where are the ph matrix elements of with respect to the state and assumed to be real. Since this calculation can be performed using the local quantities at , it has an advantage over the conventional finite difference, , with two adjacent points, and , on the collective path. Thus, we use Eq. (22) for calculation of the derivatives.
In the present ASCC method, the variable is merely a parameter to represent the collective coordinate . It should be emphasized that this is different from assuming the collective coordinate as . First of all, the potential is different, . Here, is calculated by minimization of the total energy with a constraint on . Even if the state is close to , the inertial masses for the motion along the direction can be very different from . The ASCC method guarantees a blockdiagonal form of the inertial tensor between the collective coordinate and the rest of intrinsic degrees perpendicular to . In contrast, the inertial mass tensor for the coordinate is not blockdiagonal in general. Thus, we need to adopt its diagonal element which is different from :
(23) 
where is the inertial mass tensor for the intrinsic motion. Last but not least, the inertial mass is usually calculated according to the Inglis’s cranking formula, (See Sec. III.3.1). The cranking mass cannot take into account the effects of the timeodd mean fields. In contrast, the ASCC inertial mass , which is determined from the moving RPA equation (8), reflects the presence of the timeodd mean fields. Therefore, even if the collective coordinates and are identical, the calculated inertial masses may be different. For instance, for the translational (centerofmass) motion of the nucleus, the cranking mass fails to reproduce the total mass, , when the effective mass is different from the bare nucleon mass . It is compensated by the timeodd effect in the ASCC inertial mass, that leads to the exact relation, .
ii.3 Numerical algorithm and details
Coordinatespace and mixed representation
In this paper, we adopt the BKN energy density functional Bonche et al. (1976)
for the Hamiltonian
(24)  
where is the sum of the Yukawa and the Coulomb potentials,
(25) 
We take the same parameter set as in reference Bonche et al. (1976).
For the BKN energy density functional, it is convenient to utilize the coordinatespace representation. Each singleparticle wave function is represented in the 3D grid points of the square mesh, with , where is the mesh size. Although every quantity is defined locally at , in this subsection, we omit the collective coordinate for simplicity, such as . The 3D space is a rectangular box of volume fm with mesh size fm.
We adopt the mixed representation for the moving RPA equation: The particlestate indices , are replaced by the coordinate . Thus, the generator of Eq. (4) is represented as
(26)  
where and . Since the coordinate indices contain not only the particle states but also hole states, we should remove the hole parts. Using the projection operator, , this is done by replacing by
(27) 
where . Equivalently, is replaced by . Similar modification is performed for the generator of Eq. (3).
The matrices of Eq. (12) are represented as and the same for the matrix . The hole contributions are removed in the same manner as Eq. (27). For instance, in the ph representation becomes
(28)  
which can be discretized as
(29) 
where .
Although we remove the holehole contributions in this manner,
the RPA matrices are oversize and contain redundant components.
Therefore, the diagonalization of the moving RPA equation
produces spurious solutions that consist of only holehole elements
Finite amplitude method for the moving RPA solution
Solutions of the moving RPA equation (13) determine the local generators , then is obtained from Eq. (15). To evaluate the matrix elements of in Eq. (13), we adopt the finite amplitude method (FAM) Nakatsukasa et al. (2007); Avogadro and Nakatsukasa (2011, 2013); Stoitsov et al. (2011); Liang et al. (2013); Hinohara et al. (2013); Nikšić et al. (2013); Pei et al. (2014); Kortelainen et al. (2015), especially the matrix FAM (mFAM) prescription Avogadro and Nakatsukasa (2013). The FAM requires only the calculations of the singleparticle Hamiltonian constructed with independent bra and ket states Nakatsukasa et al. (2007), providing us an efficient tool to solve the RPA problem.
Let us assume that the state is determined from the moving HF equation (6). The singleparticle states and their energies of the hole states are defined by
(30) 
where is the singleparticle Hamiltonian reduced from with . The selfconsistent density is given by . According to Ref. Avogadro and Nakatsukasa (2013), the matrix elements, can be calculated as follows:
(31) 
Here, again, the dependence is omitted for simplicity. Using a small real parameter , the mFAM provides the elements by
(32) 
where is defined as
(33) 
Note that Eq. (32) requires only the operation of the singleparticle Hamiltonian on the hole orbits. In addition, the singleparticle Hamiltonian can be replaced by the HF singleparticle Hamiltonian in Eqs. (31) and (32). It is trivial, for Eq. (32), to see the term is canceled by the subtraction. For Eq. (31), this is because the hole components are always removed from and , and the generator has only ph and hp components.
Imaginarytime method for the moving HF solution
Let us assume that the is determined from the moving RPA equation, and now we want to move to the next point on the collective path (). This can be done by solving the moving HF equation (6), with the following constraint:
(34) 
which controls the step size of the collective coordinate . Equation (34) can be understood as
(35) 
by the use of the local generator and the canonicity condition (5).
Hereafter, in this subsection, the quantities without the explicit dependence mean those defined at , such as . The moving HF equation (6) is iteratively solved using the imaginarytime method Davies et al. (1980) which is efficient in the coordinatespace representation. Each singleparticle wave function is evolved as , where
(36) 
with a small real parameter . is the singleparticle Hamiltonian calculated with the density . Here we approximate in equation (6) by , provided that is small enough. The Lagrange multiplier is determined by the constraint (34). In the first order in , is given by
(37) 
at each iteration. Here, the traces are calculated as
(38)  
(39) 
with and . Note with . In actual calculations, we also have constraints on the center of mass and the direction of the principal axis. These additional constraints are easily taken into account by extending Eqs. (36) and (37).
According to Eq. (6), in principle, we should use in Eq. (36) instead of , namely the generator at the same point . The prescription given in Sec. II.3.3 actually approximates the generator by the one at the previous point . The approximation significantly reduces the computational task. This approximation turns out to be very good as far as is small enough. After moving the state from to with small , the selfconsistency can be checked in the following way. At we calculate the generator by solving the moving RPA equations. Replacing by in Eq. (36) and changing the constraint condition Eq. (34) to , the selfconsistency between and is guaranteed if the further imaginarytime evolution of Eq. (36) keeps the state invariant. This is confirmed for the present case. The validity is also confirmed by the fact that the final result is invariant with respect to change of the step size .
Summary of the numerical algorithm
We choose the HF ground state as the starting point of the collective path, . The HF ground state is always the solution of Eq. (6) with , at which the moving RPA equation becomes identical to the conventional RPA equation. Therefore, without knowing the generator , the starting point can be determined. The procedure to construct the collective path is given as follows:

Calculate the HF ground state, .

Solve the RPA equation to obtain and .

When , , and are provided, solve the moving HF equation to obtain the state , according to the method described in Sec. II.3.3.

Solve the moving RPA equation to obtain the generators, and , according to the method described in Sec. II.3.2.
For the step 4 above, we choose the inertial mass . Then, the weak canonicity condition (5) determines the scale of as
(40) 
The scale transformation from to is performed by changing the inertial mass according to Eqs. (19) and (22).
Algorithm for fully consistent solutions
Since Be is one of the simplest cases, we also try another method to get the fully selfconsistent solutions of the , and that simultaneously satisfy Eqs. (6), (7), and (8). For Be, the conventional constrained calculation on may produce approximate solutions, . Thus, we adopt as the initial trial wave functions, and start the following iteration procedure.

Use this to solve the moving HF equation (6) with the constraint .
We also use the initial trial states prepared by the CHF calculation with constraint on the relative distance . Although the initial states are different from those obtained with the operator, after the iteration of (i)(iii) converges, they reach the same selfconsistent solutions, and .
It should be noted again that the prescription in Sec. II.3.3 is significantly easier than the present iteration (i)(iii). At every point , the selfconsistency requires us to solve the moving RPA equations many times to determine the selfconsistent state . We have confirmed that the solution of these iteration procedures (i)(iii) is practically identical to the one obtained with the algorithm in Sec. II.3.3.
Iii Numerical results
In this work the BKN energy density functional is adopted as a test for the numerical application of the ASCC method. The BKN energy density functional is rather schematic, thus, we should take the following results in a qualitative sense.
iii.1 Results of the RPA calculation at the ground state
If the frequency is positive () for Eq. (13), we may construct the normalmode excitation operator from the generators as
(41) 
For a Hermitian onebody operator , defined by Eq. (4) with replacement of , the transition matrix element is given by
(42)  
We assume that, for the coordinate representation of the operator such as or , the hole components are always projected out according to Eq. (27). The collective character of the state can be identified by choosing the onebody operator . For instance, the translational motion along axis is identified by a sizable transition matrix element of the centerofmass operator, . For the relative motion of twoalpha particles, we may choose the mass quadrupole operator (Sec. III.1.2).
Translational motion of a single particle
First, we show results for the single particle. In this case, the model space is a sphere of radius fm with various mesh sizes fm. Note that the ground state of the system is a trivial solution of the ASCC equation (6). We can clearly identify the three translational modes for , , and directions, degenerated in energy at MeV. Using smaller mesh size, the eigenfrequency of the translational motion approaches to zero. There are no lowlying excited states in the particle because of its compact and doublyclosed characters. The calculated energy of the lowest excited state is larger than 20 MeV.
Using Eqs. (19) and (22) with as the center of mass, we calculate the inertial mass of the translational motion of the particle. Figure 1 shows the results calculated with different mesh size of the 3D grid. Since this is the trivial centerofmass motion of the total system, this should equal the total mass, with . As the mesh size decreases, the total mass certainly converges to the value of . In the following, we adopt the mesh size fm.
Relative motion of two particles in Be
Figure 2 shows the calculated eigenfrequencies for the ground state of Be and the two well separated ’s at distance fm. Since the ground state of Be is deformed, there appear the rotational modes of excitation as the zero modes, in addition to the three independent modes of the translational motion. Because of the axial symmetry of the ground state, the rotation about the symmetry axis ( axis) does not appear. In Fig. 2 the calculation produces two rotational modes of excitation around 2.8 MeV with large transition matrix element of the quadrupole operator, . The finite energy of these rotational modes comes from the finite mesh size discretizing the space.
Besides these five zero modes, the lowest mode of excitation turns out to have a sizable transition strength of the quadrupole operator . This mode corresponds to the elongation of Be. The transition density is given by
(43)  
The left panels of Fig. 3 show the density profile of Be and the transition density corresponding to the lowest RPA normal mode. We can see an elongated structure along the direction in the ground state. The lowest mode of excitation corresponds to the change of its elongation (vibration).
We also perform the same calculation for the state in which two particles are located far away, at the relative distance fm. In the right panel of Fig. 3, we clearly see that the two particles are well separated, and the quadrupole mode in fact corresponds to the translational motion of the particles in the opposite directions, namely, the relative motion of two ’s. The excitation energy almost vanishes for this normal mode (Fig. 2).
iii.2 Results of the ASCC method
In Sec. III.1.2, we show that the the lowest quadrupole mode of excitation at the ground state of Be may change its character and lead to the relative motion of two ’s at the asymptotic region. We adopt this mode as the generators of the collective variables , then, construct the collective path.
Collective path, potential, and inertial mass
We successfully derive the collective path connecting the ground state of Be into the wellseparated two particles. The inertial mass is taken as unity and the collective potential is calculated according to Eq. (9). Then, according to Sec. II.2, the collective coordinate is mapped onto the relative distance with Eq. (18). Figure 4 shows the obtained potential energy along the ASCC collective path. As a reference, we also show the pure Coulomb potential between two particles at distance , , where is the calculated ground state energy of the isolated particle. Apparently, it asymptotically approaches the pure Coulomb potential. As two ’s get closer, the potential starts to deviate from the Coulomb potential at fm and finally reaches the ground state of Be. The ground state is at fm, and the top of the Coulomb barrier is at fm. Note that the path is determined selfconsistently without any a priori assumption.
With this calculated potential, we may check the selfconsistency of the ASCC potential and the eigenfrequency. If the collective path perfectly follows the direction defined by the local generators at each point of , the second derivative of the potential should coincide with the eigenfrequency of the moving RPA equation. The almost perfect agreement between these is shown in Fig. 5.
For the region of fm, there exists some discrepancy between and . In this region, the Be nucleus has even more compact shapes than the ground state, then, the coordinate and become almost orthogonal to each other, losing the onetoone correspondence between them. In other words, the states change as gets smaller, but keep almost constant. In addition, the moving RPA frequency becomes larger than the particle threshold energy, entering in the continuum. Thus, in this region of fm, the results somewhat depend on the adopted box size.
Figure 6 shows the obtained inertial mass as a function of for the scattering between two ’s As the two ’s are far away, the ASCC inertial mass asymptotically produces the exact reduced mass of . This means that the collective coordinate becomes parallel to the relative distance , even though we do not assume so. At fm, the value of inertial mass increases. This is due to the decrease of the factor in Eq. (19). Making the system even more compact than the ground state, rises up drastically, which means that the coordinates and become almost orthogonal.
Phase shift for scattering
The ASCC calculation provides us the collective Hamiltonian along the optimal reaction path. Using this, we demonstrate the calculation of nuclear phase shift. We should take this result in a qualitative sense, because of a schematic nature of the BKN energy density functional.