Exact solution of U(5)–O(6) transitional description in interacting boson model with two-particle and two-hole configuration mixing

The interacting boson model (IBM) has proven to be extremely successful in the description of both the collective valence shell [1] and multi-particle-hole [2–4] excitations in nuclei. Most noticeably, the IBM Hamiltonian without configuration mixing can be solved analytically in U(5) (vibrational), O(6) (γ-unstable), and SU(3) (rotational) limits [1], as well as in the U(5)–O(6) transitional case [5]. In contrast, configuration mixing due to multi-particle-hole excitations was considered to gain an understanding of the shape coexistence phenomena by assuming different symmetry limits of the IBM for different configurations [6–12], which has proven to be successful in describing intruder states and shape coexistence phenomena in near closed shell nuclei, typically those around proton numbers Z~50 and Z~82 [2–4]. Recently, the intruder configuration mixing schemes with $ 2n $-particle and $ 2n $-hole configurations from $ n = 0 $ up to $ n\rightarrow\infty $ in the U(5) (vibrational) and the O(6) (γ-unstable) limits of the IBM-I were proposed [13, 14], whose simple Hamiltonians, suitable to describe the intruder and normal configuration mixing, prove to be exactly solvable based on the SU(1,1) coherent states.

The configuration mixing schemes in the IBM [2–4] can be considered in both IBM-II and IBM-I with no distinction between neutron-type and proton-type bosons, as shown in [7–10, 13–16]. In this study, we demonstrate that the U(5)↔O(6) transitional Hamiltonian of the IBM-I with two-particle and two-hole configuration mixing is also exactly solvable based on the Bethe ansatz approach. The results of $ N = 2 $ and $ N = 4 $ cases are considered as examples to demonstrate the feature of the solution. To apply this theory, the model is employed to fit some low-lying level energies and B(E2) ratios of $ ^{108} {\rm{Cd}}$.

The Hamiltonian of the U(5)–O(6) transitional description in the IBM-I with two-particle and two-hole configuration mixing is expressed as [2–4]

where $ P_{N} $ and P are the projection operators, where $ P_{N} $ projects to the N-boson subspace, whereas P projects to the subspace with N and $ N+2 $ bosons,

For $ i = 1 $ and 2, these are the U(5)–O(6) transitional Hamiltonians [5], and

is the two-configuration mixing term. Here, $ S^{+} = S^{+}_{s}-S^{+}_{d} $ with $ S_{s}^{+} = {1\over{2}}s^{\dagger 2} $ and $ S^{+}_{d} = {1\over{2}}d^{+}\cdot d^{+} = {1\over{2}}\sum_{\mu}(-1)^{\mu}d^{\dagger }_{\mu}d^{+}_{-\mu} $, where $ s^{\dagger} $ (s) and $ d^{\dagger}_{\mu} $ ($ d_{\mu} $) are the creation (annihilation) operators of s- and d-bosons, respectively, $ S^{-}_{\rho} = \left(S^{+}_{\rho}\right)^{\dagger } $ for $ \rho = s $ or d, $ S^{0}_{s} = {1\over{2}}(\hat{n}_{s}+{1\over{2}}) $ and $ S^{0}_{d} = {1\over{2}}(\hat{n}_{d}+{5\over{2}}) $, respectively, with $ \hat{n}_{s} = s^{\dagger }s $ and $ \hat{n}_{d} = \sum_{\mu}d^{\dagger }_{\mu}d_{\mu} $, and $ \alpha^{(i)}_{\rho} $, $ g^{(i)} $, $ g_{s} $ , and $ g_{d} $ are real parameters. The U(5)$ \leftrightarrow $O(6) Hamiltonian (2), which is equivalent to the consistent-Q formulism of IBM but different from that used in Ref. [5], is in the U(5) (vibrational) phase when $ \alpha_{\rho}^{(i)}\neq 0 $ and $ g^{(i)} = 0 $ or in the O(6) (γ-unstable) phase when $ \alpha_{\rho}^{(i)} = 0 $ and $ g^{(i)}<0 $. Because $ 2(\alpha_{d}^{(i)}-\alpha_{s}^{(i)}) $ is the energy gap of d- and s-bosons, $ \alpha_{s}^{(1)} $ is considered to be zero. In a previous configuration mixing study, $ \Delta = 2\alpha_{d}^{(2)}-2\alpha_{d}^{(1)} = 2\alpha_{s}^{(2)}-2\alpha_{s}^{(1)} $ is considered according to the energy of the lowest intruder state to reduce the number of parameters. The shape phase within the N-boson and $ N+2 $-boson configuration controlled by the two sets of parameters $ \{\alpha_{\rho}^{(i)},\,g^{(i)}\} $ ($ i = 1,2 $) can be different in terms of their description of the shape (phase) coexistence.

