Adiabatic hyperspherical approach to large-scale nuclear dynamics

We formulate a fully microscopic approach to large-scale nuclear dynamics using a hyperradius as a collective coordinate. An adiabatic potential is defined by taking account of all possible configurations at a fixed hyperradius, and its hyperradius dependence plays a key role in governing the global nuclear motion. In order to go to larger systems beyond few-body systems, we suggest basis functions of a microscopic multicluster model, propose a method for calculating matrix elements of an adiabatic Hamiltonian with use of Fourier transforms, and test its effectiveness.

We formulate a fully microscopic approach to large-scale nuclear dynamics using a hyperradius as a collective coordinate. An adiabatic potential is defined by taking account of all possible configurations at a fixed hyperradius, and its hyperradius dependence plays a key role in governing the global nuclear motion. In order to go to larger systems beyond few-body systems, we suggest basis functions of a microscopic multicluster model, propose a method for calculating matrix elements of an adiabatic Hamiltonian with use of Fourier transforms, and test its effectiveness.

I. INTRODUCTION
Atomic nuclei present a unique example of self-bound, finite quantum many-body systems. They not only exhibit a variety of excitation modes but also decay or fission into two or a few fragments. Exploring the excitation mechanism based on single-particle, collective and clustering degrees of freedom is an interesting subject. Intrinsically different shapes such as prolate-oblate may coexist or mix at close energies, leading to the so-called large-amplitude collective motion [1]. Spontaneous fission and sub-barrier fusion are also typical examples of the collective motion that involves a large-scale change of the nuclear size [2,3]. A fully microscopic description of their dynamics is still a long-standing challenging problem.
All of the above phenomena should in principle be described starting from a Hamiltonian of the system. What is often performed is, however, to solve an equation of motion with some constraints [4] or to calculate energy surfaces assuming different shapes in order to look for a path along which the collective motion proceeds. In the case of a deep sub-barrier fusion an initial fragment decomposition is maintained for the whole fusion process and the relevant fusion potential is calculated as a function of the relative distance of the fragments. It is hard for that approach to take into account couplings with configurations corresponding to a different mass distribution of the fragments. Since the phenomena are very complicated, those approaches sound reasonable. However, neither the geometrical shape nor the relative distance between the fragments is a nuclear collective coordinate in a strict sense. A question thus arises of whether or not we can describe the large-scale dynamics by employing a true collective coordinate.
The purpose of this paper is to make use of a hyperradius as a collective coordinate, and to step forward for a consistent formulation of the large-scale dynamics together with the underlying collective potential. Most of the needed ingredients are available in the literature. The hyperradius is a global coordinate that measures matter size, and it is widely used in three-body problems [5][6][7][8]. After pioneering work with the hyperspherical approach [9], its extension to N -body systems has been proposed for solving various problems [10][11][12][13][14]. A common foundation of all the hyperspherical approaches is that a total wave function of the system is expanded in terms of a product of the hyperradial and hyperangular functions. There are two types of realization for describing the hyperangular functions. One is to use hyperspherical harmonics [10,11], and the other called an adiabatic hyperspherical approach is to employ channel wave functions that are defined by diagonalizing the hyperangular part of the Hamiltonian [13,14]. The former has the advantage that the hyperspherical harmonics are well-known eigenfunctions of the angular part of the multi-dimensional Laplacian, but its use in a real problem is fairly complicated and so far limited to fewparticle systems. Moreover, a convergence with that expansion is rather slow. See the first paper in Ref. [11] for details on the development and difficulty. The latter is widely used in atomic and molecular physics. Since an adiabatic hyperspherical potential defined there reflects the large-scale change of the system, we adopt the adiabatic hyperspherical approach in what follows.
The equation of motion in the hyperspherical method is the same independent of the number of particles in the system, which is an appealing feature of the hyperspherical approach. In spite of the various efforts, only small systems have so far been investigated mainly because calculating the matrix element of an adiabatic Hamiltonian is still not trivial and solving its eigenvalue problem is hard for general N -body systems. Correlated Gauss functions (CG) are employed for studying cold atom physics and electron-positron systems in Refs. [13,14], but their application is limited to few-body systems with the total orbital angular momentum L = 0 and 1. The use of harmonic-oscillator shell-model wave functions is discussed in Ref. [11] for calculating the matrix element needed in the hyperspherical approach. The oscillator basis is convenient for representing such one-centered configurations that are not highly excited from the ground state, but it is not flexible enough to cope with a description of large-scale change such as, for example, the clustering and fragmentation. We instead need basis functions of cluster type to describe such configurations. We attempt here an extension of microscopic multicluster wave functions [15,16] used to describe the structure of light nuclei. In the multicluster model the intrinsic fragment wave functions are described with shell-model type configurations, while the relative motion among the fragments is described with the CG [17,18]. We apply a Fourier integral for evaluating the matrix element as it is applicable for any type of many-body basis function.
The structure of the present paper is as follows. We define in Sect. II the hyperspherical coordinates and separate the kinetic energy of the system into hyperradial and hyperangular parts. In Sect. III we define the eigenvalue problem of the adiabatic Hamiltonian and present the equation of motion for hyperradial functions. In Sect. IV we discuss qualitative features of the adiabatic potential together with a separation of active and inactive degrees of freedom. In Sect. V we define basis functions of the multicluster model and give a method for calculating the matrix element integrated over the hyperangles together with examples for the overlap and kinetic energy. In Sect. VI we show how to extract the evolution of intrinsic shapes of the system as a function of the hyperradius. In Sect. VII we touch on an eigenvalue problem of the full Hamiltonian with a constraint of the mean-square matter radius in comparison with the present approach. Conclusions are drawn in Sect. VIII.

