GCM Solver ( Ver . 3 . 0 ) : A Mathematica Notebook for Diagonalization of the Geometric Collective Model ( Bohr Hamiltonian ) with Generalized Gneuss – Greiner Potential

The program diagonalizes the Geometric Collective Model (Bohr Hamiltonian) with generalized Gneuss–Greiner potential with terms up to the sixth power in β. In nuclear physics, the Bohr–Mottelson model with later extensions into the rotovibrational Collective model is an important theoretical tool with predictive power and it represents a fundamental step in the education of a nuclear physicist. Nuclear spectroscopists might find it useful for fitting experimental data, reproducing spectra, EM transitions and moments and trying theoretical predictions, while students might find it useful for learning about connections between the nuclear shape and its quantum origin. Matrix elements for the kinetic energy operator and for scalar invariants as β2 and β3 cos (3γ) have been calculated in a truncated five-dimensional harmonic oscillator basis with a different program, checked with three different methods and stored in a matrix library for the lowest values of angular momentum. These matrices are called by the program that uses them to write generalized Hamiltonians as linear combinations of certain simple operators. Energy levels and eigenfunctions are obtained as outputs of the diagonalization of these Hamiltonian operators.


Introduction and Motivation
Utilization of several important physical models is limited to specialists in the field and broad-scope programs are not generally available to a larger pool of scientists.One example of this is the Bohr-Mottelson model and its extensions and more modern variants.The original Bohr Hamiltonian [1,2] (employing a 5D harmonic oscillator) has been extended to a more complete model and it is known with a few different names: often called simply the Bohr-Mottelson model from the name of the two Nobel prize winning proposers, or it is called Geometric Collective Model (GCM) with reference to the fact that macroscopic (i.e., collective) variables are used to parameterize the geometric shape of a drop of quantum fluid and to describe the normal modes of vibrations of its surface, or it is called Roto-Vibrational model with reference to the ability of the model to describe on the same footing the quantum vibrational and rotational spectra that are observed at low excitation energy of even-even nuclei.These names merge with the so-called Generalized Collective Model of Gneuss and Greiner (Frankfurt model) [3,4] that used the simplest possible potential energy form built out of invariant operators (in turn obtained from proper tensorial couplings of the basic quantized canonical coordinates) that has anyway the property of giving all the relevant configurations: spherical, axially deformed (prolate or oblate), triaxial and γ-unstable.All these naming schemes usually refer to the model limited to the quadrupole degree of freedom, which is by far the most important, except in very high mass nuclei, but other types of motion (octupole, hexadecupole, etc.) have also been studied with generalized versions of this model.This model has the virtue of catching the essence of low-energy nuclear excitations and provides a convenient language for classifying important spectroscopic data in terms of β-and γ-bands and a convenient mathematical formalism for further calculations (quadrupole moments, electromagnetic transitions, etc.).It also has a number of underlying hypotheses, the most important of which is the starting classical continuous hydrodynamical description of oscillations of a liquid drop (that goes back to Rayleigh).This is a rather severe approximation when one considers the corpuscular nature of the nucleons that compose the nucleus, and implies some limitations, the most striking being the absence, in the model, of single-particle types of excitations and the connection with microscopic theories.Despite these limitations, the model can be profitably used to describe a subset of all possible microscopic configurations, in which all nucleons participate collectively to the motion, that is always found at low excitation energies in even-even nuclei.
The Bohr Hamiltonian and the model in all its variants and extensions has prominently returned to the cutting-edge of nuclear physics since the beginning of this century because new exact solutions of the Bohr Hamiltonian have been proposed by Iachello [5,6] as examples of critical point symmetries.These symmetries have spurred a whole new subfield in nuclear spectroscopy, both experimentally and theoretically, stimulating new experiments and new analytic or semi-analytic solutions [7,8] and formalisms to study the shape phase transitions, the corresponding critical points, the occurrence of shape coexistence and triaxiality in nuclei.Although it is clear that other models, such as the Interacting Boson Model, better address certain aspects, for example, the role of the finite number of particles that contribute to the collective motion, the Geometric Collective Model remains as a reference model for nuclear physics.Due to the complexity of the problem (five dimensions, nontrivial group theoretical constraints, and computational issues regarding basis truncation), the task of solving this model in a concise, applicable and practical way is not an easy one.It seems therefore desirable to "fill the gap" by proposing a simple, user-friendly, instructive computer program that gives a solution to the model.The only other published program that we have found is a relatively old Fortran 77 program by Trolteiner et al. [9] and the very interesting approach by Caprio [10].Of course over the years many researchers have developed their own programs for solving the same problem; the reason to publish ours is not because it is faster or better, but just because it is useful, instructive and reliable.There is also a recent Maple code [11], of which we became aware only at the time of revision of this manuscript, that uses an algebraic construction based on the Algebraic Collective Model.An alternative approach is also described in Ref. [12].
We use typeset fonts throughout to refer to program variables, functions and files, and italics to indicate mathematical symbols.
We wrote the code as a collection of Mathematica files in notebook form [13], because this is a very flexible and easy-to-use programming language with built-in graphical routines that makes it excellent as a teaching/learning tool.This work is an evolution of the Master Thesis work of one of the authors [14].
The paper is organized as follows: after this introduction, we remind basic facts and essential review papers about the Collective model in Section 2. In Section 3, we explain the workings of the code and, in Section 4, we describe the library of matrix elements that are used by the code for diagonalization.Section 5 contains examples and explanations on how to change to code to adapt it to other cases.Section 6 deals with the "Playing area", a portion of the code where one can add or change commands and make changes without touching the core part of the code.Finally some installation tips are given and conclusions are drawn.