The Hamiltonian (1) can be diagonalized in the $ N\oplus (N+2) $-boson subspace, where the complete basis vectors in each configuration can be considered as those of U(6)$\supset $U(5)$\supset $O(5)$\supset $O(3) with $ \vert N\,n_{d}\, \nu_{d}\,\eta\, L\, M\rangle $, where N is the total number of bosons, $ n_{d} $ is the number of d-bosons, $ \nu_{d} $ is the d-boson seniority number labeling the irrep of $ O(5) $, L is the angular momentum quantum number, M is the quantum number of the third component of the angular momentum, and η is an additional quantum number required to distinguish different states with the same L. Moreover, the two sets of operators $ \{{S}^{\pm}_{\rho},\; {S}^{0}_{\rho}\} $ ($ \rho = s,d $) that are two copies of the SU(1,1) algebra, satisfy the commutation relations

Equivalently, for a given N, $ n_{d} $, $ \nu_{d} $, η, L, and M, the orthonormalized basis vectors $ \vert N\,n_{d}\, \nu_{d}\,\eta\, L\, M\rangle $ can also be expressed as those of $ {{SU}}_{d} $(1,1)$ \otimes {{SU}} _{s} $(1,1) with

where $ n_{d} = 2\xi+\nu_{d} $ and $ \xi = 0,1,2,\cdots, {1\over{2}}(N-\nu_{d}-\nu_{s}) $ with $ \nu_{s} = 0 $ or 1, where the normalization constant

The conventional phase factor $ (-1)^{\xi} $ shown in Eq. (5) for $ {{SU}}_{d} $(1,1) is adopted, which is consistent with the generalized pairing operator $ S^{\pm} = S^{+}_{s}- S^{\pm}_{d} $ used in Eq. (2). The matrix representations of $ {{SU}}_{d}(1,1)\otimes $ $ {{SU}}_{s} $(1,1) under the basis vectors (5) are given by

and

The eigenstate of Eq. (1) can be written as

where $ \alpha^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $ and $ \beta^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, in general, are complex numbers to be determined, ζ labels the ζ-th set of solution $ \{x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}; y^{(\zeta)}_{1},\cdots,y^{(\zeta)}_{k+1}\} $, $ \alpha^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, $ \beta^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, $ \vert \nu_{s};\nu_{d}\,\eta\, L\,M\rangle $ is the boson pairing vacuum state satisfying $ S^{-}_{\rho}\vert \nu_{s}; \nu_{d}\,\eta\, L\,M\rangle $ for $ \rho = s $ and d, in which $ \nu_{s}( = 0 $ or 1), and

which is equivalent to the form used in Ref. [5] with a linear transformation for x, where x is the spectral parameter to be determined. Using the commutation relations (4), one can directly verify that

There are also useful identities:

or

which can be proven using a mathematical induction on k, as $ S^{+}(x) $ and $ S^{+}(y) $ are binomials of $ S_{s}^{+} $ and $ S_{d}^{+} $. Using Eq. (17) with $ x_{k+1} = -1 $, we also have