II. HYPERSPHERICAL COORDINATES
We start from defining the hyperradius for a general case consisting of K particles. The mass of the ith particle is A i in units of a suitable mass m. By denoting its position coordinate by R i , we define a set of Jacobi coordinates by where A 12...i = i j=1 A j and µ i is the reduced mass factor, µ i = A 12...i A i+1 /A 12...i+1 . The square of the hyperradius ρ is defined by which is also rewritten in several ways as where R cm is the center-of-mass (cm) coordinate of the system, Note that mρ 2 is equal to the trace of the moment of inertia tensor of the system.
It is straightforward to extend the above definition to an N -nucleon system. Protons and neutrons are assumed to have an equal mass, the nucleon mass, which is taken as m. By denoting the nucleon's position coordinate by r i , we define Jacobi coordinates as with µ i = i/(i + 1). Then ρ 2 reads We often use a matrix notation. For example, x = (x i ) stands for an N − 1-dimensional column vector or an (N − 1) × 1 matrix, andx stands for its row vector. The ρ 2 is simply written as a scalar product, ρ 2 =xx. It is clear that ρ 2 is equally defined by any coordinates that are related to x by an orthogonal transformation. In fact ρ 2 is independent of any choice of such coordinates.
A measure of the nuclear size, ρ 2 /N is an operator for the mean-square matter radius. Symmetric with respect to the nucleons' coordinates, ρ is a collective coordinate that has a unit of length. The other 3N − 4 coordinates are hyperangle coordinates denoted by Ω collectively. The volume element for integration reads where d is the dimension of the spatial coordinates excluding the cm coordinate It is well known that the volume V d of a d-dimensional hypersphere with radius ρ = ρ 0 is given by The volume element in the single-particle coordinates reads dr 1 dr 2 . . . dr N = N 3/2 dxdR cm . Let us introduce dimensionless coordinates ξ i by x i = ρξ i . They are subject to the constraint An explicit form of Ω may be constructed from the N − 1 vectors ξ i , but it is not needed in what follows. It should be noted, however, that a variety of configurations or shapes of the nucleus correspond to different functions of Ω. The total kinetic energy T of the N -nucleon system, with its cm kinetic energy T cm being subtracted, is separated into hyperradial (T ρ ) and hyperangular (T Ω ) parts: The hyperangular kinetic energy T Ω may be expressed as where K(Ω) 2 is the square of the grand angular momentum. An explicit form of K(Ω) 2 is available in a recursive way together with the definition of Ω [6,13]. Suppose that the N -nucleon system develops into K fragments or clusters each of which has N i nucleons ( . It is convenient to divide ρ 2 of Eq. (5) into two groups: where R i is the cm coordinate of the ith fragment and ρ 2 i is its squared hyperradius, The first term ρ 2 in of Eq. (12) gives a measure of the sum of the squared matter radii of the fragments (each ρ 2 i /N i is the mean-square matter radius of the ith fragment), while the second term ρ 2 rel is exactly the same as that of Eq. (3) with A i = N i , giving a measure of the spatial extension of the relative motion of the fragments. It is natural to arrange the coordinates into cluster-internal and cluster-relative to describe the motion of the K fragments. The cluster-internal coordinates, denoted by (x 1 , x 2 , . . . , x N −K ), consist of a collection of Jacobi coordinates of each fragment, and the cluster-relative coordinates denoted by (x N −K+1 , x N −K+2 , . . . , x N −1 ) are Jacobi coordinates as defined by Eq. (1). Clearly ρ 2 is independent of the number of fragments into which the N -nucleon system develops.