Physics Background
Innumerable descriptions of the model are found in nuclear physics textbooks [15], where the focus is on physical ideas rather than on mathematical and group theoretical intricacies or computational techniques.For a good student-level introduction, see, for example, Refs.[4,16,17].Several books [18][19][20] and review articles [7,21] are also available that describe the model in greater detail.An introduction to the recent developments of the model and its importance in the study of quantum shape phase transitions can be found in Refs.[22,23].Let us remind here a few basic facts that are directly needed to understand our code.
The expression for the radius of the nuclear surface in spherical coordinates is expanded in terms of spherical harmonics with certain coefficients, α λ,µ that are called collective coordinates, as: where R 0 is the radius of a sphere, λ is the angular momentum quantum number and µ is the third component.The parameters α behave as a spherical tensor of rank λ under rotation, namely where µν (θ i ) are Wigner's rotation functions.The multipolar expansion in Equation ( 1) is truncated in the model at the quadrupole level, corresponding to λ = 2, as follows with five components µ = {−2, • • • , 2}.Then, the condition α * µ = (−1) µ α −µ is used that allows restricting the number of variables from five complex to five real variables.When passing from the laboratory frame to the intrinsic (principal axes) frame, five new collective variables are defined, a 2µ , that can be in turn transformed to other three Euler angles, θ i with i = 1, 2, 3, that just describe the orientation of the principal axes with respect to the laboratory axes and two polar variables, β and γ.The former is a radial variable (positively defined, 0 ≤ β ≤ ∞) in a polar plane where the angle is the latter, 0 ≤ γ ≤ 2π.It is not necessary to consider the whole plane because certain symmetry properties [1,2,4,7,21] allow restricting the domain to the 0 ≤ γ ≤ π/3 wedge.The remaining five wedges amount to a relabeling of the axes' names.
One can write a classical Hamiltonian and then, through a quantization procedure [4], obtain the following Hamiltonian operator that is called Bohr Hamiltonian.Here, B 2 is a mass parameter, L are the components of the angular momentum in the intrinsic frame and I i are the irrotational moments of inertia with respect to the three principal axes, given by with k = 0, 1, 2. The original model of Bohr and Mottelson has a γ-independent harmonic oscillator potential in β.With such a potential, the Schrödinger equation associated to the Hamiltonian in Equation ( 4) can be separated into an equation for the β variable and one in γ and Euler angles.Cases such as this, in which the potential does not depend on the angular variable, are called γ-unstable.More in general one might want to build a potential, function of both β and γ, that is sufficiently general to be able to encompass all possible symmetry limits that generate different nuclear shapes.The five α 2µ (and also a 2µ ) form a rank-2 tensor, as well as their conjugate momentum components π 2µ .One can use these coordinates to form two fundamental scalars: [α × α] and [[α × α] 2 × α] that correspond to β 2 and β 3 cos (3γ).All higher order scalars can be built out of these two operators.In the GCM, up to the sixth order invariants in β are used to write a general collective potential.It is necessary to reach at least this order because, otherwise, it would be impossible to get a triaxial minimum (it is the presence of cos (3γ) 2 that allows it).The most general potential with these properties is therefore: where the six coefficients