Similar to the U(5)–O(6) case shown in Ref. [5], using the above commutation relations and identities, one can verify that

where the identities (18) and (12) for $ y_{j} $ within the summation over j are applied,

where the identity (19) is used.

Therefore, the eigen-equation

is fulfilled if and only if

and

Eqs. (27) and (28) are nothing but the eigen-equation

with

Thus, $ \alpha^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, $ \beta^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, and $ E^{(\zeta)}_{k,\nu_{s},\nu_{d},L} $ can be expressed in terms of the variables $ \{x^{(\zeta)}_{j}\} $ ($ j = 1,2,\cdots,k $) and $ \{y^{(\zeta)}_{j}\} $ ($ j = 1,2,\cdots,k+1 $). Once $ \alpha^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, $ \beta^{(\zeta)}_{\nu_{s},\nu_{d},\eta,L} $, and $ E^{(\zeta)}_{k,\nu_{s},\nu_{d},L} $ are expressed in terms of A, B, C, and D shown in Eq. (30), the variables $ \{x^{(\zeta)}_{j}\} $ ($ j = 1,2,\cdots,k $) and $ \{y^{(\zeta)}_{j}\} $ ($ j = 1,2,\cdots,k+1 $) are determined by Eqs. (25) and (26). It is evident that Eq. (25)–(28) become

or

when $ g_{s} = g_{d} = 0 $ without configuration mixing, which are the Bethe ansatz equations and the corresponding eigen-energy of the U(5)–O(6) transitional case for the N-boson normal states and the $ N+2 $-boson intruder states, respectively.

Similar to the results shown in Ref. [17], there are extended Heine-Stieltjes polynomials $ y^{(i)}(x) $ related to Eqs. (31) and (33) satisfying

for $ i = 1 $ or $ i = 2 $, where $ F(x) = x\,(1+x)^2 $,

and $ V^{(i)}(x) $ is a linear function of x determined by Eq. (35). The roots of Eqs. (31) or (33) are zeros of $ y^{(1)}(x) $ or $ y^{(2)}(x) $. Hence, the polynomial approach shown in Ref. [17] applies to this case as well, which can be used to obtain a solution of Eqs. (31) and (33) when there is no configuration mixing with $ g_{s} = g_{d} = 0 $.

It is evident that the pairing operators $ {\mathop \prod \nolimits_{\rho = 1}^k }S^{+}(x^{(\zeta)}_{\rho}) $ and $ {\mathop \prod \nolimits_{\rho = 1}^{k + 1} }S^{+}(y^{(\zeta)}_{\rho}) $ used in Eq. (9) are symmetric with respect to any permutation among $ \{x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}\} $ and $ \{y^{(\zeta)}_{1},\cdots,y^{(\zeta)}_{k+1}\} $. Therefore, there are $ k!(k+1)! $ identical roots of Eqs. (25) and (26), where only one is needed to construct the eigenstate (9). Once the ζ-th root $ \{x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}\} $ and $ \{y^{(\zeta)}_{1},\cdots,y^{(\zeta)}_{k+1}\} $ of Eqs. (25) and (26) are obtained, the pairing operators $ {\mathop \prod \nolimits_{\rho = 1}^k }S^{+}(x^{(\zeta)}_{\rho}) $ and $ {\mathop \prod \nolimits_{\rho = 1}^{k + 1} }S^{+}(y^{(\zeta)}_{\rho}) $ in Eq. (9) can be expressed as

where $ S^{(\mu)}(x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}) = \sum_{1\leq i_{1}<\cdots< i_{\mu}\leq k}\prod \nolimits_{q = 1}^{\mu}x^{(\zeta)}_{i_{q}} $ starting with $ S^{(0)}(x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}) = 1 $, and similarly for $ S^{(\mu)}(y^{(\zeta)}_{1},\cdots,y^{(\zeta)}_{k+1}) $, is the $ \mu $-th elementary symmetric polynomial of the k root-components $ \{x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}\} $ of Eq. (25), which is helpful for calculating matrix elements of physical quantities of the system.