III. EQUATION OF MOTION IN ADIABATIC HYPERSPHERICAL EXPANSION
To solve a Schrödinger equation for the system in the hyperspherical method, a total wave function Ψ is usually expanded in terms of a complete set of the hyperspherical harmonics or K-harmonics Y λ (Ω): Ψ = ρ −(d−1)/2 λ χ λ (ρ)Y λ (Ω) [10,11,19,20]. Here Y λ (Ω) is an eigenfunction of T Ω labeled by λ. The hyperradial functions χ λ (ρ) are determined from a set of coupled-channels equations. This method is successfully used in nuclear few-body systems [5]. However, the number of hyperspherical harmonics needed to reach a converged solution becomes very large at large ρ values [21]. Moreover, the coupling matrix elements between different Y λ (Ω) are of the same orders of magnitude as the diagonal matrix elements especially for the Coulomb interaction, which also makes the convergence slow.
In low-energy phenomena, the hyperradial motion is expected to be slow compared to the hyperangular motion. Thus an adiabatic potential that takes account of all possible hyperangular motion at a fixed hyperradius gives insight into the dynamics of the system's evolution [22]. We adopt the adiabatic hyperspherical expansion method [6,8,[12][13][14] used extensively in atomic and molecular physics.
We define an adiabatic Hamiltonian H ad by Here V is the total potential energy and H = T + V is the total Hamiltonian of the system. The nucleon-nucleon interaction of V is assumed to be an effective interaction that contains no strong short-ranged repulsion. Assuming that V contains no derivative of ρ, we solve an eigenvalue problem of the Hermitian operator H ad , to obtain channel wave functions Φ ν (ρ, Ω) and real adiabatic potentials U ν (ρ) that are labeled by ν. The quantum numbers of H such as spin-parity J π are preserved as those of H ad , and the antisymmetry requirement on Ψ applies on Φ ν (ρ, Ω) as well. Note that ρ appears parametrically in Eq. (15). At fixed ρ, all possible couplings among various hyperangular configurations are taken into account to obtain U ν (ρ) and Φ ν (ρ, Ω). The Φ ν (ρ, Ω) form a set of orthonormal functions at each ρ, where . . . Ω indicates that the integration is carried out over Ω with ρ being fixed. Actually, the Φ ν (ρ, Ω) also contain spin and isospin coordinates that have to be integrated, but they are omitted for the sake of simplicity. Apparently The Schrödinger equation, HΨ = EΨ, is solved by expanding Ψ in terms of Φ ν (ρ, Ω): The normalization of Ψ is ν ∞ 0 |f ν (ρ)| 2 dρ = 1 for a bound state. The hyperradial functions f ν (ρ) are determined by solving a set of coupled-channels equations, with non-adiabatic coupling terms Equations (15), (18), and (19) give a microscopic description of the large-scale dynamics. A unique advantage of the adiabatic hyperspherical approach is that both lower and upper bounds to the exact lowest energy of H are readily obtained [23,24]. As shown in Appendix A, we have Its differentiation with respect to ρ leads toQ whereQ is non-negative, and consequently Q νν (ρ) ≤ 0. The potential defined by always satisfies W ν (ρ) ≥ U ν (ρ). The lowest eigenvalue E obtained by truncating Eq. (18) to a single-channel equation with the lowest adiabatic potential U 0 (ρ) or W 0 (ρ) gives a lower or upper bound to the exact lowest energy of H. See Appendix A for details. Convergence of the solution of Eq. (18) is checked by increasing the number ν of channels. A time-dependent Schrödinger equation is convenient for studying how final configurations in e.g. few-body decay and sub-barrier fusion evolve from their initial states. The wave function at time t is assumed as Once f ν (ρ, 0) are given, f ν (ρ, t) for t > 0 are determined from the equation