Matrix Elements
The 5D harmonic oscillator potential allows for separation of the Schrödinger equation into a radial and a gamma-angular part.The eigenfunctions are correspondingly factorized.The radial part is of the form [25][26][27]: where L b a are associated Laguerre polynomials and 1 F 1 is the confluent hypergeometric function.The normalized gamma-angular eigenfunctions can be written as where N is the normalization constant, D are Wigner's functions and g(γ) are the functions introduced by Bés [28,29] that we calculate with the method due to Corrigan and collaborators [27]: where all the indexes and coefficients use the notation of Ref. [27].These wavefunctions are used to calculate matrix elements of the various potential terms, taking fully into account all the restrictions coming from the group theoretical classification of the subalgebras of U(5) [25,26,30,31].

Program Explanation
The program automatically finds the path to the directory in which it is stored; alternatively, the user should set up the path to locally stored files and should set the maximum number of phonons that is used for diagonalization.Some supplementary modules are read in the first cell from the file basismodules.nb,namely the basis[Nph] and subbasis[J,Nph].The first generates the full basis for a given number of phonons and the second extracts the subbasis corresponding to a given value of angular momentum.Notice that the output of these modules is not directly written into the code.A number of arrays and variables are used: for instance, dtot gives the total dimension of the full basis, while Nst and Nst2, which are used throughout the program, give the subbasis' dimension and the dimension upon trimming the last six phonon numbers.
The code reads a library of matrices, called BasMat, one type for each potential energy term and angular momentum.The matrices are always read for N max = 40.This is a practical solution, because matrices corresponding to bases of smaller dimensions, if needed, can be extracted by cutting the appropriate upper-left submatrix.The matrices are used to write the hamiltonian matrix, as a linear combination with certain real coefficients (although physically unjustified, nothing prevents us from using complex coefficients) that must be supplied by the user.The paper of Habs et al. [32] gives a list of fitted coefficientsfor a few medium mass nuclei.One can use some of these examples to run the code, keeping in mind that some of those coefficients have been fitted on old spectroscopic data and that modern experiments have discovered new states and changed some of the older assessments.Therefore, with those coefficients, one can, at most, hope to reproduce their results, while to get a better result one should select the states of interest, run again some multidimensional fitting procedure (least-squares, random steps, Monte Carlo inspired techniques, etc.) and diagonalize the matrix with new coefficients.

BasMat Library Description
The program relies on a library of matrices that have been calculated with another Mathematica notebook for a given value of the maximum number of phonons, namely N max = 40.The N40 folder contains subfolders, one for each set of values of the angular momentum, namely L = 0, 2, 3, 4, 5, 6, 7, 8 that allow going up to seniority = 5.For each pair of values {N max , L}, there is thus a subdirectory where a certain number of matrices are stored in Harwell-Boeing format (.rha corresponding to real hermitian).This is a convenient and widely used format for storing sparse matrices.The library has been generated with another Mathematica notebook (not published) on a 24 processors Intel(R) Xeon(R) E5-2440 machine with 48 Gbyte RAM and Linux O.S.Although it is not easy to estimate the timing, due to multiple users and several processes running, it has taken about 2 h for the smallest matrix (L = 3 with dim = 176) and more than a week for one of the largest (L = 8 with dim = 771).Apparently, the time scales in a difficult-to-predict way that is not only controlled by the basis dimension, but also by the time it takes to calculate to a given precision integrals involving wavefunctions (polynomials) with growing value of L. The library matrices are read by the program depending on the chosen level of convergence.Higher maximum phonon numbers insure higher convergence, but at the price of higher computational costs.Therefore, it might be better to set smaller values of N max for trials and go for the higher values only at the moment of getting results.In those cases, although all matrix elements are read, only a submatrix corresponding to a smaller phonon number is used in computations.The preliminary readings and loading of databases is done in Part-1 of the code.

