Disentangling phase transitions and critical points in the proton-neutron interacting boson model by catastrophe theory

We introduce the basic concepts of catastrophe theory needed to derive analytically the phase diagram of the proton-neutron interacting boson model (IBM-2). Previous studies [1,2,3] were based on numerical solutions. We here explain the whole IBM-2 phase diagram including the precise order of the phase transitions in terms of the cusp catastrophe.


Introduction
In the last twenty years quantum phase transitions (QPTs) (phase transitions that happen at zero temperature as a function of a control parameter) have been a subject of great interest in different areas of quantum many-body systems. In particular, in Nuclear Physics the study of shape phase transitions is a topic of current interest both theoretically and experimentally [4][5][6]. Moreover, in Molecular Physics [7,8], Quantum Optics [9,10] and Solid State Physics [11] the interest on QPTs has grown enormously in recent years.
Strictly speaking, QPTs take place for large systems in the thermodynamic limit as a discontinuity or singularity in some derivative of the ground state energy. However, finite systems like the atomic nuclei could show the precursors of a phase transition when structural changes in the ground state are observed as a function of the neutron or proton numbers (N, Z ) [12]. These transitional nuclei are characterized by specific patterns in the low lying spectrum which could be associated with a critical system. Recently, in the context of the Bohr Hamiltonian [13], F. Iachello has introduced the concept of critical point symmetry [14][15][16] that can be applied when the quantum system undergoes transitions between two phases with different shape.
A convenient model to study QPTs in nuclei is the interacting boson model (IBM) [17]. It is a symmetry-dictated model whose building blocks are ideal bosons representing nucleon pairs with angular momentum zero, s bosons, and two, d bosons. In its simplest version, IBM-1, the dynamical algebra of the model is U (6) and the common symmetry algebra O(3). IBM-1 presents three dynamical symmetries (SU (5), O (6), and SU(3)) corresponding to well defined nuclear shapes (spherical, γ -unstable, and axially deformed, respectively). The presence of different dynamical symmetries in the model is a key ingredient to study QPTs, since they appear when two different symmetries are mixed in the Hamiltonian through a control parameter: . For a particular value of the control parameter, ξ c , the system undergoes a structural QPT from symmetry 1 to symmetry 2.
The geometry and shape phase transitions of the IBM-1 were studied long ago [18][19][20] and also more recently the many facets of QPTs in IBM have been analyzed [5]. A similar kind of study could be extended to the more realistic version of the IBM that includes the proton-neutron degree of freedom, known as IBM-2 [21]. Most of the studies [1][2][3] rely on the numerical diagonalization of the IBM-2 Hamiltonian which is restricted to systems with, as maximum, 10 proton and 10 neutron bosons, or in the numerical treatments of semiclassical approximations. These numerical studies, that provided a very accurate description of the  phase diagrams, could not state in an unambiguous way the classification of phase transition orders. In this Letter we present an analytic study of the IBM-2 phase diagram that is able to determine unambiguously the order of the phase transitions present in the two-fluid IBM-2 model making use of the catastrophe theory (CT) [22][23][24]. CT allows analyzing energy surfaces depending on various parameters, i.e., families of potentials, that contain several shape variables, providing information on the nature and numbers of minima (also maxima, if they exist), i.e., stability, depth, etc. This analysis leads to a partition of the parameter space into different regions where the energy surface has different qualitative properties (number and nature of maxima and minima) and, in particular, it serves to study and classify the existing QPTs. Similar studies were carried out for IBM-1 [25], for IBM including configuration mixing [26,27], for three-component thermodynamic systems [28,29], and for elementary chemical reactions [30].