IV. HYPERRADIUS DEPENDENCE OF ADIABATIC POTENTIAL
The ρ-dependence of U ν (ρ) or W ν (ρ) governs how the nucleus responds to its change of size. The kinetic energy and the Coulomb potential respectively give 1/ρ 2 and 1/ρ contributions to U ν (ρ) at large ρ values. Short-range pairwise nuclear interactions give a ρ −n (n ≥ 3) contribution [25]. Let us focus on the lowest adiabatic potential with the same spin-parity J π as that of the ground state. U 0 (ρ) has a minimum at ρ ≈ ρ min corresponding to the matter size of the ground state. As ρ decreases from ρ min , U 0 (ρ) rises because of a loss of nuclear potential energy as well as an increase in the kinetic energy. As ρ increases from ρ min , various configurations contribute to determining Φ 0 (ρ, Ω). Here, deformations, shell effects, couplings with different modes and so on participate in determining the adiabatic potential. U 0 (ρ) reaches a peak at some ρ value or may even have a couple of local peaks at different ρ values. As ρ increases further, U 0 (ρ) approaches the lowest decay threshold of the nucleus.
The above global feature of the adiabatic potential well corresponds to the decomposition of ρ 2 in conformity with a formation of fragments or clusters. As shown in Eq. (12), the different decomposition of the fragments can be treated on an equal footing in the hyperspherical approach, which makes it possible to assess what configurations play an important role in determining the adiabatic potential. If one instead calculates a sort of adiabatic potential or potential energy surface as a function of the relative distance between two fragments, there is no way to compare such potentials for different fragment decompositions because their relative distances have a different meaning.
What fragment decompositions or configurations contribute to the lowest adiabatic potential clearly depends on ρ. The expectation value of H is a major contribution to the adiabatic potential (see Eq. (14)). We rewrite H according to the fragment decomposition: where T i + V i is the intrinsic Hamiltonian of the ith fragment, T rel the kinetic energy of the relative motion among the fragments, and V rel denotes the potential energies acting between the nucleons belonging to the different fragments. V rel depends on both cluster-internal and cluster-relative coordinates, thus causing a coupling of the relative motion among the fragments with their intrinsic motion. When ρ rel is so large compared to ρ in that the nucleon-nucleon interactions of V rel can be neglected and only the leading term of the Coulomb potentials of V rel is retained, V rel reduces to where Z i e is the charge of the ith fragment and ρ rel and Ω rel are the hyperradius and hyperangles constructed from the cluster-relative coordinates (x N −K+1 , x N −K+2 , . . . , x N −1 ). With increasing ρ the intrinsic motion of each fragment is stabilized toward its own ground state, while the configurations responsible for the relative motion are decoupled from the intrinsic motion. Both the coupling and decoupling of various degrees of freedom are naturally taken into account in the hyperspherical approach.
FIG. 1: The 10 lowest adiabatic potential curves of the three-α system with J π = 0 + that are taken from Ref. [26]. The hyperradius denoted R here is defined by When there are several thresholds corresponding to different fragment decompositions, avoided crossings of the adiabatic potential energy curves may occur. As an example, we show the case of 12 C that is described with a cluster model of three α-particles [26]. The eigenvalue problem (15) for H ad is solved accurately, and an analysis of the adiabatic potentials clarifies how the contributions of the hyperangular kinetic energy, the nuclear potential and the Coulomb potential change as a function of ρ. Figure 1, taken from Fig. 2 of Ref. [26], displays the 10 lowest adiabatic potential curves for J π = 0 + . The lowest potential U 0 (R) has a minimum at R ≈ 3.5 fm, which is deep enough to support a bound state, that is, the ground state of 12 C. Furthermore, the lowest potential reaches a broad peak around 12 fm, corresponding to the second 0 + state of 12 C, the Hoyle resonance state. The adiabatic potential indicated by the solid line is dominated by the two-body 8 Be+α state and approaches the 8 Be+α threshold at large R, while the other potentials indicated by the dashed lines are all dominated by the three-α continuum states. As seen in Fig. 1 (b), an avoided crossing begins to occur at R ≈ 140 fm, which is because the three-α continuum state comes down closely to the two-body 8 Be+α state. Since the avoided crossing actually occurs within a small range of R, it may be hard to see it in the figure. Refer to Fig. 3 of Ref. [26] to confirm the crossing clearly. Since the 8 Be+α threshold is higher than the three-α threshold, a number of avoided crossings successively appear below the adiabatic potential indicated by the solid line. As is well known, the non-adiabatic coupling terms (19) may be singular especially when the avoided crossing is sharp, namely it occurs within a small range of ρ. In that case, a diabatic procedure is proposed for accurately solving Eq. (18) [6,27,28]. The slow variable discretization method combined with a complex absorbing potential makes it possible to solve Eq. (18) and to reproduce the energy and width of the Hoyle resonance in good agreement with experiment [26].
Let us speculate concerning the adiabatic potential curves of 252 98 Cf that are crucially important for determining its decay mode. The ground state of 252 Cf decays mostly by an α-particle emission. The rest is a spontaneous fission (SF), emitting 3.7 neutrons on average. To make things simple, we approximate the SF as occurring through a single channel of 140 54 Xe + 108 44 Ru + 4n. The two decay modes contain different numbers of fragments, two in α + 248 96 Cm and six in the SF, but the hyperspherical approach can treat both in a unified way. The threshold of α+ 248 Cm is 6.2 MeV below the ground state of 252 Cf, whereas that of the SF is 200.4 MeV lower than the ground state. See the schematic diagram of Fig. 2. The lowest adiabatic potential U 0 (ρ) approaches the SF threshold at large ρ. Above that threshold many U ν (ρ) curves, not drawn in Fig. 2, show up corresponding to the continuum states of the SF mode. A unique U ν (ρ) with the two-body α+ 248 Cm character appears high above the SF threshold. When moving inward from this asymptotic region, the Coulomb potential (27) produces a distinct difference between the two decay modes. The charge factor Z 1 Z 2 of the SF mode is more than ten times larger than that of the α channel. Thus those U ν (ρ) curves that are dominantly contributed by the SF configurations rise up rapidly, while the U ν (ρ) curve of the α channel increases much more slowly. At the avoided crossing point ρ ac , the lowest curve U 0 (ρ) comes very close to that of the α curve, and for ρ < ρ ac the α channel makes a dominant contribution to U 0 (ρ). With further decrease of ρ many different configurations begin to mix due to an increasing role of the nuclear interaction V . The U 0 (ρ) reaches a barrier top around some point and reaches its minimum at ρ min corresponding to the matter radius of the ground state of 252 Cf. Though much more complicated than the 12 C case, the gross feature of the adiabatic potential curves of 252 98 Cf should have some similarity to those of 12 C, and the decay branch of 252 98 Cf will be determined by solving Eq. (18).