The solution of Eq. (1) can be derived using the extended Heine-Stieltjes polynomial approach shown in Eq. (35). When there is no configuration mixing with $ g_{s} = g_{d} = 0 $, the roots $ \{x^{(\zeta)}_{01},\cdots,x^{(\zeta)}_{0k}\} $ and $ \{y^{(\zeta)}_{01},\cdots,y^{(\zeta)}_{0k+1}\} $ can be obtained, where the total number of roots is $ 2k+3 $, for which $ N = 2k+v_{d}+v_{s} $ for the allowed angular momentum quantum number L determined by the reduction rule $ (\nu_{d}\,0)\downarrow L $ of O(5)$\supset $O(3). Then, for small values of $ g_{s} $ and $ g_{d} $, the root $ \{x^{(\zeta)}_{1},\cdots,x^{(\zeta)}_{k}\} $ and $ \{y^{(\zeta)}_{1},\cdots,y^{(\zeta)}_{k+1}\} $ of Eqs. (25) and (26) for the given ζ is obtained using $ \{x^{(\zeta)}_{01},\cdots,x^{(\zeta)}_{0k}\} $ and $ \{y^{(\zeta)}_{01},\cdots,y^{(\zeta)}_{0k+1}\} $ as the initial root to determine the solution. Repeating this procedure, one can identify $ 2k+3 $ sets of the roots for any real values $ g_{s} $ and $ g_{d} $. Notably, all roots are real when $ g_{s} = g_{d} = 0 $, which is a common feature of the general SU(1,1) Gaudin models without previous study of the configuration mixing [5, 17, 18]. However, for $ g_{s}\neq 0 $ and $ g_{d}\neq0 $, complex roots occur in the middle part of the spectrum when the mixing of the N-boson and $ N+2 $-boson configurations is relatively strong, particularly when k is large, which is common when the configuration mixing strengths $ g_{s} $ and $ g_{d} $ become sufficiently large. Because -1 and 0 are singular points of Eqs. (25) and (26), the real part of the root components lies in the union $ (-\infty,-1)\bigcup(-1,0)\bigcup(0,\infty) $, where no pair of the root components is the same. In fact, similar to the solution of the pairing model [19], all root components are always symmetric with respect to the real axis on the complex plane; namely, if a root component is complex, the conjugate root component must be involved.