Control of "Border Effects" Due to the Truncation of the Matrices
There is an artificial "border effect" due to the truncation of the basis to a finite number of phonons.Let us consider the action of any operator expressed as a certain power of β on a given state labeled by N. Since β can be expressed as a linear combination of phonon creation and annihilation operators, where ω 2 = √ C 1 /B 2 (having used the notation introduced in Equations ( 4) and ( 6) (not to be confused with the fact that in the original Bohr Hamiltonian C 2 is usually the coefficient of the harmonic oscillator potential)), then the exponent of β k corresponds either to the (maximum) number of phonon added or subtracted to the state.This gives an immediate selection rule on the number of phonons: every pair of states that differ for more than k will give a null matrix element.The selection rule allows an important optimization, because the calculation of many matrix elements can be skipped at once.On the other side, it clarifies why a naive diagonalization using matrices of the same dimension of the basis is affected by numerical errors.On the last few rows and columns (depending on the order of the operators we are considering), the matrix elements will have many connections with states with a higher number of bosons.Therefore, it is a good practice to enlarge the number of basis states that are used in computations up to what is needed according to the selection rule from N max to N max + 6 so that the N max × N max matrix will have all elements correctly calculated.We normally take six more phonons because this value is valid for all operators appearing in the potential (Equation ( 6)).Table 1 shows the dimensions of the extended and trimmed subbases for each angular momentum used in the program.

Scaling Factor in the Radial Variable
As shown in Ref. [27] and subsequently used in various calculations, it is convenient to introduce a scaling parameter, s that minimizes the eigenvalues.The expression for the radial wavefunctions leads to Consequently, the probability density can be squeezed or stretched along β with a suitable value of s.In particular, one generally wants to attain small values of β 2 and therefore a somewhat high value of s is desirable.This can be achieved because of the property of the matrix elements [33]: When N max is small, the dependence of the eigenvalues on s is sizable, while it disappears in the limit N max → ∞.The best estimate for an energy eigenvalue is obtained by minimizing E i (s) with respect to s, but using a different scaling parameter for each eigenvalue would be impractical.A solution is found in searching for an optimal scaling parameter that minimizes the difference between the wavefunctions of the original Hamiltonian and the solution in a truncated basis.Our strategy for calculating s is based on Ref. [33] and implemented in the sfind[] function that takes, as input, some of the coefficients.
The rescaled Hamiltonian is written as: where C 0 = h2 2B .

Model Calculations and Examples
The two files coeff-habs.nb and coeff-fortunato.nbcontain "best" sets of coefficients obtained from literature in the first case or from extensive trials with this program on state-of-the-art high-performance cpu's in the second case.Of course, one can never be sure that the algorithm has hit the absolute minimum and these sets can be probably improved by running more and more random walk fitting procedures (for example, with a faster computer), changing the initial conditions, etc.Another source of discrepancies might be that some of the energies have been slightly re-determined since the times when the paper by Habs et al. was written.
The file data-sets.nbcontains a few sets of experimental data that can be called in connection with the coefficients of the previous paragraph.
As a ready-made example of calculation, we propose, in Part-2 of the code, the spectrum of 128 Xe as obtained from the list of coefficients.This can be easily changed or adjusted by the user.After reading the input parameters from the database, the program computes the spectrum and gives the output reported in Figure 1 Here, as a function of J π , we can see absolute and relative energies (in MeV) and reduced energies, obtained by normalizing the first 2 + state to 1 with the following formula: where E(J π ) are the absolute energies in MeV of the state with spin-parity J π .The program also allows calculating wavefunctions (this is rather slow, because it is a summation over many basis states) and transitions rates and can be easily generalized depending on the needs.