IBM-2 Hamiltonian
In this work we use a simplified form of the IBM-2 Hamiltonian, known as the Consistent-Q Hamiltonian [31], which retains all the main ingredients of the full Hamiltonian where μ and ρ = π, ν. N = N π + N ν is the total number of bosons representing the number of valence nucleon pairs.
We study the phase diagram of the IBM-2 in the semi-classical or mean field formalism. In this approach, which is exact in the thermodynamic limit, the ground state wavefunction is a product of a proton condensate times a neutron condensate [32,33], |g = |N π , N ν , β π , γ π , β ν , γ ν where |0 is the boson vacuum and Γ † ρ is the creation operator for a coherent π or ν boson, defined as The equilibrium values of the structure parameters (β π , γ π , β ν , γ ν ) and the energy of the system for given values of the control parameters (ξ , χ π , χ ν ) in the Hamiltonian (1) can be obtained by minimizing its expectation value in the intrinsic state (2): δ g|H|g = 0. In the limit N π , N ν → ∞ the energy per boson can be obtained in a straightforward way as E(β π , γ π , β ν , γ ν ; χ π , χ ν , ξ) Notice that Hamiltonian (1) can only lead to energy surfaces where proton and neutron ellipsoids are axially symmetric with symmetry axis either parallel or perpendicular, or to triaxial shapes in both ellipsoids (see [34] for details). Therefore, it is not required to include explicitly the Euler angles.
In order to analyze the energy surface it is convenient to introduce new variables, In Fig. 1 we depict the phase diagram of the model, which has been obtained numerically [1][2][3]. The different phases correspond to the spherical region S, the prolate axially deformed region P , the oblate axially deformed region O , and the region with triaxial shapes T . It is interesting to note that the IBM-1 phase diagram is recovered for χ = 0. The phase diagram could also be extended to positive values of χ , by reflection in the horizontal plane.
The orders of the different phase transitions were inferred from the numerical calculations in [1][2][3]. All references coincide in the following classification: • The x-x * -e-x surface: (equally the surface x-x * -e-x in the extended diagram to include oblate shapes) is first order, except for the x * − e line, which was proposed to be second order; In what follows we will introduce the basic ingredients of CT to determine in an unambiguous way the properties of the QPTs along the critical lines/surfaces in the IBM-2 phase diagram.