To demonstrate the feature of the solution, we consider an example with $ \alpha_{s}^{(1)} = 0 $, $ \alpha_{d}^{(1)} = 0.3 $ MeV, $ g^{(1)} = -0.5 $ MeV, $ \alpha_{s}^{(2)} = 1.5 $ MeV, $ \alpha_{d}^{(2)} = 1.8 $ MeV, $ g^{(2)} = -0.2 $ MeV, and $ g_{s} = g_{d} = 0.2 $ MeV. Only $ \nu_{s} = \nu_{d} = 0 $is exemplified in the following. As shown in Table 1, all roots in this case are real for $ N = 2 $. It is obvious that the first two roots mainly lie in the $ N = 2 $ configuration, as indicated by the small $ \beta^{(\zeta)}/\alpha^{(\zeta)} $ values, whereas the last three roots mainly lie in the $ N+2 = 4 $ configuration, as indicated by the relatively larger $ \beta^{(\zeta)}/\alpha^{(\zeta)} $ values. For $ N = 4 $, the pattern of the roots is similar. The first three roots mainly lie in the $ N = 4 $ configuration, whereas the last four roots mainly lie in the $ N+2 = 6 $ configuration. As shown in Table 2, the third root components with $ {y^{(3)}_{1} = y^{(3)*}_{2}} $ and the fifth root components with $ x^{(5)}_{1} = x^{(5)*}_{2} $ are complex in this case. In fact, with further increasing of the configuration mixing strengths $ g_{s} $ and $ g_{d} $, complex roots also occur for the $ N = 2 $ case. Because the two roots with relatively small $ \beta^{(\zeta)}/\alpha^{(\zeta)} $ values mainly lie in the $ N = 2 $ configuration, whereas the three roots with relatively larger $ \beta^{(\zeta)}/\alpha^{(\zeta)} $ values mainly lie in the $ N+2 = 4 $ configuration, the root of the ground state in this case becomes complex when $ g_{s} $ or $ g_{d} $ become sufficiently large. For example, if we keep other parameters the same as those shown in Table 1, the root of the ground state becomes complex when $ g_{s} = g_{d}\geqslant 6.06492 $ MeV with $ x_{1}^{(1)} = 0.56635 $, $y^{(1)}_{1} = 0.72673+ $$ 0.00027\,I $, $y^{(1)}_{2} = $$ 0.72673- 0.00027\,I $, $ \beta^{(1)}/\alpha^{(1)} = -0.3067 $, $ E^{(1)} = -14.8053 $ MeV for $ g_{s} = g_{d} = 6.06492 $ MeV. Therefore, the complex solution is likely to occur when the configuration mixing is sufficiently strong.

	$ x^{(\zeta)}_{1} $	$ y^{(\zeta)}_{1} $	$ y^{(\zeta)}_{2} $	$ \beta^{(\zeta)}/\alpha^{(\zeta)} $	$ E^{(\zeta)} $
$ \zeta=1 $	$ -1.21628 $	$ -1.23643 $	$ 0.92428 $	$ -0.03496 $	$ -0.92692 $
$ \zeta=2 $	$ 4.13257 $	$ 1.27972 $	$ 3.02724 $	$ -0.04105 $	$ 0.38300 $
$ \zeta=3 $	$ 0.38403 $	$ -2.88886 $	$ -1.23803 $	$ 4.99472 $	$ 4.45235 $
$ \zeta=4 $	$ -0.83687 $	$ -1.38586 $	$ 1.96679 $	$ 3.92177 $	$ 4.93409 $
$ \zeta=5 $	$ 3.65856 $	$ 0.66637 $	$ 9.05257 $	$ 2.33688 $	$ 5.88248 $

Table 1.  Roots of Eqs. (25) and (26), ratio $ \beta^{(\zeta)}/\alpha^{(\zeta)} $ used in eigen-state (9) and corresponding eigen-energy $ E^{(\zeta)} $ (in MeV) of the model for $ N = 2 $, where the parameters of Eq. (1) are considered as $ \alpha_{s}^{(1)} = 0 $, $ \alpha_{d}^{(1)} = 0.3 $ MeV, $ g^{(1)} = -0.5 $ MeV, $ \alpha_{s}^{(2)} = 1.5 $ MeV, $ \alpha_{d}^{(2)} = 1.8 $ MeV, $ g^{(2)} = -0.2 $ MeV, and $ g_{s} = g_{d} = 0.2 $ MeV.

	$ x^{(\zeta)}_{1} $	$ x^{(\zeta)}_{2} $	$ y^{(\zeta)}_{1} $	$ y^{(\zeta)}_{2} $	$ y^{(\zeta)}_{3} $	$ \beta^{(\zeta)}/\alpha^{(\zeta)} $	$ E^{(\zeta)} $
