1. Introduction
Most atomic nuclei are deformed in their ground states [1, 2]. The nuclear shapes can be described by the surface multipole expansion
$\begin{eqnarray}R(\theta ,\varphi )={R}_{0}\left[1+{\beta }_{00}+\sum _{\lambda =2}^{\infty }\sum _{\mu =-\lambda }^{\lambda }{\beta }_{\lambda \mu }^{* }{Y}_{\lambda \mu }(\theta ,\varphi )\right],\end{eqnarray}$
where R is the distance from the center of mass to the nuclear surface, and λ and μ are integers characterizing different kinds of multipole deformations. For example, βλ0 describes axially symmetric deformations, the triaxial deformations correspond to βλμ with even λ, and the tetrahedral shape means β32 ≠ 0 while all other βλμ vanish. These shape parameters can be viewed as generalized coordinates, from which we can build dynamical theories describing the rotations and vibrations of the nuclei. These collective motions are foundations of nuclear spectroscopy.The nuclear spectra are the fingerprints of the underlying nuclear shapes. The most studied deformations are the axially symmetric quadrupole deformation with only β20 ≠ 0, which manifests itself in the typical spectra of a symmetric rigid rotor, E ∝ I(I + 1). A more interesting case is the triaxial deformation with non-zero β20 and β22, which results in the wobbling motion [2] and the chiral doublet bands [3–5]. The multiple chiral doublet (MχD) in nuclei with different triaxial configurations was predicted by Meng et al [6] and recently confirmed experimentally in several nuclei [7–11]. A number of studies have shown that the low-spin signature inversion indicates the existence of triaxiality in certain regions of the nuclide chart [12–17]. On the other hand, the parity-breaking octupole deformations β3μ ≠ 0 are also found in many experiments [18–20]. For example, the low-lying alternative parity bands found in actinides and some rare-earth nuclei hint the possible existence of reflection asymmetric shapes [21–28]. Direct evidence of a static reflection asymmetric shape was found in 224Ra [29], 144Ba [30] and 228Th [31].
In recent years various exotic nuclear deformations have attracted a lot of theoretical interest. For example, reflection asymmetric deformations combined with triaxial deformations have been proposed to influence the nuclear fission barrier [32–36]. Nuclear clustering has been proved to be important in light nuclei, especially in the states near the nα binding threshold [37–40]. Nuclear clustering configurations with various arrangements of α clusters are naturally accompanied by exotic shape deformations. It implies that a proper theoretical description requires simultaneous inclusion of various shape degrees of freedom. It has been demonstrated that the density functional theories are able to describe the clustering structure in light nuclei [41–49]. The study in [41] discussed the clustering structure of light nuclei with N = Z using the axially deformed relativistic mean-field model. Using the cranking covariant density functional theory (CDFT), rod-shaped carbon isotopes at high spins were studied [48, 50, 51]. Studies based on the algebraic model [52] and the antisymmetrized molecular dynamics model [53] suggested that the α clustering can appear in the low-lying states of 12C and other light even–even nuclei [54–56]. To study how nuclei cluster and how the clustering influences nuclear properties, we need theories that can describe both the excited states and the ground state, and include as many shape degrees of freedom as possible.
The CDFT gives a universal description of the nuclei all over the nuclide chart [57–65]. In CDFT, the Lorentz invariance nature of the Lagrangian automatically determines the spin–orbit interaction, and reveals the origin of the pseudospin symmetry [64, 66–71] in the Fermi sea and the spin symmetry in the Dirac sea [72, 73]. An advantage of studying nuclear clustering with CDFT is that it does not presume the existence of α clusters a priori but generates such configurations from microscopic interactions. However, the mean-field approximation used in CDFT brings the problem that the intrinsic symmetries and conservation laws, such as the translational invariance, rotational invariance, parity and particle number conservations, are usually broken [74]. One solution is to restore these symmetries using the projection method; examples are the angular momentum projection (AMP) and the parity projection (PP) methods [62]. By combining AMP and PP with the generator coordinate method, we can calculate the rotational–vibrational spectra, which can be directly compared with the experiments [62, 75–78].
In recent years, multidimensionally constrained CDFTs (MDCCDFTs) [33, 79–81] have been established and applied successfully in several topics. In these models, the reflection asymmetric and axially asymmetric deformations can be considered simultaneously, with only a spatial V4 symmetry remaining unbroken for the densities and potentials. Depending on how to deal with the pairing correlation, two variants were developed: the multidimensionally constrained relativistic mean-field model [33, 79], and the multidimensionally constrained relativistic Hartree–Bogoliubov model (MDCRHB) [80]. MDCCDFTs have been applied to study the potential energy surfaces (PES) and fission barriers in actinides [32–34, 79] and superheavy nuclei [82], the third minima on PESs of light actinides [35], higher-order deformations in superheavy nuclei [83], non-axial octupole Y32 correlations in N = 150 isotones [84] and Zr isotopes [80] and axial octupole Y30 correlations in MχD [11]. Fission dynamics in actinides were studied based on PESs calculated with the MDCCDFTs [85–89]. Furthermore, by including the strangeness degrees of freedom, MDCCDFTs were employed to study the shape evolution [90, 91] and hyperon pairing correlation [92] in hypernuclei, and to build new nucleon–hyperon effective interactions [93].
In this work we focus on the MDCRHB model and its extensions. The mean-field state $\left|\phi (q)\right\rangle $ can be obtained by solving the relativistic Hartree–Bogoliubov (RHB) equation in the intrinsic frame. Here, q is a collection of quantum numbers specifying the deformation of the state or other intrinsic properties. Several intrinsic symmetries may be broken due to the mean-field approximation and need to be restored using the symmetry projection methods. For example, the AMP for a relativistic mean-field model has been implemented upon a triaxially deformed ground state [94] and an axially deformed reflection asymmetric ground state [95, 96]. The generating coordinate method (GCM) has also been included to treat soft PESs [75, 76]. In this work, we incorporate the AMP and PP methods into the MDCRHB model. Our implementation differs from the previous works in two aspects. First, we apply a unified framework to treat both the triaxial and octupole deformations. While in [94] the triaxial deformations are considered with an anisotropic 3D harmonic oscillator (3DHO) basis, in this work we employ an axially symmetric harmonic oscillator basis in all cases. This setting allows us to simultaneously break as many symmetries as possible and consider exotic shapes like the triangle or tetrahedron. Second, we use the Bogoliubov transformation together with a separable pairing force of finite range to treat the pairing effect. This makes our results more reliable for extreme deformations for which the pairing force can play key roles in cluster formations.
With the projected MDCRHB model (p-MDCRHB), we are able to explore various exotic nuclear shapes. One direct application involves searching for possible α-clustering structures in light nuclei. As the first step toward this direction, in this work we only consider the AMP and PP methods. The particle number projection method and GCM also play essential roles in some circumstances and will be considered in future works. In section 2 , we present the theoretical framework of the p-MDCRHB model. In section 3 , we present the numerical check and the results for the α-conjugate nucleus 12C. We summarize this work in section 4 .
2. Theoretical Framework
CDFTs solve the nuclear many-body problems within the framework of the density functional theories. We design a general energy density functional E[ρ(r)], where ρ(r) is the intrinsic nucleon densities, then the ground state energy and densities are given by the variational principle δE/δρ = 0. In most CDFTs this framework is formally equivalent to a mean-field approximation for an effective Lagrangian, in which the nucleons are Dirac fermions interacting indirectly via meson exchange (ME) or directly through zero-range point-coupling interactions [97, 98]. To give the correct saturation properties of the symmetric nuclear matter, the nonlinear self-coupling terms [99–101] or the density dependence of the coupling constants [102, 103] were introduced and applied extensively. Accordingly, the framework of CDFT can be formulated using four representations: effective interactions with ME, or point-coupling (PC) combined with nonlinear self-coupling (NL) or density dependence of coupling constants (DD). In MDCCDFT we implemented all four possibilities. In this work, we illustrate the formalism of the projected MDCRHB model using the density-dependent point-coupling interaction DD-PC1 [104]. The formalisms for other parameters are similar and omitted for brevity.
2.1. Relativistic Hartree–Bogoliubov model
In this section we briefly introduce the theoretical framework of the MDCRHB model with density-dependent point-coupling interactions. We start with the effective Lagrangian
$\begin{eqnarray}\begin{array}{l}{ \mathcal L }=\bar{\psi }({\rm{i}}{\gamma }_{\mu }{\partial }^{\mu }-M)\psi \\ -\displaystyle \frac{1}{2}{\alpha }_{S}(\rho ){\rho }_{S}^{2}-\displaystyle \frac{1}{2}{\alpha }_{V}(\rho ){j}_{V}^{2}-\displaystyle \frac{1}{2}{\alpha }_{{TV}}(\rho ){\vec{j}}_{{TV}}^{2}\\ -\displaystyle \frac{1}{2}{\delta }_{S}({\partial }_{\nu }{\rho }_{S})({\partial }^{\nu }{\rho }_{S})\\ -{{eA}}_{\mu }{j}_{C}^{\mu }-\displaystyle \frac{1}{4}{F}^{\mu \nu }{F}_{\mu \nu },\end{array}\end{eqnarray}$
where M is the nucleon mass, and αS(ρ), αV(ρ) and αTV(ρ) are density-dependent coupling constants in scalar, vector and isospin-vector channels, respectively. Here, ρS, jV and ${\vec{j}}_{{TV}}$ are the corresponding densities and four-currents, δS is the coupling constant for the derivative term, ${j}_{C}^{\mu }$ is the proton current, Aμ is the electromagnetic field and Fμν = ∂μAν − ∂νAμ is the field strength tensor. For the DD-PC1 functional [104] the functions αS(ρ), αV(ρ) and αTV(ρ) share the same ansatz $\begin{eqnarray}{\alpha }_{i}(\rho )={a}_{i}+({b}_{i}+{c}_{i}x){e}^{-{d}_{i}x}\quad (i=S,V,{TV}),\end{eqnarray}$
with x = ρ/ρsat, where ρsat denotes the saturation density of the symmetric nuclear matter. The constants ai, bi and ci are parameters fitted to the experiments.The pairing correlations are included with the Bogoliubov transformation. According to the mean-field approximation, we take the Bogoliubov vacuum as the trial wave function and solve the variational problem to find the ground state. To prevent the system from collapsing to the negative energy Dirac sea, we apply the ‘no-sea' approximation, i.e. only the positive energy states are kept when calculating the densities and currents. The resulting RHB equation [58, 105] in coordinate space is7 ) and what follows we use the symbols with and without prime to denote the quantities for the initial and final states, respectively. As a standard procedure in RHB models [106], we only include the large component of the Dirac spinors when calculating the pairing tensor.
$\begin{eqnarray}\begin{array}{l}\displaystyle \int {{\rm{d}}}^{3}{\boldsymbol{r}}^{\prime} \left(\begin{array}{cc}h-\lambda & {\rm{\Delta }}\\ -{{\rm{\Delta }}}^{* } & -h+\lambda \end{array}\right)\left(\begin{array}{c}{U}_{k}\\ {V}_{k}\end{array}\right)\\ ={E}_{k}\left(\begin{array}{c}{U}_{k}\\ {V}_{k}\end{array}\right),\end{array}\end{eqnarray}$
where Ek is the quasi-particle energy and λ is the chemical potential. Here, ${U}_{k}={U}_{k}({{\boldsymbol{r}}}^{{\prime} }{\sigma }^{{\prime} })$ and ${V}_{k}={V}_{k}({{\boldsymbol{r}}}^{{\prime} }{\sigma }^{{\prime} })$ are quasi-particle wave functions; ${\rm{\Delta }}={\rm{\Delta }}({\boldsymbol{r}}\sigma ,{{\boldsymbol{r}}}^{{\prime} }{\sigma }^{{\prime} })$ is the pairing potential and h is the single particle Hamiltonian $\begin{eqnarray}\begin{array}{l}h={\boldsymbol{\alpha }}\cdot [{\boldsymbol{p}}-{\boldsymbol{V}}({\boldsymbol{r}})]+\beta [M+S({\boldsymbol{r}})]\\ +{V}_{0}({\boldsymbol{r}})+{{\rm{\Sigma }}}_{R}({\boldsymbol{r}}),\end{array}\end{eqnarray}$
where the scalar potential S(r), vector potential Vμ(r) = (V0(r), V(r)) and the rearrangement term ΣR(r) are expressed with various nuclear densities and currents, $\begin{eqnarray}S={\alpha }_{S}(\rho ){\rho }_{S}+{\delta }_{S}{\rm{\Delta }}{\rho }_{S},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{V}^{\mu }={\alpha }_{V}(\rho ){j}_{V}^{\mu }+{\alpha }_{{TV}}(\rho ){j}_{{TV}}^{\mu }\cdot \vec{\tau } +e\displaystyle \frac{1-{\tau }_{3}}{2}{A}^{\mu },\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{{\rm{\Sigma }}}_{R}=\displaystyle \frac{1}{2}\displaystyle \frac{\partial {\alpha }_{S}}{\partial \rho }{\rho }_{S}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{\partial {\alpha }_{V}}{\partial \rho }{j}_{V}^{2}+\displaystyle \frac{1}{2}\displaystyle \frac{\partial {\alpha }_{{TV}}}{\partial \rho }{\vec{j}}_{{TV}}^{2}.\end{array}\end{eqnarray}$
The pairing potential Δ is written as $\begin{eqnarray}\begin{array}{l}{{\rm{\Delta }}}_{{p}_{1}{p}_{2}}({{\boldsymbol{r}}}_{1}{\sigma }_{1},{{\boldsymbol{r}}}_{2}{\sigma }_{2})\\ =\displaystyle \int {{\rm{d}}}^{3}{{r}}^{{\prime} }_{1}{{\rm{d}}}^{3}{{r}}^{{\prime} }_{2}\displaystyle \sum _{\sigma {{\prime} }_{1}\sigma {{\prime} }_{2}}^{p{{\prime} }_{1}p{{\prime} }_{2}}{V}_{{p}_{1}{p}_{2},{p}_{1}^{{\prime} }{p}_{2}^{{\prime} }}^{\mathrm{pp}}({{\boldsymbol{r}}}_{1}{\sigma }_{1},{{\boldsymbol{r}}}_{2}{\sigma }_{2},{\boldsymbol{r}}{{\prime} }_{1}\sigma {{\prime} }_{1},{\boldsymbol{r}}{{\prime} }_{2}\sigma {{\prime} }_{2})\\ \ \,\times {\kappa }_{{p}_{1}^{{\prime} }{p}_{2}^{{\prime} }}({\boldsymbol{r}}{{\prime} }_{1}\sigma {{\prime} }_{1},{\boldsymbol{r}}{{\prime} }_{2}\sigma {{\prime} }_{2}),\end{array}\end{eqnarray}$
with the pairing tensor $\kappa =\langle \hat{a}({{\boldsymbol{r}}}_{1}{\sigma }_{1})\hat{a}({{\boldsymbol{r}}}_{2}{\sigma }_{2})\rangle $. In equation (In this work we use a separable pairing force of finite range [107, 108]
$\begin{eqnarray}\begin{array}{l}V({{\boldsymbol{r}}}_{1}{\sigma }_{1},{{\boldsymbol{r}}}_{2}{\sigma }_{2},{{r}}^{{\prime} }_{1}\sigma ^{{\prime} }_{1},{{r}}^{{\prime} }_{2}\sigma ^{{\prime} }_{2})\\ \ \,=-G\delta ({\boldsymbol{R}}-{\boldsymbol{R}}^{\prime} )P(r)P(r^{\prime} )\displaystyle \frac{1}{2}(1-{P}_{\sigma }),\end{array}\end{eqnarray}$
where G is the pairing strength, (1 − Pσ)/2 is the projector onto the total spin S = 0, R = (r1 + r2)/2 and r = r1 − r2 are the center of mass coordinate and relative coordinate, respectively. Here, P(r) is a normalized Gaussian function $P(r)={\left(4\pi {a}^{2}\right)}^{-3/2}\exp \left(-{r}^{2}/(4{a}^{2})\right)$, where a is an adjustable parameter specifying the range of the pairing force. The details for calculating the pairing tensor and pairing potential in the axially deformed harmonic oscillator (ADHO) basis can be found in [80]. In this illustrative work we fix the pairing strength to G = 1.1G0 with G0 = 728.0 MeV·fm3 and the range to a = 0.644 fm, in which G0 and a are determined by fitting to the pairing gap in the nuclear matter [107, 109].The RHB equation (4 ) can be solved self-consistently starting from a properly chosen initial state. The resulting minima represent the ground state or the shape isomers. To obtain the PESs, E(β), we use a modified linear constraint method [33]. The Routhian with deformation constraints reads2 ) and the pairing force (8 ). The variable Cλμ varies during the iteration. βλμ is the nuclear deformation
$\begin{eqnarray}E^{\prime} ={E}_{\mathrm{MF}}+\sum _{\lambda \mu }\displaystyle \frac{1}{2}{C}_{\lambda \mu }{\beta }_{\lambda \mu },\end{eqnarray}$
where EMF is the energy functional obtained with the Lagrangian ( $\begin{eqnarray}{\beta }_{\lambda \mu }=\displaystyle \frac{4\pi }{3{{AR}}^{\lambda }}\langle {Q}_{\lambda \mu }\rangle ,\end{eqnarray}$
where R = 1.2A1/3 fm is the approximate nuclear radius, A is the number of nucleons and Qλμ = rλYλμ(Ω), where Yλμ is the spherical harmonics, is the multipole moment operator.2.2. Axially symmetric harmonic oscillator basis
In the MDCRHB model, we solve the RHB equation (4 ) by expanding the spinor wave functions Uk(rσ) and Vk(rσ) in an ADHO basis [33, 80, 110],
$\begin{eqnarray}{U}_{k}({\boldsymbol{r}}\sigma )=\left(\begin{array}{c}\sum _{\alpha }{f}_{U}^{k\alpha }{{\rm{\Phi }}}_{\alpha }({\boldsymbol{r}}\sigma )\\ \sum _{\alpha }{g}_{U}^{k\alpha }{{\rm{\Phi }}}_{\alpha }({\boldsymbol{r}}\sigma )\end{array}\right),\end{eqnarray}$
$\begin{eqnarray}{V}_{k}({\boldsymbol{r}}\sigma )=\left(\begin{array}{c}\sum _{\alpha }{f}_{V}^{k\alpha }{{\rm{\Phi }}}_{\alpha }({\boldsymbol{r}}\sigma )\\ \sum _{\alpha }{g}_{V}^{k\alpha }{{\rm{\Phi }}}_{\alpha }({\boldsymbol{r}}\sigma )\end{array}\right),\end{eqnarray}$
where Φα(rσ) is the basis wave function generated by solving the Schrödinger equation with an ADHO potential. Here, we impose the V4 symmetry, which means that the densities and potentials are invariant under the mirror reflections $\begin{eqnarray}\begin{array}{rcl}{\hat{S}}_{x}\phi (x,y,z) & = & \phi (-x,y,z),\\ {\hat{S}}_{y}\phi (x,y,z) & = & \phi (x,-y,z),\end{array}\end{eqnarray}$
where φ represents any scalar field or spatial distribution in the theory. It immediately follows that the rotation ${\hat{S}}_{z}(\pi )={\hat{S}}_{x}{\hat{S}}_{y}$ around the z-axis by π is also a symmetry operation.In the ADHO basis, the reflection asymmetric shapes can be included by mixing the basis states with different parities. While the axially deformed shapes preserve the z-axis AMP K, the triaxial shapes mix different K blocks and we need to calculate the corresponding mixing matrix elements. In this case, we expand the densities and potentials in terms of the Fourier series of the azimuthal coordinate φ in cylindrical coordinates. Under the restriction of the V4 symmetry, the expansion for any field f writes
$\begin{eqnarray}\begin{array}{l}f(\rho ,\varphi ,z)={f}_{0}(\rho ,z)\displaystyle \frac{1}{2\sqrt{\pi }}\\ +\sum _{n=1}^{\infty }{f}_{n}(\rho ,z)\displaystyle \frac{1}{\sqrt{\pi }}\cos (2n\varphi ).\end{array}\end{eqnarray}$
The ADHO basis is also eigenfunctions of the simplex operator ${\hat{S}}_{z}(\pi )$ with eigenvalues S = ± 1. The basis with S = − 1 can be calculated with a time-reversal operation on those with S = 1. For a system with time-reversal symmetry, the single particle energy levels have Kramer's degeneracy and we only need to diagonalize the Hamiltonian in the S = 1 subspace.In mean-field calculations we need to calculate the potentials on a cylindrical lattice. The number of mesh points in the z- and ρ- directions are nz = 24 and nρ = 12, respectively. In MDCCDFT, the triaxial deformation is included by the Fourier transformation given in equation (13 ). Choosing z as the principal axis, only the n = 0 Fourier component survives in the axially symmetric case, which simplifies the calculation, regardless of φ. A triangular deformed system lying on the x-z plane is not axially symmetric about the intrinsic z-axis. In such a system, n > 1 components in the Fourier transformation need to be taken into account. Due to the V4 symmetry restriction, we only need to calculate the components with even n in the Fourier transformation. The V4 symmetry can simplify the calculation into 1/4 of the full space with respect to φ and we choose the number of points for the φ integral as nφ = 12 for 0 < φ < π/2.
2.3. Angular momentum and parity projections
The AMP operator for a general wave function breaking both the axial symmetry and reflection symmetry can be written as
$\begin{eqnarray}{\hat{P}}_{{MK}}^{I}=\displaystyle \frac{2I+1}{8{\pi }^{2}}\int {\rm{d}}{\rm{\Omega }}{D}_{{MK}}^{I* }({\rm{\Omega }})\hat{R}({\rm{\Omega }}),\end{eqnarray}$
where ${D}_{{MK}}^{I}({\rm{\Omega }})={{\rm{e}}}^{-{\rm{i}}M\alpha }{d}_{{MK}}^{I}(\beta ){{\rm{e}}}^{-{\rm{i}}K\gamma }$ is the Wigner function with the Euler angle Ω ≡ {α, β, γ} and $\hat{R}({\rm{\Omega }})={{\rm{e}}}^{-{\rm{i}}\alpha {\hat{I}}_{z}}{{\rm{e}}}^{-{\rm{i}}\beta {\hat{I}}_{y}}{{\rm{e}}}^{-{\rm{i}}\gamma {\hat{I}}_{z}}$ is the rotational operator applied to the many-body wave function. I is the total angular momentum, M and K are angular momentum components along the z-axis in the laboratory frame and intrinsic frame, respectively. For nuclei with reflection asymmetric shapes, the PP should also be included to build states with definite parities [74]. The PP operator writes $\begin{eqnarray}{\hat{P}}^{\pi }=\displaystyle \frac{1}{2}(1+\pi \hat{P}),\end{eqnarray}$
where $\hat{P}$ is the spatial reflection operation and π = ± 1 is the parity.With ${\hat{P}}_{{MK}}^{I}$ and ${\hat{P}}^{\pi }$, we can project out states with good angular momentum I, M and parity π from an RHB mean-field wave function ∣Φ⟩,16 ) as the trial wave function and solving the variational problem against the variables ${f}_{K}^{I\pi }$, we end up with the Hill–Wheeler–-Griffin equation [74]
$\begin{eqnarray}\begin{array}{l}\left|{{\rm{\Psi }}}_{M}^{I\pi }\right\rangle =\sum _{K}{f}_{K}^{I\pi }{\hat{P}}^{\pi }{\hat{P}}_{{MK}}^{I}\left|{\rm{\Phi }}\right\rangle \\ =\sum _{K}{f}_{K}^{I\pi }\left|{IMK}\pi \right\rangle ,\end{array}\end{eqnarray}$
where ${f}_{K}^{I\pi }$ are some undetermined constants. Taking equation ( $\begin{eqnarray}\sum _{K^{\prime} }[{{ \mathcal H }}_{{KK}^{\prime} }^{I\pi }-{E}_{\alpha }^{I\pi }{{ \mathcal N }}_{{KK}^{\prime} }^{I\pi }]{f}_{K^{\prime} \alpha }^{I\pi }=0,\end{eqnarray}$
where the Hamiltonian kernel ${ \mathcal H }$ and the norm kernel ${ \mathcal N }$ are defined as $\begin{eqnarray}\begin{array}{l}{{ \mathcal H }}_{{KK}^{\prime} }^{I\pi }=\left\langle {\rm{\Phi }}\right|{{HP}}^{\pi }{P}_{{KK}^{\prime} }^{I}\left|{\rm{\Phi }}\right\rangle \\ =\displaystyle \frac{2I+1}{16{\pi }^{2}}\displaystyle \int \left[{D}_{{KK}^{\prime} }^{{I}^{* }}({\rm{\Omega }})+\pi {\left(-1\right)}^{I}{D}_{K-K^{\prime} }^{I}({\rm{\Omega }})\right]\\ \times \left\langle {\rm{\Phi }}\right|\hat{H}\hat{R}({\rm{\Omega }})\left|{\rm{\Phi }}\right\rangle {\rm{d}}{\rm{\Omega }},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{{ \mathcal N }}_{{KK}^{\prime} }^{I\pi }=\left\langle {\rm{\Phi }}\right|{P}^{\pi }{P}_{{KK}^{\prime} }^{I}\left|{\rm{\Phi }}\right\rangle \\ =\displaystyle \frac{2I+1}{16{\pi }^{2}}\displaystyle \int \left[{D}_{{KK}^{\prime} }^{{I}^{* }}({\rm{\Omega }})+\pi {\left(-1\right)}^{I}{D}_{K-K^{\prime} }^{I}({\rm{\Omega }})\right]\\ \times \left\langle {\rm{\Phi }}\right|\hat{R}({\rm{\Omega }})\left|{\rm{\Phi }}\right\rangle {\rm{d}}{\rm{\Omega }}.\end{array}\end{eqnarray}$
Note that the expressions of the kernels are derived by considering the V4 symmetry. By discretizing the group space, the problem turns to the calculation of the norm overlap n(Ω) and the Hamiltonian overlap h(Ω) as functions of the Euler angle Ω, $\begin{eqnarray}n({\rm{\Omega }})=\left\langle {\rm{\Phi }}\right|\hat{R}({\rm{\Omega }})\left|{\rm{\Phi }}\right\rangle ,\end{eqnarray}$
$\begin{eqnarray}h({\rm{\Omega }})=\left\langle {\rm{\Phi }}\right|\hat{H}\hat{R}({\rm{\Omega }})\left|{\rm{\Phi }}\right\rangle /n({\rm{\Omega }}).\end{eqnarray}$
We adopt the method introduced in [94] to calculate the overlaps in the ADHO basis using the generalized Wick's theorem.We follow the standard way of solving the Hill–Wheeler–Griffin equation (17 ), as discussed in [74]. The first step is to diagonalize the norm kernel
$\begin{eqnarray}\sum _{K^{\prime} }{{ \mathcal N }}_{{KK}^{\prime} }^{I\pi }{u}_{m}^{{IK}^{\prime} \pi }={n}_{m}^{I\pi }{u}_{m}^{{IK}\pi }.\end{eqnarray}$
The normalized vectors can be used to build the orthonormalized basis $\begin{eqnarray}\left|m\right\rangle =\displaystyle \frac{1}{\sqrt{{n}_{m}^{I\pi }}}\sum _{K}{u}_{m}^{{IK}\pi }{\hat{P}}^{\pi }{\hat{P}}_{{MK}}^{I}\left|{\rm{\Phi }}\right\rangle .\end{eqnarray}$
Eigenvectors with small eigenvalues ${n}_{m}^{I\pi }$ may induce numerical instabilities. To solve this problem, we introduce a cutoff χ = 10−3 and only consider eigenstates with eigenvalues ${n}_{m}^{I\pi }\gt \chi $. In the othornormalized basis ∣m⟩ the Hill–Wheeler–Griffin equation becomes a standard eigenvalue problem, $\begin{eqnarray}\sum _{m^{\prime} }\left\langle m| \hat{H}| m^{\prime} \right\rangle {g}_{m^{\prime} }^{I\pi }={E}^{I\pi }{g}_{m}^{I\pi }.\end{eqnarray}$
The eigen energies ${E}_{\alpha }^{I\pi }$ can be solved by matrix diagonalization and the wave function components ${f}_{K}^{I\pi }$ can be obtained from the corresponding eigenvectors ${g}_{m}^{I\pi }$, $\begin{eqnarray}{f}_{K}^{I\pi }=\sum _{m,{n}_{m}^{I\pi }\ne 0}\displaystyle \frac{{g}_{m}^{I\pi }}{\sqrt{{n}_{m}^{I\pi }}}{u}_{m}^{{IK}\pi }.\end{eqnarray}$
With the projected wave functions, the reduced transition probabilities B(Eλ) from the Ji state to the Jf state writes $\begin{eqnarray}\begin{array}{l}B(E\lambda ;{J}_{i}\to {J}_{f})=\displaystyle \frac{1}{2{J}_{i}+1}\\ \,\times \,\sum _{{M}_{i}{M}_{f}\mu }{\left|\left\langle {J}_{f},{M}_{f}\right|{\hat{Q}}_{\lambda \mu }\left|{J}_{i},{M}_{i}\right\rangle \right|}^{2}\\ =\displaystyle \frac{1}{2{J}_{i}+1}{\left|\left\langle {J}_{f}| \right|{\hat{Q}}_{\lambda }\left|| {J}_{i}\right\rangle \right|}^{2}.\end{array}\end{eqnarray}$
The reduced matrix element has the following expression $\begin{eqnarray}\left\langle {J}_{f}| \right|{\hat{Q}}_{\lambda }\left|| {J}_{i}\right\rangle =\sum _{{K}_{i}{K}_{f}}{f}_{{K}_{f}}^{{I}_{f}{\pi }_{f}}{f}_{{K}_{i}}^{{I}_{i}{\pi }_{i}}{\left[{Q}_{\lambda }\right]}_{{K}_{f}{K}_{i}}^{{J}_{f}{J}_{i}},\end{eqnarray}$
where $\begin{eqnarray}\begin{array}{l}{\left[{Q}_{\lambda }\right]}_{{K}_{f}{K}_{i}}^{{J}_{f}{J}_{i}}=\displaystyle \frac{(2{J}_{i}+1)(2{J}_{f}+1)}{8{\pi }^{2}}{\left(-1\right)}^{{J}_{f}-{K}_{f}}\\ \,\times \,\sum _{\mu M}\left(\begin{array}{ccc}{J}_{f} & \lambda & {J}_{i}\\ -{K}_{f} & \mu & M\end{array}\right){Q}_{\lambda \mu },\end{array}\end{eqnarray}$
in which $\begin{eqnarray}{Q}_{\lambda \mu }=\int {\rm{d}}{\rm{\Omega }}{D}_{{{MK}}_{i}}^{{J}_{i}* }({\rm{\Omega }})\left\langle {\rm{\Phi }}\right|{\hat{Q}}_{\lambda \mu }\hat{R}({\rm{\Omega }})\left|{\rm{\Phi }}\right\rangle .\end{eqnarray}$
3. Results and discussion
In this section, we apply the p-MDCRHB model to study the low-lying excitations of the nucleus 12C, which is a typical nucleus with possible α clustering structures. Among the proposed clustering states, the second 0+ state of 12C, namely the Hoyle state, is the most interesting because it plays a crucial role in the stellar nuclear synthesis [111]. The Hoyle state is a resonance consisting of three interacting α clusters lying above the three-α breakup threshold [112]. The energy and width of this resonance play an essential role in synthesizing three α's into 12C in the main sequence stars. The abundance of the carbon elements would be insufficient for supporting life in our universe if the Hoyle state has a different energy or does not exist at all. Furthermore, several recent studies of the ground state band indicate that the 3-α structure may already appear in the ground state of 12C [53, 113]. A 5− state was measured in the ground state band [54], providing evidence of a triangular D3h symmetry.
In the MDCRHB model and its extensions, the system respects the mirror symmetries about the x-z and y-z planes. One interesting case in which these models specialize in is the clustering structure in the light α-conjugate nuclei. Here, the α clusters can form various molecular-like structures and the possible shapes include rod, triangle, tetrahedron, octahedron, etc. Most of these shapes cannot be described by a conventional multipole expansion method that only includes either triaxial or octupole deformations. Further, the structure might be different for the ground state and the excited state, a mean-field description is too simplified and we must take into account the symmetry restorations with projection techniques.
In this section, we apply the p-MDCRHB model to study the ground state and low-lying excitations of a 12C nucleus with an intrinsic triangular shape. To generate such a configuration, we define a triangular deformation constraint in the mean-field calculations. In analogy with the usual multipole deformation constraint equation (10 ), we use the ansatz29 ) has the desired symmetry. If rotated around the y-axis, the operator S3 has a period 2π/3 and the maxima point to the directions along which the nucleus is stretched. The parameter λ3 is the counterpart of the multipole deformation parameters βλμ. When constrained to a large λ3, the operator S3 elongates the nucleus in the directions θ = 0, 2π/3, 4π/3 and compresses it in other directions. The resulting shape would be a regular triangle with three 4He at the vertices. In the MDCRHB model, the triangular constraint can be implemented in exactly the same way with the multipole moment constraint. In figure 1 we show the density profiles of 12C constrained to λ3 = 0.8. We present the contours in all three coordinate planes x = 0, y = 0, z = 0 and a 3D isosurface. The densities are consistent with what we expect from the structure of the triangular moment constraint. The energetically preferred directions are distributed evenly in the x-z plane. With much stronger constraining potential, individual α clusters can form and be squeezed out of the nucleus.
$\begin{eqnarray}{\lambda }_{3}=\displaystyle \frac{4\pi }{3{{NR}}^{3}}\langle {S}_{3}\rangle ,\end{eqnarray}$
where the ‘triangular moment operator' is defined as $\begin{eqnarray}{S}_{3}={\left({z}^{2}+{\rho }^{2}{\cos }^{2}\varphi \right)}^{3/2}\cos (3\theta ),\end{eqnarray}$
where ρ, φ and z are the cylindrical coordinates and $\theta $ is the polar angle on the x-z plane. We can verify that equation (
Figure 1. Density profiles of 12C constrained for λ3 = 0.8. The contour plots display the total density ρtotal in the x-z, y-z and x-y planes (top left, top right and bottom left, respectively). The 3D plot (bottom right) shows the isosurface of ρtotal = 0.08 fm−3. |
In mean-field calculations, the Dirac equation is solved by expanding the wave functions in an ADHO basis with a major shell truncation Nf. For Dirac spinors we use Nf shells to expand the large component and Nf + 1 shells to expand the small component to remove the spurious states. Firstly, we examine the convergence of the mean-field results against the basis truncation Nf. In figure 2 we present the PESs of 12C calculated with the triangular constraint using different values of Nf. We found that all energy minima locate at λ3 = 0 corresponding to a spherical shape. There is no explicit triangular correlation in the mean-field level. As Nf increases, differences between adjacent curves become smaller, signifying a clear convergence pattern. At the ground state λ3 = 0, the calculated binding energy already converges at Nf = 8, with the truncation error smaller than 0.1 MeV estimated by comparing to the result at Nf = 12. For a larger λ3, the mean-field potential assumes a triangular shape, differing from that of the harmonic oscillator basis. In this case we observe a slower convergence. For example, the mean-field energies for Nf = 10 and Nf = 12 differ by about 1 MeV for λ3 ≈ 3.0. Nevertheless, the qualitative behaviour of the PES is not changed when Nf is large enough. Concerning the big computational cost of the 3D AMP, in this illustrative work we use Nf = 8 in all of the following calculations. For more realistic calculations we need to either employ a triangular harmonic oscillator basis or increase Nf until all truncation errors are smaller than the desired precision.
Figure 2. The energy of 12C as functions of the triangular shape parameter λ3 calculated with different basis truncations in the MDCRHB model. The squares, circles, up triangles and down triangles denote the results obtained with major shell truncation Nf = 6, 8, 10 and 12, respectively. The center of mass correction to the binding energy is not included. |
For angular momentum and parity projections, the Hamiltonian kernel equation (18a ) and the norm kernel equation (18b ) are calculated by integrating over the SO(3) group space parameterized by the Euler angles (α, β, γ). These integrals are performed numerically by the Gauss–Legendre integral with numbers of mesh points (nα, nβ, nγ) in the α, β and γ directions, respectively. As most of the computational cost comes from the calculation of the Hamiltonian overlap and norm overlap for each Euler angle, it is desirable to reduce the total number of integral mesh points for a given precision. Here, we make a convergence check against the number of mesh points in the group space. In figure 3, we show the projected energies as functions of nα, nβ and nγ. On the left panel, we set nβ = 8 and study the projected energies for 6 ≤ nα = nγ ≤ 14, where EJπ denotes the energy of the lowest $\left|{J}^{\pi }\right\rangle $ state with the corresponding Euler angle mesh size, and ${E}_{0}^{J\pi }$ represents the energy with the largest mesh size. We see that the influence of nα and nγ on the projected energies for J < 6 is less than 10 keV, which is negligible compared with the excitation energies of the order of MeV. In contrast, the right panel shows that the number nβ has a rather large influence on the projected energies. All projected energies up to 6+ achieve relative precisions less than 10−6 when nβ ≥ 10. In this work, we choose the mesh sizes nα = nγ = 12 and nβ = 14, with which the numerical integral in the group space can be evaluated rather precisely for J ≤ 6.
Figure 3. The projected energies of 12C calculated with different mesh sizes for the SO(3) group integral. For each curve, the energies are relative to the corresponding converged values E0J calculated at large mesh sizes. All calculations are performed with a mean-field wave function constrained to λ3 = 1.2. Here, nα, nβ and nγ are the number of mesh points for integrating the Euler angles α, β and γ, respectively. Squares, circles, up triangles, down triangles, diamonds and left triangles denote results projected to Jπ = 2+, 4+, 6+, 3−, 4−, 5−, respectively. |
The constraint equation (29 ) produces a triangular shape with the 3-fold axis pointing along the y-axis and one of the 2-fold axes pointing along the z-axis. For AMP we rotate the mean-field wave function to Euler angles (α, β, γ) and calculate its overlap with the original wave function. These norm overlaps should exhibit periodicity if the rotations are about the symmetry axes of the intrinsic shape, which can be examined numerically. In figure 4 we show the norm overlaps for different Euler angles. All results are calculated with a mean-field wave function constrained to λ3 = 1.2. In the left panels we set β = γ = 0 and rotate the nucleus around the z-axis by increasing α from 0° to 360°. As the z-axis is a 2-fold axis of the triangle, we found a period of 180° for both the norm overlap and Hamiltonian overlap. In the right panels we set α = γ = 0; then variation of β corresponds to a rotation around the y-axis. In this case we observe a period of 120° for both overlaps, which is consistent with the fact that the y-axis is a 3-fold axis of the triangle. In all subplots, the overlaps return to their original values at α, β, γ = 0 after rotating by the corresponding period. This is the direct consequence of the invariance of the nuclear shape under these symmetry operations.
Figure 4. The norm and Hamiltonian overlaps of 12C as functions of the Euler angles α and β, calculated with a mean-field wave function constrained to λ3 = 1.2. β = γ = 0 on the left panels and α = γ = 0 on the right panels. |
The spectrum of a triangular rotor can be mostly determined using the group theory analysis. The triangle shape corresponds to a point group D3h, which has 12 elements and six irreducible representations (irreps) ${A}_{1}^{{\prime} }$, ${A}_{2}^{{\prime} }$, ${E}^{{\prime} }$, ${A}_{1}^{{\prime\prime} }$, ${A}_{2}^{{\prime\prime} }$ and E″. As D3h is a subgroup of the O(3) group, we can decompose the irreps of the O(3) group into irreps of the D3h group [114–116]. The O(3) group has infinitely many irreps corresponding to different angular momentum I and parity π. In table 1 we present the decomposition of the O(3) irreps into the D3h irreps. In each column we show the number of each D3h irrep appearing in the decomposition. For example, the identity representation corresponds to 0+ for the O(3) group and is marked by ${A}_{1}^{{\prime} }$ for the D3h group, while the irrep of O(3) with Iπ = 1+ can be written as a direct sum of ${A}_{2}^{{\prime} }$ and E″ representations.
Table 1. The decomposition of the irreps of the O(3) group with quantum numbers Iπ into irreps of the D3h group. In each column the numbers are the degeneracies. |
| Iπ | 0+ | 1+ | 2+ | 3+ | 4+ | 5+ | 6+ | 7+ | 0− | 1− | 2− | 3− | 4− | 5− | 6− | 7− |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| $A{{\prime} }_{1}$ | 1 | 0 | 1 | 0 | 1 | 0 | 2 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 |
| $A{{\prime} }_{2}$ | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 |
| $E^{\prime} $ | 0 | 0 | 1 | 1 | 2 | 2 | 2 | 2 | 0 | 1 | 1 | 1 | 1 | 2 | 2 | 3 |
| A″1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 2 | 1 |
| A″2 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 2 |
| E″ | 0 | 1 | 1 | 1 | 1 | 2 | 2 | 3 | 0 | 0 | 1 | 1 | 2 | 2 | 2 | 2 |
For any state of the triangular rotor with quantum numbers Iπ, the wave function transforms under the corresponding irrep of the O(3) group. If we rotate the rotor with an element of the D3h group in the intrinsic frame, nothing is actually changed; thus the wave functions should be scalars for these rotations. In the language of the group theory, the allowed states all belong to the ${A}_{1}^{{\prime} }$ irrep of the intrinsic D3h group. If an angular momentum and parity combination Iπ does not contain any ${A}_{1}^{{\prime} }$ irrep, it should not appear in the spectrum. In other words, the Iπ states with non-zero entries in the first row of table 1 are all states allowed by the symmetry analysis. In summary, the lowest energy levels of a triangular rigid rotor should be Iπ = 0+, 2+, 3−, 4+, 4−, 5−, 6+, 6−, … Other angular momentum and parity combinations should be consequences of the intrinsic excitations.
On the one hand, the group theory analysis is strict and can be generalized to any other intrinsic shapes. For example, we can predict the low-lying states of a tetrahedral rotor and compare the results with that of the 16O nucleus, which was predicted to exhibit a four-α cluster structure. On the other hand, there are also less abstract ways of understanding the essentiality of the nuclear rotational spectrum. For a triangular rotor, a stable rotation can only occur around its symmetry axes. The wave function for angular momentum I is $\psi (\varphi )=\exp ({\rm{i}}I\varphi )$, where φ is the rotational angle. For an n-fold axis we must have ψ(φ + 2π/n) = ψ(φ); thus, among all integer values of I, only integer multiples of n are allowed. For example, the rotations about a 2-fold axis generate angular momenta I = 0, 2, 4, 6, …, while the rotations about a 3-fold axis give the angular momenta I = 0, 3, 6, … Further, the parity of the rotational wave function can be determined by $\pi ={\left(-1\right)}^{I}$. The total rotational energy of a triangular rotor is the sum of that around the 2-fold and 3-fold axes. If there is no rotation around the 3-fold axis, we have the spectrum Iπ = 0+, 2+, 4+, 6+, … If there is one unit of quantized angular momentum around the third axis, we have the states Iπ = 3−, 4−, 5−, 6−, … Here, the states 4− and 5− can be understood as a rotation around the axis in between the 2-fold axis and the 3-fold axis, which is already a complex motion in the classical level.
Next, we turn to the projected PES. We first generate mean-field configurations with λ3 varying from 0.0 to 5.0. The results describe the disintegration of 12C into three α clusters. In figure 5 we show the projected energies as well as mean-field energies for all λ3 in this interval with a step size Δλ3 = 0.2. Here, we only show the lowest energy levels with Iπ = 0+, 2+, 3−, 4± and 5−. The mean-field PES gives a ground state at the spherical shape λ3 = 0. The energy increases quickly as three α's are separated. For λ3 = 1.0 the deformation energy is more than 5 MeV, preventing the dissociation of the nucleus. By projecting out the spurious rotational energies, we see that the PES for the 0+ state is quite different from the mean-field PES. The 0+ PES is lower than the mean-field PES and becomes very flat for λ3 < 2.0. We find a minimum at λ3 = 1.2, which is very shallow in that the energy at λ3 = 0 is only about 1 MeV higher. In the insets of figure 5 we show the density profiles from the mean-field calculations at λ3 = 1.2 and λ3 = 3.8. The first inset corresponds to the ground state identified with the 0+ PES. Although, at this point there is no emerging clustering structure, we can see a clear triangular shape from the total densities. At λ3 = 3.8 the deformation is so strong that the clustering structure is very pronounced. Theoretically, for extremely large λ3, the 12C nucleus would break up into three α particles and the energy should converge to three times the 4He binding energy. However, for such light nuclei we might need a very large center of mass correction and the mean-field model might not be a good approximation. Even though there are still a lot of difficulties, the extremely large λ3 deformation still deserves some further investigations. This three-α disintegration can be viewed as an inverse reaction of the famous 3α → 12C reaction and may be related to the Hoyle state. Whether the mean-field model is adequate for describing such a process and what is lacking is thus an interesting topic. Although in recent years ab initio calculations have been more popular for studying clustering structures, mean-field models have clear physical pictures and have proven successful on this topic. The p-MDCRHB model we presented here can describe various α cluster configurations and can serve as a starting point for beyond-mean-field calculations in this direction.
Figure 5. The mean-field and the projected potential energies of 12C versus the triangular deformation λ3. The squares, circles, down triangles, up triangles, diamonds and left triangles represent Iπ = 0+, 2+, 3−, 4+, 4− and 5−, respectively. The mean-field energies are displayed as a black solid line. The insets show the mean-field density profiles at λ3 = 1.2 and 3.8. |
The PESs for higher angular momenta tell us how the α clusters move about the rotational axes. At λ3 = 0 the system is spherical and there is no rotation. For larger λ3, the 2+ and 4+ states describe rotations around the 2-fold axes, and the corresponding PESs have mimina at relatively small λ3. The 3− and 4− states contain rotations around the 3-fold axis; thus, the centrifugal force prefers a larger λ3, which also means a larger moment of inertia.
The energy spectra of 12C at large λ3 can be understood by comparing them to the spectra of a triangular rigid rotor discussed above. By taking such a rotor with equal masses on the three vertices, we can compute the moment of inertia ${ \mathcal M }$ along every principle axis. For the orientation shown in figure 1, the principle axes are simply the three coordinate axes x, y and z. It is straightforward to prove that $2{{ \mathcal M }}_{x}$ = $2{{ \mathcal M }}_{z}$ = ${{ \mathcal M }}_{y}$ = ${ \mathcal M }$. By only considering rotations about the center of mass, we can write the Hamiltonian as30 ),
$\begin{eqnarray}\hat{H}=\displaystyle \frac{{\hat{I}}_{x}^{2}}{{ \mathcal M }}+\displaystyle \frac{{\hat{I}}_{y}^{2}}{2{ \mathcal M }}+\displaystyle \frac{{\hat{I}}_{z}^{2}}{{ \mathcal M }}=\displaystyle \frac{{\hat{I}}^{2}}{{ \mathcal M }}-\displaystyle \frac{{\hat{I}}_{y}^{2}}{2{ \mathcal M }},\end{eqnarray}$
where ${\hat{I}}_{x}$, ${\hat{I}}_{y}$ and ${\hat{I}}_{z}$ are the angular momentum components in the intrinsic frame and ${\hat{I}}^{2}$ is the total angular momentum. As $\hat{H}$, ${\hat{I}}^{2}$ and ${\hat{I}}_{y}^{2}$ commute with each other, we can take the total energy E, total angular momentum I and the angular momentum component K as a group of good quantum numbers specifying the eigenstates. Here, I(I + 1) is the eigenvalue of ${\hat{I}}^{2}$ and K is the eigenvalue of ${\hat{I}}_{y}$. Rotation around the intrinsic y-axis by 2π/3 takes the nucleus back to its original position, while it produces a phase $\exp ({\rm{i}}2\pi K/3)$. This phase should be unit one for any physically allowed state; thus we have K = 3n, where n is an integer. For a state ∣IKα⟩ with α quantum numbers other than I and K, we can calculate the energy according to equation ( $\begin{eqnarray}\begin{array}{rcl}E & = & \langle {IK}\alpha | \hat{H}| {IK}\alpha \rangle \\ & = & \displaystyle \frac{1}{{ \mathcal M }}[I(I+1)-{K}^{2}/2].\end{array}\end{eqnarray}$
The spatial reflection P combined with a rotation around the 3-fold axis by 180° resets the system. The operator P produces a factor of parity π and the rotation gives an additional phase of ${\left(-1\right)}^{K}$; thus we have $\pi ={\left(-1\right)}^{K}$ for all physically allowed states.The energy spectrum of equation (31 ) can be classified according to the quantum number K. For K = 0 we have the ground state band 0+, 2+, 4+, ..., with positive parity. For K = 3 we have a negative parity band 3−, 4−, 5−, ... The excitation energies relative to the 0+ state are
$\begin{eqnarray}\begin{array}{l}{E}_{x}({4}^{+}):{E}_{x}({4}^{-}):{E}_{x}({3}^{-}):{E}_{x}({2}^{+})\\ \ \ =3.333:2.583:1.25:1,\end{array}\end{eqnarray}$
up to an overall multiplicative factor. One important feature of this spectrum is that the 4+ energy is higher than the 4− energy. The reason for this is that the 4+ state is a pure rotation around the 2-fold axes, while the 4− state contains rotations around the 3-fold axis. The moment of inertia for the latter rotation is relatively larger. However, in figure 5 we see that the 4+ state is much lower than 3− and 4− for small λ3; only for extremely large λ3 can the projected energies be explained by the rigid rotor model. The reason is twofold. First, the relation ${{ \mathcal M }}_{y}$ = $2{{ \mathcal M }}_{z}$ used to derive the rigid rotor spectrum is only approximately satisfied for large deformations where the masses of α clusters are concentrated at the vertices. For small λ3 we have a triangular density distribution without clustering and the moments of inertia along different axis are not so different. Second, the angular momentum and parity projections capture the full quantum mechanical elements. The impact of the intrinsic states is reflected by the calculated norm kernel and Hamiltonian kernel. Only for large λ3 can we omit the intrinsic structures and use the Gaussian overlap approximation to treat the α clusters as point particles.In table 2 and table 3 we compare the low-lying spectrum and reduced transition probabilities B(E2) and B(E3) of 12C with the experiments. The results are calculated by projecting a mean-field state constrained to λ3 = 1.2, which corresponds to the minimum of the 0+ PES in Fig 5. The experimental values are well reproduced at this level of approximations. These results can be improved from various aspects. First, as the PESs are soft at small λ3, in figure 5 we see that the λ3 value minimizing the projected energy has a strong dependence on Iπ, which should be considered when calculating the energy spectrum. Second, the soft PESs suggest we introduce the GCM to mix different deformations. Lastly, other configurations like the rod shape might also play an important role in the low-lying excitations and should be considered when building the low-lying states with good quantum numbers.
Table 2. Reduced E2 transition probabilities from states ${I}_{i}^{\pi }$ to states ${I}_{f}^{\pi }$ and the corresponding excitation energies of 12C. The calculation is made for λ3 = 1.2, corresponding to the energy minimum of the 0+ PES shown in figure 5. The experimental B(E2) value is taken from [117]. |
| ${I}_{i}^{\pi }$ | ${E}_{{xi}}^{\mathrm{Th}}$ (MeV) | ${E}_{{xi}}^{\mathrm{Expt}}$(MeV) | ${I}_{f}^{\pi }$ | ${E}_{{xf}}^{\mathrm{Th}}$ (MeV) | $B{\left(E2\right)}^{\mathrm{Th}}$ (e2 fm4) | $B{\left(E2\right)}^{\mathrm{Expt}}$ (e2 fm4) |
|---|---|---|---|---|---|---|
| 2+ | 1.951 | 4.440 | 0+ | 0.000 | 9.72 | 7.6 ± 0.4 [117] |
| 4+ | 9.015 | 14.079 | 2+ | 1.951 | 16.28 | |
| 4− | 15.253 | 3− | 12.864 | 20.95 | ||
| 5− | 19.915 | 3− | 12.864 | 6.78 | ||
| 4− | 15.253 | 18.90 | ||||
| ${6}_{1}^{+}$ | 26.734 | 4+ | 9.015 | 26.04 | ||
| 6− | 25.131 | 4− | 15.253 | 13.38 | ||
| 5− | 19.915 | 15.01 |
Table 3. Reduced E3 transition probabilities from states ${I}_{i}^{\pi }$ to states ${I}_{f}^{\pi }$ and the corresponding excitation energies of 12C. The calculation is made for λ3 = 1.2, corresponding to the energy minimum of the 0+ PES shown in figure 5. The experimental B(E3) value is taken from [117]. |
| ${I}_{i}^{\pi }$ | ${E}_{{xi}}^{{\rm{Th}}}$(MeV) | ${E}_{{xi}}^{\mathrm{Expt}}$(MeV) | ${I}_{f}^{\pi }$ | ${E}_{{xf}}^{{\rm{Th}}}$(MeV) | $B{\left(E3\right)}^{\mathrm{Th}}$ (e3 fm6) | $B{\left(E3\right)}^{\mathrm{Expt}}$ (e3 fm6) |
|---|---|---|---|---|---|---|
| 3− | 12.864 | 9.641 | 0+ | 0.000 | 145.81 | 103 ± 17 [117] |
| 2+ | 1.951 | 260.25 | ||||
| 4+ | 9.015 | 53.53 | ||||
| 4− | 15.253 | 13.35 | 2+ | 1.951 | 268.69 | |
| 4+ | 9.015 | 198.55 | ||||
| 5− | 19.915 | 22.4(2) | 2+ | 1.951 | 96.28 | |
| 4+ | 9.015 | 318.35 | ||||
| ${6}_{1}^{+}$ | 26.734 | 4− | 15.253 | 33.95 | ||
| 5− | 19.915 | 129.86 | ||||
| 6− | 25.131 | 250.19 | ||||
| ${6}_{2}^{+}$ | 27.234 | 3− | 12.864 | 480.28 | ||
| 4− | 15.253 | 219.26 | ||||
| 5− | 19.915 | 90.68 | ||||
| 6− | 25.131 | 30.99 | ||||
| 6− | 25.131 | 4+ | 9.015 | 274.90 |
In [52], the α arrangement in 12C is assumed to be an equilateral triangle. In [118], the author performed calculations using the orthogonality condition model and the shallow potential model, and obtained different configurations for each model. Due to insufficient experimental data, whether the intrinsic configuration of the ground state of 12C is a regular triangle remains an open question [118]. It is necessary to explore this point by performing multidimensional constrained calculations that include all important shape configurations. In figure 6 we show the projected 0+ PES of 12C with respect to the deformation parameters β20 and λ3. Here, we constrain λ3 from 0.0 to 3.0 with a step size Δλ3 = 0.2, and constrain β20 from −0.2 to 1.0 with a step size Δβ20 = 0.1. The projected energies are much lower than the mean-field energies for non-zero deformations. The rotational energy correction, defined as the difference between the 0+ energy and the mean-field energy, plays a key role in determining the ground state configurations. On the 0+ PES, we found two energy minima. The lowest one locates at λ3 = 0.8 and β20 = 0.6, corresponding to an isosceles acute triangular shape, while the shape isomer at λ3 = 1.2 and β20 = − 0.1 corresponds to an isosceles obtuse triangular shape.
Figure 6. The 2D PES of 12C projected to the 0+ state. The λ3 and β20 constraints are imposed with 0 < λ3 < 3.0, Δλ3 = 0.2 and − 0.2 < β20 < 1.0; Δβ20 = 0.1. |
4. Summary
We have developed a p-MDCRHB model. In this model, the shape of the system is only restricted as mirror reflection symmetric about the x-z and y-z planes. In the mean-field calculation, the RHB equation is solved in the ADHO basis. The parity and angular momentum projections are performed to restore the intrinsic symmetries.
This model is well suited for exploring exotic shapes like triangles and tetrahedrons found in both the nuclear ground state and low-lying excited states. One important scenario where these shapes play essential roles is nuclear clustering. The most important stellar nuclear reactions, such as the production of carbon and oxygen, depend heavily on the clustering structures of the reactant. In this work we benchmark the p-MDCRHB model by studying the clustering states in 12C. We define a triangular moment constraint and study the corresponding triangular deformation λ3. This constraint generates a clustering state with three α clusters arranged on the vertices of a regular triangle. We found that with the growth of λ3, the spectra approach that of a regular triangular rigid rotor. This was explained by symmetry analysis and the rigid rotor Hamiltonian. We also performed 2D constrained calculations of 12C. Two minima are found on the projected 0+ PES. The softness of the PESs indicates that the configuration mixing GCM calculation is needed for a deeper understanding of the structure of 12C and the related nuclear reactions occurring in the stellar environment.
Acknowledgments
Helpful discussions with Xiang-Xiang Sun and Yu-Ting Rong are gratefully acknowledged. We thank Xiang-Xiang Sun and Shan-Gui Zhou for reading the manuscript and valuable suggestions. This work has been supported by NSAF (Grant No. U1930403), the National Key R&D Program of China (Grant No. 2018YFA0404402), the National Natural Science Foundation of China (Grant No. 11 525 524, No. 12 070 131 001, No. 12 047 503, and No. 11 961 141 004), the Key Research Program of Frontier Sciences of Chinese Academy of Sciences (Grant No. QYZDB-SSWSYS013) and the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB34010000 and No. XDPB15). The results described in this paper are obtained on the High-performance Computing Cluster of ITP-CAS and the ScGrid of the Supercomputing Center, the Computer Network Information Center of the Chinese Academy of Sciences.
Appendix: Calculation of the Rotation Matrix
Since the rotation matrix is naturally written as the Wigner D-function with spherical harmonic oscillators, the calculation of the rotation matrix can be simplified by expanding the ADHO basis with the spherical harmonic oscillator (SHO) basis.
An ADHO basis $\left|\alpha \right\rangle $ with quantum numbers {nz, nρ, ml, ms} reads
$\begin{eqnarray}\left|\alpha \right\rangle =\left|{n}_{z}{n}_{\rho }{m}_{l}{m}_{s}\right\rangle ={C}_{L}{\phi }_{{n}_{z}}(z){R}_{{m}_{l}}^{{n}_{\rho }}(\rho ,\varphi ){\chi }_{{m}_{s}},\end{eqnarray}$
with $\begin{eqnarray}{\phi }_{{n}_{z}}(z)=\displaystyle \frac{1}{\sqrt{{b}_{z}}}\displaystyle \frac{1}{{\pi }^{1/4}\sqrt{{2}^{{n}_{z}}{n}_{z}!}}{H}_{{n}_{z}}\left(\displaystyle \frac{z}{{b}_{z}}\right){{\rm{e}}}^{-\tfrac{{z}^{2}}{2{b}_{z}}},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{R}_{{m}_{l}}^{{n}_{\rho }}(\rho ,\varphi )=\frac{1}{{b}_{\rho }}\sqrt{\frac{2{n}_{\rho }!}{({n}_{\rho }+| m| )!}}\frac{1}{\sqrt{2\pi }}{{\rm{e}}}^{-{\rho }^{2}/(2{b}_{\rho }^{2})}\\ {\,}^{\times }{L}_{{n}_{\rho }}^{| m| }\left(\frac{{\rho }^{2}}{{b}_{\rho }^{2}}\right){{\rm{e}}}^{{\rm{i}}m\varphi }.\end{array}\end{eqnarray}$
Here, bz and bρ are the harmonic oscillator lengths along the z and ρ directions, which are defined as ${b}_{z}=1/\sqrt{M{\omega }_{z}}$ and ${b}_{\rho }=1/\sqrt{M{\omega }_{\rho }}$.The SHO basis function reads34 ) and (36 ), one can calculate the transformation matrix element, which has also been given in [119]
$\begin{eqnarray}\left|{n}_{r}{ljm}\right\rangle =\sum _{{m}_{s}}{C}_{{m}_{l}{m}_{s}{jm}}^{l\tfrac{1}{2}}{\phi }_{{{lm}}_{l}}^{{n}_{r}}(r,{\rm{\Omega }}){\chi }_{{m}_{s}},\end{eqnarray}$
where ${C}_{{m}_{l}{m}_{s}{jm}}^{l\tfrac{1}{2}}$ is the Clebsch–Gordan coefficient and $\begin{eqnarray}\begin{array}{l}{\phi }_{{{lm}}_{l}}^{{n}_{r}}(r,{\rm{\Omega }})=\displaystyle \frac{1}{{b}_{r}^{3/2}}\sqrt{\displaystyle \frac{2{n}_{r}!}{{\rm{\Gamma }}({n}_{r}+l+3/2)}}{\left(\displaystyle \frac{r}{{b}_{r}}\right)}^{l}\\ \times \,\exp \left(-\displaystyle \frac{{r}^{2}}{2{b}_{r}^{2}}\right){L}_{{n}_{r}}^{l+1/2}\left(\displaystyle \frac{{r}^{2}}{{b}_{r}^{2}}\right){Y}_{{lm}}({\rm{\Omega }}).\end{array}\end{eqnarray}$
The transformation matrix can be written as $\begin{eqnarray}\begin{array}{l}\left\langle {n}_{z}{n}_{\rho }{m}_{l}{m}_{s}| {n}_{r}{{ljm}}_{j}\right\rangle \\ ={C}_{L}^{* }{C}_{{m}_{l}{m}_{s}{{jm}}_{j}}^{l\tfrac{1}{2}}{\delta }_{{m}_{l}+{m}_{s},{m}_{j}}\left\langle {\phi }_{{n}_{z}}{R}_{{m}_{l}}^{{n}_{\rho }}| {\phi }_{{{lm}}_{l}}^{{n}_{r}}\right\rangle \\ ={C}_{L}^{* }{C}_{{m}_{l}{m}_{s}{{jm}}_{j}}^{l\tfrac{1}{2}}{\delta }_{{m}_{l}+{m}_{s},{m}_{j}}\left\langle {n}_{z}{n}_{\rho }{m}_{l}| {n}_{r}{{lm}}_{l}\right\rangle .\end{array}\end{eqnarray}$
From equations ( $\begin{eqnarray}\begin{array}{l}\left\langle {n}_{r}{lm}| {n}_{z}{n}_{\rho }m\right\rangle \\ =\sum _{\lambda }\displaystyle \frac{{n}_{\rho }!(2l+1){2}^{{n}_{\rho }-{n}_{r}-\lambda }{\left(-1\right)}^{\lambda }({n}_{r}+2\lambda )!}{\lambda !({n}_{\rho }-\lambda )!(\lambda -{n}_{\rho }+{n}_{r})!(2{n}_{r}+2\lambda -2{n}_{\rho }+2l+1)!!}.\end{array}\end{eqnarray}$
Thus, one obtains the relation between ADHO and SHO bases, $\begin{eqnarray}\left|\alpha \right\rangle =\sum _{{n}_{r}{ljm}}{C}_{{n}_{r}{ljm}}^{\alpha }\left|{n}_{r}{ljm}\right\rangle .\end{eqnarray}$
The time-reversal states of the SHO basis are obtained by applying the time-reversal operator $\hat{T}$
$\begin{eqnarray}\hat{T}\left|{n}_{r}{{ljm}}_{j}\right\rangle ={\left(-1\right)}^{l+{m}_{j}+j}\left|{n}_{r}{lj}{\bar{m}}_{j}\right\rangle .\end{eqnarray}$
The rotation matrix in the ADHO basis is calculated as41 ), while the other half can be calculated directly from the relation (42 ).
$\begin{eqnarray}\left\langle \beta \right|\hat{R}({\rm{\Omega }})\left|\alpha \right\rangle =\sum _{{n}_{r}{lj}}\sum _{{mm}^{\prime} }{C}_{{n}_{r}{ljm}}^{\beta * }{D}_{{mm}^{\prime} }^{j}({\rm{\Omega }}){C}_{{n}_{r}{ljm}^{\prime} }^{\alpha },\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\left\langle \bar{\beta }\right|\hat{R}({\rm{\Omega }})\left|\alpha \right\rangle \\ =\displaystyle \sum _{{n}_{r}{lj}}\displaystyle \sum _{{mm}^{\prime} }{\left(-1\right)}^{l+m-j}{C}_{{n}_{r}{ljm}}^{\beta }{D}_{{mm}^{\prime} }^{j}({\rm{\Omega }}){C}_{{n}_{r}{ljm}^{\prime} }^{\alpha }\end{array}\end{eqnarray}$
By applying the properties of the Wigner D-function, one finds the relation $\begin{eqnarray}{R}_{\bar{\beta }\bar{\alpha }}={R}_{\beta \alpha }^{* },\quad {R}_{\beta \bar{\alpha }}=-{R}_{\bar{\beta }\alpha }^{* },\quad \mathrm{for}\,\alpha ,\beta \gt 0,\end{eqnarray}$
that means only half of the R matrix elements are needed in equation (