Catastrophe theory program
Once the IBM-2 phase diagram is known through a numerical study [1][2][3], for the conclusive determination of the order of the phase transitions associated to the obtained critical lines/surfaces it is required an analytic study. CT is specially suited to carry out such analysis.
The basic goal of CT is the study of a potential or, more generically, a family of potentials. As we will see, these potentials comprise three types of points. To begin with, let us assume a system described by a real family of smooth potentials: where x ∈ n stand for the state (order) variables and λ ∈ r are the control parameters. The majority of points x have a nonzero gradient and are called regular points. The points where the gradient vanishes are called stationary or critical points and they can be classified in two groups: (i) the points where the determinant of the Hessian matrix is different from zero, called isolated, non-degenerated or Morse points, and (ii) the points where the determinant of the Hessian matrix is zero, called non-isolated, degenerated or non-Morse points. In summary, points of a family of smooth potentials can be classified according to their gradient and Hessian matrix H as: Morse theorem [23,24] guarantees that around a Morse point, a smooth potential is equivalent to a quadratic form, thanks to a smooth non-linear change of variables. Therefore, the stability of the potential under small perturbation in the parameters is guaranteed in Morse points. At non-Morse points the potential cannot be written as a quadratic form because the Hessian matrix has at least one zero eigenvalue. It is at non-Morse points where CT shows all its power. Let us illustrate the CT program starting with the following expression (see Ref. [24] for a more detailed description), where x 0 is a degenerated critical point (non-Morse) for the control parameter λ 0 . We perform a Taylor expansion in the order parameters till k-order. The problem of determinacy consists in getting the Taylor series that can be truncated without loss of substantial information with respect to the original function. The issue of finding the most general family of functions with the smallest dimension, d, which contains the original function, is known as unfolding. The number of parameters appearing in this unfolding is called codimension or number of essential parameters. This number is connected with the number of lowest order terms in the Taylor expansion that can be canceled out. The so-called catastrophe germ is obtained when all the unfolding terms go to zero, i.e., all possible terms in the Taylor expansion vanish, The final concept to be introduced is transversality. One says that the original function, V is k transversal when it is isomorphic to the canonical form of the unfolding [23].
Thom's splitting lemma [22] guarantees that a smooth potential at non-Morse points can be written as the sum of a quadratic form, associated to the subspace with nonzero eigenvalues, plus a function containing the variables associated to the zero eigenvalues of the Hessian matrix. The non-Morse part of the latter is a canonical form called catastrophe function. This function is composed by the catastrophe germ, which only depends on the number of vanishing eigenvalues and on the number of control parameters, and by a universal perturbation that removes the degeneracy and makes the potential structurally stable. The catastrophe germs and the related universal perturbations were listed by Thom for potentials up to two variables and up to five parameters [23,24]. The transformation into the canonical form only exists for a function which is k transversal. If this is not the case, the function cannot be treated with CT.
The first step in the CT program is to find out the critical points of the potential (∇ V = 0). Among them, the most important is the most degenerate one. This point is the fundamental root taking place at a definite value of the control parameters which we will call critical values. We next proceed making use of a Taylor expansion of the potential around the fundamental root. A Taylor expansion around such a point is also valid for the critical points that arise from the fundamental root when the degeneracy is broken. Depending on the degeneracy of the fundamental root the number of extremes that can be analyzed simultaneously will change.
It is important to note here that if the original function is k determined, then the k-order Taylor expansion can be transformed, under the appropriate non-linear change of variables, into a finite polynomial that is valid in the neighborhood of the critical point and critical control parameter, (11) where x i stand for the original and x i for the transformed variables. Finally, this polynomial can be written in a canonical form composed by the catastrophe germ plus a generic perturbation. Once the properties of the catastrophe function have been established, the stability of the potential is completely determined. In CT the set of non-Morse points is known as bifurcation set while the ensemble of critical points with equal energy are known as the Maxwell set. If the potential is in the neighborhood of a Maxwell set, it possesses a first order phase transition. Conversely, if the potential is close to a bifurcation set, it has a second order phase transition.
When the potential depends on several variables it is important to find out which are the variables involved in the phase transition, since they play the role of order parameters. These variables are associated to the subspace with vanishing Hessian eigenvalues, called bad or essential variables, while there is another set of variables related to the non-vanishing Hessian eigenvalues, called good or non-essential variables. The potential could be separated into a part depending on the essential variables and into another part depending on the non-essential ones by rewriting it in terms of the eigenvectors of the Hessian matrix. As a consequence of the splitting lemma the potential is separated into a function of the essential variables and a sum of quadratic terms associated to the non-essential variables [24]. Therefore, the appearance of critical phenomena is associated exclusively with the behavior of the essential variables.
In next section this program is applied to the two-fluid nuclear IBM-2. It will be shown that the relevant elementary catastrophe for this model is the cusp catastrophe ( A +3 ) which is characterized by one state variable (z) and two-control parameters (a, b). The germ of this catastrophe is z 4 and the perturbation az + bz 2 .

Application of the catastrophe theory program to IBM-2
Even for the restricted Hamiltonian (1) it is not possible to carry out a general analysis of the whole phase diagram due to the large number of shape variables (four in the case of the IBM-2). Moreover, it is a nontrivial task the identification of the appropriate order parameters. In order to proceed with the analysis we will concentrate in the critical surfaces depicted in Fig. 1, which were already studied numerically in [1][2][3].