$ \zeta=1 $	$ -1.18956 $	$ -1.25809 $	$ -1.29710 $	$ -1.18038 $	$ 0.95864 $	$ -0.02389 $	$ -3.19779 $
$ \zeta=2 $	$ -1.12636 $	$ 3.76299 $	$ -1.14255 $	$ 1.17349 $	$ 3.10565 $	$ -0.02637 $	$ -1.94532 $
$ \zeta=3 $	$ 0.78722 $	$ \; 11.78400 $	$ 0.79059-0.5954{I} $	$ 0.79059+0.5954{I} $	$ \; 10.41310 $	$ -0.03314 $	$ 0.48860 $
$ \zeta=4 $	$ -2.51113 $	$ 1.06583 $	$ -3.9343 $	$ -1.55612 $	$ -1.05363 $	$ 3.72052 $	$ 4.72435 $
$ \zeta=5 $	$ -0.38146-1.1214I $	$ -0.38146+1.1214I $	$ -1.73495 $	$ -1.11698 $	$ 1.49812 $	$ 3.10755 $	$ 5.25399 $
$ \zeta=6 $	$ -1.13035 $	$ 3.19011 $	$ -1.24602 $	$ 0.613486 $	$ 8.04287 $	$ 2.37150 $	$ 6.15037 $
$ \zeta=7 $	$ 0.73826 $	$ \; 10.59300 $	$ 0.331345 $	$ 1.70441 $	$ 18.8986 $	$ 1.51102 $	$ 7.55080 $

Table 2.  Same as Table 1, but for $ N = 4 $, where $ I = \sqrt{-1}. $

To apply this theory, the Hamiltonian (1) is employed to describe the low-lying spectrum of $ ^{108} {\rm{Cd}}$ with $ N = 6 $ bosons, where the term $ \hat{H}_{L} = f L\cdot L $ is added to Eq. (1) to lift the degeneracy of the levels with the same seniority but different angular momentum quantum numbers. The E2 operator is chosen as

with which the B(E2) values are given by

where $ q_{2} $ and $ q^{\prime}_{2} $ are effective charge parameters of the normal and intruder configurations, respectively, and the reduced matrix element is defined in terms of the CG coefficient, such that $ \langle\zeta_{f},\nu^{\prime}_{d}\,\eta^{\prime}\,\eta\,L_{f}\vert\vert \hat{I}\vert\vert \zeta_{i},\nu_{d}\,L_{i}\rangle \!=\! \delta_{ \zeta_{f},\zeta_{i}}\delta_{ \nu^{\prime}_{d},\nu_{d}} \delta_{L_{f},L_{i}} $ with unit identity operator $ \hat{I} $.