V. MULTICLUSTER APPROXIMATION AND INTEGRATION OVER HYPERANGLES
Solving Eq. (15) is of vital importance in the adiabatic hyperspherical approach. Its accurate solution is obviously very hard except for few-body system. The difficulty is enhanced by the fact that the matrix element has to be calculated by integrating over Ω only. Some efforts have been made for extending to larger systems [11,13,14]. We take up this problem assuming the use of many-body wave functions that contain all the coordinates.
Before discussing the eigenvalue problem (15), we note that a usual approach defines an adiabatic potential barrier or energy surface at a given 'collective' coordinate by searching for a minimum of V for various parameters that characterize the nuclear density or shape [29]. This makes sense in that V is a major part of H ad , and because, since V is a function of ρ and Ω, its minimum gives information on the most important Ω values contributing to the lowest adiabatic potential. As mentioned before, the adiabatic hyperspherical approach can go beyond that by taking account of various couplings with different degrees of freedom.
Let us assume that the channel wave function Φ ν (ρ, Ω) at a given ρ is expanded in terms of suitable basis functions φ i (x): Equation (15) is then reduced to the following generalized eigenvalue equation for determining the coefficients C νi (ρ) and the adiabatic potential U ν (ρ): where H ij (ρ) and B ij (ρ) are adiabatic Hamiltonian and overlap matrices defined by We include only those basis functions that give a c-number ρ 2 for the expectation value of the squared hyperradius operatorxx: We face two problems. One is what basis functions we use for φ i (x). The other is how to calculate the matrix element in Eq. (30). The first one is crucially important for assessing the quality of Φ ν (ρ, Ω) and U ν (ρ). Though it is difficult to give a general answer, our ansatz is to employ a microscopic multicluster approximation [15,16]. This is because, as mentioned in Sects. I and IV, the structure change we are interested in includes a variety of configurations ranging from one-centered shell-model wave functions to those with a few fragments or subsystems. A general form of the multicluster wave function containing K fragments reads where A is an antisymmetrizer, Ψ i an antisymmetrized intrinsic state of the ith fragment containing N i nucleons and χ is the relative motion function for the fragments. The cluster-internal coordinates (x 1 , x 2 , . . . , x N −K ) are abbreviated as (z 1 , z 2 , . . . , z K ), where e.g. z 1 stands for the first N 1 − 1 Jacobi coordinates (x 1 , x 2 , . . . , x N1−1 ). The spin-isospin coordinates are again suppressed. In general Ψ i may represent not only the ground state of the fragment but also its excited state. The quantum numbers for characterizing Ψ i are omitted. The coupling of the angular momenta of the Ψ i s and χ to a total angular momentum JM is implicitly understood in Eq. (32). We presume φ i (x) to belong to the space spanned by Note that any states in {φ (K) (x)} are in general nonorthogonal to each other even when they belong to different K subspaces. The questions of what intrinsic states of the fragments are important and what K subspaces have to be included depend on a given system and energy range of interest. To proceed further, we assume that Ψ i is approximated by harmonic-oscillator shell-model wave functions, while χ is described well with a superposition of Gauss functions [17,18,30] as developed in few-body problems.
We have to calculate a matrix element for some operator O(x), by integrating over Ω at fixed ρ, say ρ 0 . The calculation of the matrix element of T ρ in H ad can be aided with use of the identity See Appendix B for an example. In calculating the matrix element of H, the cluster-intrinsic term (26)) may be replaced by using the observed energy E i of Ψ i . This approximation looks reasonable and practically useful because any nuclear interaction can not satisfactorily reproduce the saturation property of nuclear binding energies despite the fact that reproducing the threshold energy for the fragment decomposition is important in the present approach.
The second problem has so far been examined using integral transform techniques [14,31]. We use a δ function technique as in Ref. [14]. Using the expression for Dirac δ function we can express O(ρ 0 ) as a Fourier transform of F ρ0 (ω): Note that x i is changed to ρ 0 ξ i with a dimensionless variable ξ i . In Eq. (39) dξ stands for dξ 1 dξ 2 . . . dξ N −1 , where the integration range of each ξ i covers the whole three-dimensional space. Since e −iωξξ = N −1 k=1 e −iωξ 2 k results in a simple modification of the basis function, F ρ0 (ω) can be calculated with a technique developed in microscopic cluster models [32,33].
In some cases the Fourier integral (38) can easily be obtained by Cauchy's integral formula that reduces to a residue calculation. Whether or not we have a practical means for evaluating Eq. (34) for a general case depends on how fast and accurately the Fourier integral is computed. For this aim we test the Whittaker cardinal series or the Whittaker-Shannon interpolation formula [34]: where sinc x is the sinc function, sin x/x, and ω n = nh (n = 0, ±1, ±2, . . .) is the grid of the sampling points. The series (40) is known to converge if F ρ0 (ω) is a band-limited function. Because sinc nπ = δ n,0 , the series is exact at all the sampling points. It is in fact an expansion in terms of orthogonal functions {sinc π h (ω − ω n )} that have the properties: The third equation called the Dirichlet integral leads to an approximation for O(ρ 0 ): which is nothing but a trapezoidal rule for the integration. This result is due to the fact that the Fourier transform of the sinc function is the rectangular function and vice versa. To determine M , we need to know how fast F ρ0 (ω) decreases as a function of ω. The mesh size h (h < π) is determined by examining how accurate the expansion is at, e.g. ω = (n + 1 2 )h, the midpoint of ω n and ω n+1 . Other interpolations, e.g. a cubic spline interpolation may also be worthwhile testing because it leads to a simple expression for Eq. (38) and in addition the mesh size can be taken as piecewise variable. Once dF ρ0 (ω)/dω values at both boundaries of the interpolation are calculated, we can completely fix the interpolating function of the cubic spline.
Since ω-dependence of F ρ0 (ω) is of practical importance, we examine it for the diagonal matrix elements (φ i (x) = φ j (x)) of O(x) = 1 and T Ω in a very schematic model. As the model, we employ CG ignoring the antisymmetry requirement of the wave function and focus only on its spatial part. See Appendix B for some basic matrix elements with the CG. For a spherical CG, exp(− 1 2x Ax), the positive-definite symmetric matrix A is set to TrA −1 = 2 3 ρ 2 0 because of Eqs. (31) and (B2). We may choose A to be diagonal, A = (a i δ ij ), as far as the diagonal matrix element of O(ρ 0 ) is concerned. Our first choice for A is a uniform nuclear expansion, a i = a, leading to a hyperscalar Gaussian, exp(− 1 2x Ax) = exp(− 1 2 aρ 2 ). This function is totally symmetric and Ω-independent. By taking a as (N − 1)/a = 2 3 ρ 2 0 , the overlap matrix element is (see Eq. (B3)) Clearly |F ρ0 (ω)| becomes very small if ω is significantly larger than d/2. The Fourier transform (38) can be rigorously computed in this case. If d/2 is an integer, F ρ0 (ω) has a pole of order d/2 at ω = id/2, so that the integral is reduced to a residue calculation, yielding Even when d/2 is a half integer, we can derive the above result as follows. By the change of the integration variable, By changing the integration path to the Hankel contour and using Hankel's integral representation and Euler's reflection formula for the gamma function, we find the above integral to be 2π/iΓ(d/2). The result (44) is in fact trivial thanks to Eq. (8) if we note that the hyperscalar Gaussian at ρ = ρ 0 is exp(− 1 2 aρ 2 0 ) = e −d/4 and hence O(ρ 0 ) must be e −d/2 dΩ. We note that O(ρ 0 ) for O(x) = T Ω vanishes because the hyperscalar Gaussian is Ω-independent and, when acted on by T Ω , vanishes. This is also confirmed by using Eq. (B6). The next example is a 'symmetric fission', that is, the nucleus fissions into two identical fragments with mass number N/2 and only the relative motion between them expands with increasing ρ 0 . Let R rms (N ) denote the rootmean-square radius of a nucleus with mass number N , and set it equal to 3/5r 0 N 1/3 (r 0 = 1.1 fm). The matrix A for the symmetric fission is chosen as a 1 = a 2 = . . . = a N −2 = a, and a and a N −1 are determined by the condition where ρ f is fixed to N/2R rms (N/2). The mass number N is changed to 4, 40, and 240, and for each N ρ 0 is taken as ρ 0 = η √ N R rms (N ) (η = 1, 3, 5). Figure 3 displays log 10 r(ω) for the overlap, where r(ω) = |F ρ0 (ω)/F ρ0 (0)|. In the case of N = 4, the fall-off of r(ω) is slow with increasing ω and η. For example, log 10 r(ω) at ω = 5000 is −13.7, −11.4, −10.1 for η = 1, 3, 5, respectively. For N = 40, r(ω) rapidly drops to 10 −15 as a function of ω, but its decrease becomes slower for η = 5. This behavior is also valid for N = 240, and the decrease in ω becomes even slower with increasing η. As shown in Fig. 4, the ratio log 10 r(ω) for T Ω is very similar to that of the overlap.