Spherical-axially deformed energy surface
This surface is marked with points x-e-x-x * -x in Fig. 1. It corresponds to a situation in which γ π ,ν can be assumed to be zero, because in the spherical phase the energy does not depend on γ , and in the axially deformed phase γ = 0 (or γ = π/3 in the oblate side). This assumption is not valid in the case of the line e − x * as it will be explained later on.
Following the procedure described above, we use the fundamental root corresponding to β π = β ν = 0 and γ π = γ ν = 0 for χ < 0 (γ π = γ ν = π/3 for χ > 0), to construct the Hessian matrix associated to Eq. (4) The two eigenvalues are 5ξ − 4 and ξ , and the corresponding eigenvectors are The eigenvalue associated to β 1 vanishes for ξ = 4/5 while the one associated to β 2 only vanishes for the trivial case ξ = 0. Therefore the essential variable turns out to be β 1 , while β 2 becomes the non-essential one. If we make an expansion of the energy in terms of β 1 and β 2 we get: where the terms Θ(β 5 1 ) and Θ(β 1 β 2 2 , β 2 β 2 1 ) can be canceled through a nonlinear transformation (11) in the non-essential variable. Notice that the expansion (15) has no quadratic term proportional to β 1 β 2 since β 1 and β 2 are eigenvectors of the Hessian matrix. Therefore, E simplifies to In order to keep the notation simple we keep the variable β 2 , although there, it corresponds to the transformed variable (see Eq. (11)). The most salient feature of Eq. (16) is the existence of a cubic term, which guarantees that the phase transition around ξ = 4/5 will always be of first order if χ = χ π +χ ν 2 = 0 [24].
Eq. (16) can be transformed into the cusp catastrophe for which a first order phase transition exists if the linear term is different from zero, which is indeed the case for χ = 0.
Following Ref. [25], the critical value of the control parameter ξ c is the solution of the equation: where and which leads to the solution: This expression gives the well known values ξ c = 4/5 for χ = 0 and ξ c = 9/11 for χ = ± √ 7/2, which are also valid for IBM-1. It is important to note that (20) is independent on χ which implies that the first order IBM-1 critical line propagates vertically generating a first order critical surface separating spherical and axially deformed shapes in IBM-2, as was already established in [3] using different arguments.

The e-x * line
This line is a limiting case of the previous energy surface, defined by χ = 0 (χ π = −χ ν ). The line has some important differences, which deserve a particular analysis. Along this line β π = β ν = β. Moreover, γ π = π/3 − γ ν = γ cannot be zero because one of the phases is triaxial while the other is γ independent.
These conditions define the fundamental root be β π = β ν = β = 0, γ π = γ ν = γ = π/6. The Hessian matrix evaluated at the fundamental root can be written as The eigenvalue associated to the variable β is 5ξ − 4 and it is canceled for ξ = 4/5. All derivatives in γ at the fundamental root vanish. Thus β will be the essential variable and γ the non-essential one.
The expansion around the fundamental root gives where the odd powers vanish. Since Eq. (22) does not depend on χ and it has no cubic term, the whole line e-x * will be second order [24].

Axially deformed-triaxial surface
This surface is delimited by the points e-O (6)-y-x * -e (e-O (6)y-x * -e in the oblate side) in Fig. 1. It represents the most complex situation because the four shape variables (β π , γ π , β ν , γ ν ) have to be treated simultaneously. With the exception for the y point, the y-O (6) line, and the e-O (6) line, which will be treated separately, no simplification is possible.
The fundamental root corresponds to γ π = γ ν = 0 (γ π = γ ν = π/3 for χ > 0). The critical values of β π and β ν depend on the value of the Hamiltonian parameters, therefore we will proceed to study the Hessian matrix for γ π = γ ν = 0. The main feature of this matrix is that for γ π = γ ν = 0 (γ π = γ ν = π/3 for χ > 0) where ρ and ρ stand for π, ν. This equation implies that β's are decoupled from the angular variables. As shown in Ref. [1] the behavior of β variables when crossing the axially deformed-triaxial surface is smooth and they are related to the non-vanishing eigenvalues of the Hessian matrix. Therefore, the β's are non-essential variables. According to the splitting lemma the energy surface can be expanded as a quadratic form in β's (except in the line e-x * already discussed) plus an expansion in the γ variables. Note that the coefficients will depend on the equilibrium values of β π and β ν , i.e., β 0 π and β 0 ν ,