Playing Area and Fitting Subroutines
In Part-3 of the code, one finds a testing function, i.e., a subroutine that calculates the sum of square of deviations between theory and input experimental data to be found in the data-sets.nbfile.If a particular dataset is wanted, but not contained in this file, it can be easily added by the user.The type and number of states to be compared can easily be changed by adapting the last line of the f[] function.In the present version, the first 0 + , 2 + , 4 + and 6 + of the ground state band;, the 2 + , 3 + and 5 + of the gamma band; and the 0 + of the beta band are used.Notice that the 3 and 5 states are rather isolated (other occurrences of these values of J appear only at high energy), which is why they are used in place for example of the 4 + 2 that is amidst or very close to other 4 + states and might be highly mixed with them.Often, one has to guess which excited 0 + state might be the bandhead of the beta band, and it might happen that is not always the lowest.Then, a random fitting procedure is given: starting from initial foul's parameters, it performs a random walk in parameter space trying to find the minimum.Normall,y a few runs are enough to identify deep local minima, although with this technique one cannot be insured to have hit the lowest absolute minimum.The suggestion here is to run the routine from scratch a number of times (by experience, about 10 or 20 times) and try to seek the lowest minimum.When this is done, one can run a second minimization routine with a smaller variation of parameters to slowly minimize the sum of squared deviations.The final results are usually very good, although there might be cases where the whole procedure does not work.Finally, one can draw the potential energy surface as a contour plot or in 3D.We show in Figure 2 two examples, one for a spherical minimum and one with two almost coexisting minima in 180 Pt (notice that the deformed minimum is lower).The program adjusts itself to plot things around the minimum, but of course the user can force any range upon the plotting routine.The program also locates the minimum of the energy surface, more than one minimum is present they must be searched after inspection, by changing the starting point in the FindMinimum[] command line.

Installation Tips
Just unfold the package with the files GCMsolver_3.nb,basismodules.nb,data-sets.nb,coeff-habs.nb and coeff-fortunato.nb(Supplementary Materials) in a directory and run the first under Mathematica.There might be an issue in setting directory and subdirectory names and paths, namely certain systems use the slashes and others backslashes.The user should adjust that according to the system in use.

Conclusions
We wrote [14] an easy-to-use program that serves two scopes: (1) as a didactic tool to learn the basics of the Bohr-Mottelson model, by solving the Bohr Hamiltonian in the formulation of Gneuss and Greiner (i.e., with potential terms coming from the scalars that can be built with the quadrupole tensors up to the sixth order in β); and (2) as a research tool, because the Bohr Hamiltonian has returned prominently to the forefront of nuclear physics thanks to the new solutions of Iachello that give a paradigm for nuclei sitting at the critical point of shape phase transitions between different shapes.The code was written in Mathematica [13].Because the interpretation of the code is immediate, it allows ample space for improvements and customization and the built-in visualization routines easily allow generating and exporting graphs and data.

Figure 2 .
Figure 2. Output of the contour plot drawing subroutine for 128 Xe from the habs[] database and 180 Pt from the fort[] database.The plots are given as a function of (β, γ).
[7,24]l the functional form of the potential energy surface.This potential is at the core of the Gneuss-Greiner generalized Collective model.Clearly, certain special cases are recovered by setting to zero certain coefficients: for example when all coefficients except the first are zero, we recover the original Bohr-Mottelson model for vibrational nuclei.The pure quartic solution, which plays a prominent role at the critical point of spherical to γ-unstable shape phase transitions[7,24], is found when only C 3 is non-null.Axial prolate or oblate shapes require C 2 or C 4 different from zero.Triaxial shapes (with a minimum occurring for 0

Table 1 .
Dimension of extended subbasis (second line) and trimmed subbasis (last line) for angular momentum used in the program (first line).
. Output of the diagonalization routine for 128 Xe.The first panel shows absolute energies in MeV, the second shows energies relative to the ground state and the third shows reduced energies, i.e., the energies of the second panel normalized with the energy of the first 2 + state.