VI. EVOLUTION OF INTRINSIC SHAPES
It is interesting to know how an intrinsic shape of the nucleus changes with increasing ρ. When a decay or an SF is considered as a tunneling through a barrier, the shape will give insight into where the fragments are formed and how they evolve during the passing through the barrier. The barrier is conventionally calculated by assuming some density distribution constrained with shape or deformation parameters such as quadrupole and octupole [2,29]. Such deformations are not observable, however. Our view is to reverse this approach. Since the nucleus should in principle preserve the total angular momentum, it is not trivial to imagine the intrinsic shape in the space-fixed frame. For example, any state with L = 0 is spherical in that frame, but it can happen that such state is intrinsically deformed and rotates. As shown in Ref. [35], the intrinsic two-α structure of the rotational state of 8 Be emerges from the wave function obtained by a quantum Monte Carlo calculation.
Following the procedure of Ref. [35], we can get the intrinsic density or deformation indicated by e.g. the lowest channel wave function Φ 0 (ρ, Ω). Since Φ 0 (ρ, Ω) is normalized as in Eq. (16), its square, P ρ (Ω) = |Φ 0 (ρ, Ω)| 2 , gives the probability density as a function of Ω at a given ρ. First, we generate many sampling points Ω 1 , Ω 2 , . . . , Ω M according to the distribution of P ρ (Ω) using the Metropolis-Hastings algorithm. Secondly, we define a body-fixed intrinsic frame for each Ω j = (ξ j 1 , ξ j 2 , . . . , ξ j N −1 ) as follows. By using Eq. (4) together with R cm = 0, Jacobi coordinates x j i = ρξ j i (i = 1, 2, . . . , N − 1) specify the positions of N nucleons (r j 1 , r j 2 , . . . , r j N ) in the space-fixed frame. From these position vectors, we calculate the moment of inertia tensor where r j i α (α = x, y, z) is the Cartesian component of r j i . Diagonalizing the 3 × 3 symmetric matrix I j determines the principal moments of inertia, which define the axes of the intrinsic frame. For example, the axis is called x ′ , y ′ , z ′ in increasing order of the principal moment of inertia. The direction of the axis also has to be chosen consistently. By reading (r j 1 , r j 2 , . . . , r j N ) as (r ′ j 1 , r ′ j 2 , . . . , r ′ j N ) in reference to the intrinsic frame, we obtain the desired position coordinates of N nucleons in the intrinsic frame. Finally, accumulating these position coordinates over j = 1, 2, . . . , M leads to the intrinsic single-particle density at ρ. Once the intrinsic density is obtained, it is easy to extract multipole deformations.