Similar to that in Ref. [13], the level energies up to the three-phonon states in the normal bands and the intruder states $ 0^{+}_{1}({ i}) $, $ 2^{+}_{1}({ i}) $, and $ 4^{+}_{1}({ i}) $ of $ ^{108} {\rm{Cd}}$ deduced in Ref. [20] are considered. The model parameters are produced by a best global fit to the experimental level energies alone, where we obtain $ \alpha_{s}^{(1)} = 0 $, $ \alpha_{d}^{(1)} = 1.261 $ MeV, $ g^{(1)} = -1 $ keV, $ \alpha_{s}^{(2)} = 300 $ keV, $ \alpha_{d}^{(2)} = 1.416 $ MeV, $ g^{(2)} = -51 $ keV, $ g_{s} = 220 $ keV, $ g_{d} = 200 $ keV, $ f = 5 $ keV, and $ q'_{2}/q_{2} = -0.38 $. Then, the experimentally measured B(E2) ratios, $ R( L_{i}\rightarrow L_{f}) = B({E}2; L_{i}\rightarrow L_{f})/B({E}2; 2^{+}_{1}\rightarrow 0^{+}_{g} ) $, provided in Ref. [20] are fitted by only adjusting the ratio $ q_{2}^{\prime}/q_{2} $. The fitted low-lying level energies and B(E2) ratios are shown in Table 3, where the corresponding results of the $ 2n $-particle and $ 2n $-hole configuration mixing from $ n = 0 $ to $ n\rightarrow\infty $ in the U(5) limit of the IBM (CM5) [13] are likewise provided. The ratio $ q_{2}^{\prime}/q_{2} = 2.9 $ is mainly determined according to the lower limit of the experimental ratio $ R( 4^{+}_{1}({ i})\rightarrow 2^{+}_{1}({ i})) $. Regarding the level energies, the U(5)–O(6) transitional description is slightly better than the CM5, whereas the B(E2) ratios generated in the two models are quite the same for the transitions among normal states. However, although the E2 decays out of the intruder band are still weaker [20], the B(E2) ratio $ R( 2^{+}_{1}({ i})\rightarrow 0^{+}_{1}({ i})) $ and those for the transitions from the intruder states to the normal states predicted in this model are far larger than those of the CM5, as shown in III. Arguably, these values can be reduced when configuration mixing with $ N+2n $ bosons for $ n\geqslant 2 $ is considered. Because the E2 operator is simply selected as the generator of O(6), as shown in Eq. (38) and in the CM5 case [13], the E2 selection rules are similar to those given in the U(5) or the O(6) limit without configuration mixing, which are given by $ \Delta\nu_{d} = \pm 1 $. Therefore, the $ \Delta\nu_{d} = \pm 2 $ transitions, such as $ B({E}2; 2^{+}_{2}\rightarrow 0^{+}_{ g}) $, and the $ \Delta\nu_{d} = 0 $ transitions, such as $ B({E}2; 2^{+}_{1}({ i})\rightarrow 2^{+}_1) $ and $ B({E}2; 4^{+}_{1}({ i})\rightarrow 4^{+}_{1}) $, are consistently zero. To improve the theory, the O(6) symmetry breaking terms, such as the $ (d^{\dagger }\tilde{d})^{2}_{\mu} $ term, which allows $ \Delta\nu_{d} = 0 $ transitions, must be added in the E2 operator. Alternatively, high order interactions, such as those adopted in Ref. [14], may be considered.

Level energy/MeV	This study	Exp. [20, 21]	CM5 [13]	$ R( L_{i}\rightarrow L_{f}) $	This study	Exp. [20, 21]	CM5 [13]
$ E(2^{+}_{1}) $	0.789	$ 0.632 $	$ 0.718 $	$ R(4^{+}_{1}\rightarrow 2^{+}_1) $	$ 1.475 $	$ 1.5639 $	$ 1.6688 $
$ E(4^{+}_{1}) $	1.574	$ 1.508 $	$ 1.484 $	$ R(2^{+}_{2}\rightarrow 2^{+}_1) $	$ 1.475 $	$ 0.6579 $	$ 1.6688 $
$ E(2^{+}_{2}) $	1.504	$ 1.601 $	$ 1.470 $	$ R(2^{+}_{2}\rightarrow 0^{+}_{ g}) $	$ 0 $	$ 0.0676 $	$ 0 $
$ E(0^{+}_{2}) $	1.350	$ 1.913 $	$ 1.264 $	$ R(2^{+}_{1}({ i})\rightarrow 0^{+}_{1}({ i})) $	1.378	$ \geq0.338 $	$ 0.3428 $
$ E(3^{+}_{1}) $	2.234	$ 2.146 $	$ 2.268 $	$ R(4^{+}_{1}({ i})\rightarrow 2^{+}_{1}({ i})) $	0.213	$ \geq0.226 $	$ 0.5901 $
$ E(4^{+}_{2}) $	2.274	$ 2.239 $	$ 2.276 $	$ R(2^{+}_{1}({ i})\rightarrow 0^{+}_{ g}) $	0.338	$ \geq0.002 $	$ 0.0029 $
$ E(0^{+}_{3}) $	2.174	$ 2.375 $	$ 1.896 $	$ R(4^{+}_{1}({ i})\rightarrow 2^{+}_1) $	0.234	$ \geq0.005 $	$ 0.0043 $
$ E(2^{+}_{3}) $	2.059	$ 2.486 $	$ 1.982 $	$ R(2^{+}_{1}({ i})\rightarrow 2^{+}_1) $	$ 0 $	$ \geq0.015 $	$ 0 $
$ E(6^{+}_{1}) $	2.384	$ 2.541 $	$ 2.298 $	$ R(4^{+}_{1}({ i})\rightarrow 4^{+}_{1}) $	$ 0 $	$ \geq0.015 $	$ 0 $
$ E(0^{+}_{1}({ i})) $	1.755	$ 1.720 $	$ 1.720 $
$ E(2^{+}_{1}({ i})) $	2.367	$ 2.163 $	$ 2.438 $
$ E(4^{+}_{1}({ i})) $	2.797	$ 2.739 $	$ 3.084 $
$ E(2^{+}_{2}({ i})) $	2.727	$ 2.366 $	$ 3.070 $
$ E(0^{+}_{2}({ i})) $	2.639	$ 2.740{*} $	$ 2.984 $
$ E(0^{+}_{3}({ i})) $	2.856	$ 2.936{*} $	$ 2.984 $