VII. EIGENVALUE PROBLEM OF HAMILTONIAN WITH RADIUS CONSTRAINT
It looks as though the adiabatic hyperspherical approach has some relationship to an eigenvalue problem of the Hamiltonian with a constraint [4]. Let us attempt to find a solution of the Schrödinger equation by first constraining the expectation value of the squared hyperradius to a fixed ρ 2 value. Suppose that the solution is expanded in terms of some basis functions: The constraint (31) is not necessarily imposed on φ i (x) itself, but we demand the solution we are looking for to satisfy the condition Here κ is a label related to the Lagrange multiplier. The unknown coefficients G κi (ρ) and the energy eigenvalue E κ (ρ) are determined from the following equation where H, B, and Q are matrices defined by Unlike Eq. (30), the above matrices are obtained by integrating over the whole coordinates. To determine the coefficients G κj (ρ) from Eq. (50), the value of κ has to be given. Actually κ should be such that both Eqs. (49) and (50) are simultaneously met. Apparently Ψ ρ ′ κ ′ (x) and Ψ ρκ (x) are not orthogonal to each other even for ρ ′ = ρ. The next step is to use the generator coordinate method in which a solution Ψ for the Schrödinger equation is assumed as The coefficients C κ (ρ) are determined from the Hill-Wheeler equation which should be satisfied for any ρ ′ and κ ′ values. An approximate solution to the Hill-Wheeler equation gives an upper bound to the ground-state energy. Note that the adiabatic hyperspherical approach gives both lower and upper bounds as discussed in Sect. III.
We refer to two interesting calculations with a constraint in comparison to the adiabatic hyperspherical approach. One is a Hartree-Fock-Bogoliubov calculation performed by constraining the mean-square radius, N i=1 r 2 i /N , to study how self-conjugate nuclei fragment into α clusters [36]. As Eq. (5) indicates, this constraint is equivalent to that of ρ 2 provided the contribution of R 2 cm to the squared radius remains a constant. The treatment of the cm motion in Ref. [36] does not satisfy this condition as usual in a mean-filed model. It would be a challenge for the mean-field approximation to cope with such diverse structure at large distances that is composed of different numbers of fragments. What should be further pursued at this moment is to establish the essential relationship between the adiabatic hyperspherical approach and 'beyond mean-field' calculations or configuration interaction calculations that constrain the mean-square matter radius.
Another is a simultaneous study of both α+ 6 He reactions and the structure change of 10 Be in a microscopic α+α+n+n model [37], in which a distance parameter between the two α-clusters is constrained. Since the motion of the two neutrons is restricted to either molecular or atomic orbits around the α-clusters, the main configurations included are α+ 6 He and 5 He+ 5 He two-body types. The adiabatic energy surfaces are calculated within that approximation. An avoided crossing is treated by the generator coordinate method. As noted in Sect. IV, the relative distance of the fragments is not a collective coordinate. If one constrains ρ 2 as the generator coordinate, it would be possible in the same four-body model to take account of possible couplings with the 9 Be+n channel that is the lowest threshold of 10 Be as well as the three-and four-body channels, 8 Be+n + n and α + α + n + n, that are open in the energy region treated in Ref. [37].

VIII. CONCLUSION
Stressing that the hyperradius is a collective coordinate, we have formulated a fully microscopic adiabatic hyperspherical approach to large-scale nuclear dynamics. The equation of motion for hyperradial functions is universal, independent of the number of nucleons, and enables one to consistently treat the dynamics from confined nuclear motion to relative motion among fragments in their asymptotic region. It is possible to describe in a unified way cases where the nucleus fragments into several channels. No spurious center-of-mass motion appears and couplings with different degrees of freedom can naturally be taken into account. These properties are due to the fact that both the squared hyperradius and the kinetic energy are flexibly decomposed into cluster-internal and cluster-relative quantities responding to the fragment formation.
The adiabatic potential as a function of the hyperradius plays a key role in the present approach. It is unambiguously defined solely by the Hamiltonian of the system, and there is no need to assume specific geometrical shapes or deformations to compute it. Conversely the shape or intrinsic density, if necessary, comes out after the adiabatic potential is obtained or the equation of motion for the hyperradial functions is solved. The calculation of the adiabatic potential involves the integration over all the coordinates but the hyperradius. Expecting that a microscopic multicluster model is a promising candidate for applying the present approach to larger systems, we have discussed the use of Fourier transforms for evaluating the matrix elements needed to obtain the adiabatic potential. The matrix elements can be obtained in exactly the same way as the usual matrix elements needed in nuclear many-body calculations. A merit of the Fourier transform technique is its simplicity, and test calculations indicate that accurate evaluations of the matrix elements are feasible.
Although the calculation of the adiabatic potential still requires much computer time for large systems, a real challenge is whether we can provide large enough basis functions to cover important configurations for fixed ρ. Further developments are certainly indispensable for a microscopic, realistic description of large-scale nuclear dynamics.
The hyperradial function f (ρ) has to vanish at ρ = 0. The ground-state energy reads where From the normalization condition of Φ(ρ, Ω), we obtain Thus P (ρ) must be pure imaginary or zero. If P (ρ) is not zero but pure imaginary, f * (ρ)df (ρ)/dρ in Eq. (A3) must also be pure imaginary because E exact is real. With f (ρ) = g(ρ) + ih(ρ), where g(ρ) and h(ρ) are real functions, f * (ρ)df (ρ)/dρ reads which leads to d{g(ρ) 2 + h(ρ) 2 }/dρ = 0. Thus g(ρ) 2 + h(ρ) 2 is a constant, and it must be zero because of f (0) = 0. Namely, f (ρ) vanishes identically, which can not be accepted. Using P (ρ) = 0 in Eq. (A3) leads to Suppose that for Φ(ρ, Ω) we take the Φ 0 (ρ, Ω) that gives the lowest adiabatic potential. The corresponding quantities Q(ρ) and U (ρ) in Eq. (A4) are denoted by Q 0 (ρ) and U 0 (ρ), respectively. It follows from the Ritz variational principle that If f (ρ) is chosen to be the solution of the equation (the adiabatic approximation), with the lowest eigenvalue E U , E U turns out to be an upper bound of E exact : E exact ≤ E U . Differentiating P (ρ) = 0 with respect to ρ leads to Since the last term in the square brackets is non-negative, we obtain By using the inequality U (ρ) ≥ U 0 (ρ) and choosing f (ρ) to be the solution of the equation (the Born-Oppenheimer approximation), with the lowest eigenvalue E L , we obtain a lower bound of E exact as E exact ≥ E L . If we calculate the expectation value of H for the wave function Ψ = ρ −(d−1)/2 f (ρ)Φ ν (ρ, Ω) with the νth channel wave function, we confirm Eq. (20) using the same argument as above.

Appendix B: Matrix elements with correlated Gaussians
In this appendix we calculate F ρ0 (ω), Eq. (39), using as φ i (x) the generating function g(s; A, x) [17,18,38,39] of the CG: where A is an (N − 1) × (N − 1) symmetric, positive-definite matrix and s = (s i ) is an (N − 1)-dimensional column vector to describe motion with non-zero orbital angular momentum. They are both parameters that characterize the CG. The constraint (31) reads Note that for the special case that A is diagonal, A = (a i δ i,j ), g(s; A, x) reduces to a product of Gaussian wave packets: g(s; A, x) = We present formulas for F ρ0 (ω) calculated between g(s; A, x) and g(s ′ ; A ′ , x). See Ref. [18] for details. The case with s = s ′ = 0 is given in Ref. [40].

Overlap
The function F ρ0 (ω) for O(x) = 1 is given by Here I is the (N − 1) × (N − 1) identity matrix. Since A + A ′ can be diagonalized by an orthogonal matrix, the matrix B can be diagonalized as well.
This integral can be performed analytically. In the case of s = s ′ = 0, we obtain (B10) Potential energy The matrix element for O(x) = V is conveniently calculated by expressing the distance vector of two nucleons as a combination of Jacobi coordinates where ζ is an (N − 1)-dimensional column vector determined by i and j. For a Gauss potential, O(x) = e −τ 2 (ri−rj ) 2 = e −τ 2x ζζx , F ρ0 (ω) reduces to that of the overlap. For s = s ′ = 0, we obtain F ρ0 (ω) = (2π) N −1 det(B + 2τ 2 ρ 2 0 ζζ) In the last step Sherman-Morrison formula, det(B + cζζ) = (1 + cζB −1 ζ) detB, is used, where c is a constant. By comparing this result with Eq. (B3) and by noting that B −1 is constrained by the condition (B2), we expect that the contribution of the potential energy V to the adiabatic potential behaves as (ρ/τ ) −3 for large ρ values. The ρ −3 dependence was found for a three-body system [25], but our result suggests that it is valid for many-body systems as well.
As an important application of Eq. (B12), we calculate the matrix element for the Coulomb potential, O(x) = 1/|r i − r j |. Using and the integral we obtain the matrix element for the Coulomb potential as As expected, the inverse ρ-dependence appears naturally.