Table 3.  Some low-lying level energies and B(E2) ratios $ R( L_{i}\rightarrow L_{f}) = B({E}2; L_{i}\rightarrow L_{f})/B({E}2; 2^{+}_{1}\rightarrow 0^{+}_{g} ) $ of $ ^{108} {\rm{Cd}}$, where $ {*} $ indicates that the corresponding spin assignment is not fully confirmed. The spin of both $ 0^{+}_{2}({ i}) $ and $ 0^{+}_{3}({ i}) $ states was assigned with $ (0^{+},1^{+}, 2^{+}) $, as shown in Ref. [21]. Model parameters are considered as $ \alpha_{s}^{(1)} = 0 $, $ \alpha_{d}^{(1)} = 1.261 $ MeV, $ g^{(1)} = -1 $ keV, $ \alpha_{s}^{(2)} = 300 $ keV, $ \alpha_{d}^{(2)} = 1.416 $ MeV, $ g^{(2)} = -51 $ keV, $ g_{s} = 220 $ keV, $ g_{d} = 200 $ keV, $ f = 5 $ keV, and $ q'_{2}/q_{2} = 2.9 $.

In this study, we demonstrate that the U(5)-O(6) transitional Hamiltonian of the interacting boson model with two-particle and two-hole configuration mixing is exactly solvable. An exact solution is derived based on the Bethe ansatz approach, where the Bethe ansatz equations are provided to determine the eigenstates and the corresponding eigen-energies . The solution features are numerically exemplified by the $ N = 2 $ and $ N = 4 $ cases. As an example of application, some low-lying level energies and B(E2) ratios of $ ^{108} {\rm{Cd}}$ are fitted and compared with the corresponding experimental data.

Because the solution of the Hamiltonian without configuration mixing can be easily derived using the extended Heine-Stieltjes polynomials, the roots of the Bethe ansatz equations for cases with small mixing parameters are approximately found using the roots of the equations without configuration mixing as the initial values. Therefore, a progressive approach can be established to obtain the solution of the model with arbitrary mixing parameters. Although the solution is only demonstrated for the $ N\oplus (N+2) $ configuration mixing, it is expected that the model with $ 2n $-particle and $ 2n $-hole configuration mixing for $ n = 0 $ up to a finite n is also exactly solvable by using the identities and procedures shown in Sec. 2. This is because the eigenstates of the model can always be expressed in terms of binomials of s- and d-boson pair operators, although the equations involved become more complicated. A similar extension to the IBM-II case is likewise straightforward. Moreover, a chain of isotopes or isotones in the vibrational to γ-soft transitional region may be analyzed using the model to reveal their shape phase coexistence and evolution, for example, as the analysis for Cd isotopes shown in [22–24], which will be considered in our future study.

One of the authors (F. Pan) is grateful to Professor P. Van Isacker for stimulating discussions on the subject and his suggestion on conducting this work.