1. Introduction
With the development of radioactive ion beam facilities and detection techniques around the world, the study of nuclei far from the β stability line has attracted wide attention [1–5]. New and exotic phenomena have been observed in nuclei near the drip lines, such as halo nuclei [6], changes in nuclear magic numbers [7], the island of inversion [8], pygmy resonances [9], etc. Several exotic nuclear phenomena have also been predicted, for example, giant halos [10, 11], shape decoupling between the core and halo [12, 13], etc. In exotic nuclei, the extremely weak binding of valence nucleons leads to many new features, such as coupling between bound states and continuum due to pairing correlations and the very extended spatial density distributions. The description of exotic nuclei which are weakly bound and have very extended spatial density distributions has become a very challenging topic in nuclear physics.
The density functional theory (DFT) with a small number of parameters is widely used to study the ground- and excited-state properties of the nuclei all over the nuclear chart [14]. In particular, the covariant density functional theory (CDFT), incorporating the Lorentz symmetry in a self-consistent way, has received wide attention for its successful description of many nuclear phenomena [15–19]. It includes naturally the nucleonic spin degree of freedom and automatically results in the nuclear spin–orbit potential with empirical strength in a covariant way. It can reproduce well the isotopic shifts in the Pb region [20], and give naturally the origin of the pseudospin symmetry [21, 22] and the spin symmetry in the anti-nucleon spectrum [23, 24]. Furthermore, it can include nuclear magnetism [25], that is, a consistent description of currents and time-odd fields, which plays an important role in nuclear magnetic moments [26–28] and nuclear rotations [29–32]. Besides, in the framework of CDFT, a self-consistent quasiparticle random phase approximation approach is developed and applied to investigations on β decay successfully [33, 34].
In the framework of CDFT, by extending the relativistic mean field theory with the Bogoliubov transformation in the coordinate representation, the relativistic continuum Hartree–Bogoliubov (RCHB) theory was developed [35, 36] and it provides a proper treatment of pairing correlations and mean-field potentials in the presence of the continuum. With the RCHB theory, the first microscopic self-consistent description of halo in 11Li has been provided [35] and the giant halos in light and medium-heavy nuclei have been predicted [10, 11]. The RCHB theory has been generalized to treat the odd nucleon system [37], and together with the Glauber model, the charge-changing cross sections from carbon (C) to fluorine (F) isotopes on a carbon target have been well reproduced [38]. Focusing on the continuum effects on the nuclear landscape, a systematic calculation for the nuclear landscape in the framework of RCHB theory with PC-PK1 is performed [39, 40]. The RCHB mass table has been applied to investigate α-decay energy [41], shell evolution [42] and proton radioactivity [43].
For deformed nuclei, a deformed relativistic Hartree–Bogoliubov theory in continuum (DRHBc) theory has been developed [12, 13]. By solving the relativistic Hartree–Bogoliubov (RHB) equation with the Dirac Woods–Saxon basis [44], the DRHBc theory could treat the axial deformation, continuum and pairing correlations self-consistently. The DRHBc theory is useful to answer questions concerning exotic nuclear phenomena in nuclei near the drip line, such as whether there are deformed halos or not and what new features can be expected in deformed exotic nuclei. As a first application, halo phenomena in deformed nuclei were investigated within the DRHBc theory, and an interesting shape decoupling between the core and the halo was predicted[12, 13]. Later, the DRHBc theory has been extended to incorporate the density-dependent meson-nucleon couplings [45]. In [46], the DRHBc theory is extended to incorporate the blocking effect due to an odd nucleon. The DRHBc theory has achieved success in resolving the puzzles concerning the radius and configuration of valence neutrons in 22C [47] and studying particles in the classically forbidden regions for magnesium isotopes [48]. Recently, the mass table based on the DRHBc theory is performing [49]
The triaxial deformation of the nucleus has been a subject of much interest in the study of nuclear structure for a long time. It plays important roles in investigations on γ-band [50], chiral doublet bands [51–53], wobbling motion [50, 54, 55], signature inversion [56], anomalous signature splitting [57] and superdeformed triaxial bands [58]. Especially, the discovery of chiral doublet band and wobbling motion, which are evidence of triaxial deformation, has stimulated the study of triaxial nuclei [59, 60]. The inclusion of triaxiality can dramatically reduce the barrier separating prolate and oblate minima, leading to structures that are soft or unstable for triaxial distortions [61]. Triaxial deformation is essential for the investigations on the fission barrier [62–64]. Furthermore, the role of triaxiality in deformed halo nuclei is discussed [65].
Theoretically, the systematical calculation based on the macroscopic-microscopic finite-range droplet model has illustrated that specific combinations of single-particle orbitals near the Fermi surface and the additional binding energy from non-axial degrees of freedom can enhance the tendency to form nuclei with triaxial shapes, and some islands of triaxiality have been revealed throughout the nuclear landscape [66]. Using non-relativistic DFT, [67–71] study the triaxial deformation in the nuclei with mass numbers A ∼ 100, 130 and 190. In the framework of CDFT, based on the relativistic Hartree–Bogoliubov model and relativistic mean-field model, the ground-state and low-lying excited state properties of nuclei with A ∼ 80, 100 and 130 are investigated, and the roles of triaxial deformation are discussed [72–80]. However, in these calculations the harmonic-oscillator (HO) basis is used in solving the RHB equation or Dirac equation, due to the incorrect asymptotic behavior of the HO wave functions, the expansion in a localized HO basis is not appropriate for the description of exotic nuclei [35, 36, 81, 82].
Considering the success of RCHB theory and DRHBc theory in the description of exotic nuclei, the present work developed triaxially deformed relativistic Hartree–Bogoliubov theory in Woods–Saxon basis (TDRHBWS) which could treat the triaxial deformation, pairing correlations and continuum.
The paper is organized as follows. In section 2 , the theoretical framework for the triaxially deformed relativistic Hartree–Bogoliubov theory in Woods–Saxon basis is briefly presented. The numerical check and details are given in section 3 . In section 4 , the separation energy, deformation, radii, single-particle levels and potential energy surface calculated by TDRHBWS are presented. A short summary is given in section 5 .
2. Theoretical framework
The starting point of covariant DFT is a Lagrangian density of form
$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L } & = & \bar{\psi }({\rm{i}}{\gamma }_{\mu }{\partial }^{\mu }-M)\psi -\displaystyle \frac{1}{2}{\alpha }_{S}(\bar{\psi }\psi )(\bar{\psi }\psi )\\ & & -\displaystyle \frac{1}{2}{\alpha }_{V}(\bar{\psi }{\gamma }_{\mu }\psi )(\bar{\psi }{\gamma }^{\mu }\psi )-\displaystyle \frac{1}{2}{\alpha }_{{TV}}(\bar{\psi }\vec{\tau }{\gamma }_{\mu }\psi )(\bar{\psi }\vec{\tau }{\gamma }^{\mu }\psi )\\ & & -\displaystyle \frac{1}{2}{\alpha }_{{TS}}(\bar{\psi }\vec{\tau }\psi )(\bar{\psi }\vec{\tau }\psi )-\displaystyle \frac{1}{3}{\beta }_{S}{\left(\bar{\psi }\psi \right)}^{3}\\ & & -\displaystyle \frac{1}{4}{\gamma }_{S}{\left(\bar{\psi }\psi \right)}^{4}-\displaystyle \frac{1}{4}{\gamma }_{V}{[(\bar{\psi }{\gamma }_{\mu }\psi )(\bar{\psi }{\gamma }^{\mu }\psi )]}^{2}\\ & & -\displaystyle \frac{1}{2}{\delta }_{S}{\partial }_{\nu }(\bar{\psi }\psi ){\partial }^{\nu }(\bar{\psi }\psi )-\displaystyle \frac{1}{2}{\delta }_{V}{\partial }_{\nu }(\bar{\psi }{\gamma }_{\mu }\psi ){\partial }^{\nu }(\bar{\psi }{\gamma }^{\mu }\psi )\\ & & -\displaystyle \frac{1}{2}{\delta }_{{TV}}{\partial }_{\nu }(\bar{\psi }\vec{\tau }{\gamma }_{\mu }\psi ){\partial }^{\nu }(\bar{\psi }\vec{\tau }{\gamma }_{\mu }\psi )\\ & & -\displaystyle \frac{1}{2}{\delta }_{{TS}}{\partial }_{\nu }(\bar{\psi }\vec{\tau }\psi ){\partial }^{\nu }(\bar{\psi }\vec{\tau }\psi )\\ & & -\displaystyle \frac{1}{4}{F}^{\mu \nu }{F}_{\mu \nu }-e\displaystyle \frac{1-{\tau }_{3}}{2}\bar{\psi }{\gamma }^{\mu }\psi {A}_{\mu },\end{array}\end{eqnarray}$
where M is the nucleon mass; Aμ and Fμν are respectively the four-vector potential and field strength tensor of the electromagnetic field. Here αS, αV, αTS and αTV represent the coulpling constants for four-fermion point-coulping terms, βS, γS and γV are those for the higher-order terms which are responsible for the effects of medium dependence, δS, δV, δTS and δTV refer to those for the gradient terms which are included to simulate the finite range effects. The subscripts S, V and T stand for scalar, vector and isovector, respectively.Pairing correlations are crucial in the description of open-shell nuclei. For exotic nuclei, the conventional BCS approach turns out to be only a poor approximation [81, 83]. Starting from the Lagrangian density (1), a relativistic theory of pairing correlations in nuclei has been developed by Kucharek and Ring [84]. If we neglect the Fock terms as it is usually done in the covariant DFT, the RHB equation for the nucleons reads
$\begin{eqnarray}\left(\begin{array}{cc}{\hat{h}}_{D}-\lambda & \hat{{\rm{\Delta }}}\\ -{\hat{{\rm{\Delta }}}}^{* } & -{\hat{h}}_{D}+\lambda \end{array}\right)\left(\displaystyle \genfrac{}{}{0em}{}{{U}_{k}}{{V}_{k}}\right)={E}_{k}\left(\displaystyle \genfrac{}{}{0em}{}{{U}_{k}}{{V}_{k}}\right),\end{eqnarray}$
where Ek is the quasiparticle energy, λ the chemical potential, and hD is the Dirac Hamiltonian $\begin{eqnarray}{h}_{D}=[{\boldsymbol{\alpha }}\cdot {\boldsymbol{p}}+V({\boldsymbol{r}})+\beta (M+S({\boldsymbol{r}}))].\end{eqnarray}$
The scalar potential S(r) and vector potential V(r) $\begin{eqnarray}\begin{array}{rcl}S({\boldsymbol{r}}) & = & {\alpha }_{S}{\rho }_{S}+{\beta }_{S}{\rho }_{S}^{2}+{\gamma }_{S}{\rho }_{S}^{3}+{\delta }_{S}\bigtriangleup {\rho }_{S},\\ {V}^{\mu }({\boldsymbol{r}}) & = & {\alpha }_{V}{j}^{\mu }+{\gamma }_{V}{\left({j}^{\mu }\right)}^{3}+{\delta }_{V}\bigtriangleup {j}^{\mu }\\ & & +{{eA}}^{\mu }+{\alpha }_{{TV}}{\tau }_{3}{\vec{j}}_{{TV}}^{\mu }+{\delta }_{{TV}}{\tau }_{3}\bigtriangleup {\vec{j}}_{{TV}}^{\mu },\end{array}\end{eqnarray}$
depend on the scalar density ρS, vector density ρV and iso-vector density ρTV, and the densities can be constructed by quasiparticle wave functions, $\begin{eqnarray}\begin{array}{rcl}{\rho }_{S}({\boldsymbol{r}}) & = & \displaystyle \sum _{k\gt 0}{\bar{V}}_{k}({\boldsymbol{r}}){V}_{k}({\boldsymbol{r}}),\\ {\rho }_{V}({\boldsymbol{r}}) & = & \displaystyle \sum _{k\gt 0}{V}_{k}^{\dagger }({\boldsymbol{r}}){V}_{k}({\boldsymbol{r}}),\\ {\rho }_{{TV}}({\boldsymbol{r}}) & = & \displaystyle \sum _{k\gt 0}{V}_{k}^{\dagger }({\boldsymbol{r}}){\tau }_{3}{V}_{k}({\boldsymbol{r}}).\end{array}\end{eqnarray}$
The pairing potential is determined by
$\begin{eqnarray}\begin{array}{l}{\rm{\Delta }}({{\boldsymbol{r}}}_{1}{s}_{1}{p}_{1},{{\boldsymbol{r}}}_{2}{s}_{2}{p}_{2})\\ \quad =\displaystyle \sum _{s{{\prime} }_{1}p{{\prime} }_{1}s{{\prime} }_{2}p{{\prime} }_{2}}{V}^{\mathrm{pp}}({{\boldsymbol{r}}}_{1},{{\boldsymbol{r}}}_{2};{s}_{1}{p}_{1},{s}_{2}{p}_{2},s{{\prime} }_{1}p{{\prime} }_{1},s{{\prime} }_{2}p{{\prime} }_{2})\\ \quad \times \kappa ({{\boldsymbol{r}}}_{1}s{{\prime} }_{1}p{{\prime} }_{1},{{\boldsymbol{r}}}_{2}s{{\prime} }_{2}p{{\prime} }_{2}),\end{array}\end{eqnarray}$
where p = 1, 2 is used to represent the large and small components of the Dirac spinors, s is used to represent the spin. Vpp is the effective pairing interaction and $\kappa ({{\boldsymbol{r}}}_{1}s{{\prime} }_{1}p{{\prime} }_{1},{{\boldsymbol{r}}}_{2}s{{\prime} }_{2}p{{\prime} }_{2})$ is the pairing tensor.For triaxially deformed nuclei, we expand the potentials in equation (4 ) and densities in equation (5 ) in terms of spherical harmonic functions,
$\begin{eqnarray}\begin{array}{rcl}f(r,\theta ,\phi ) & = & \displaystyle \sum _{\lambda \mu }{f}_{\lambda \mu }(r){Y}_{\lambda \mu }(\theta ,\phi ),\\ \lambda & = & 0,2,4...,\,\mu =-\lambda \,,-\lambda +2...\,\lambda ,\end{array}\end{eqnarray}$
with $\begin{eqnarray}{f}_{\lambda \mu }(r)=\int {\rm{d}}{\rm{\Omega }}{Y}_{\lambda \mu }^{* }(\theta ,\phi )f({\boldsymbol{r}}).\end{eqnarray}$
For triaxially deformed nuclei, the densities and potentials remain the same under spatial reflection transformation, (xy), (yz) and (zx) specular reflection transformations, i.e. V(r, θ, φ) = V(r, π − θ, π + φ), V(x, y, z) = V(x, y, −z), V(x, y, z) = V( − x, y, z) and V(x, y, z) = V(x, −y, z). For spatial reflection transformation7 ). For (xy) specular reflection transformations7 ). Similarly, (yz) and (zx) specular reflection symmetries require fλμ = fλ−μ in equation (7 ).
$\begin{eqnarray}\begin{array}{rcl}\hat{P}V(r,\theta ,\phi ) & = & \hat{P}\displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){Y}_{\lambda \mu }(\theta ,\phi )\\ & = & \displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){Y}_{\lambda \mu }(\pi -\theta ,\pi +\phi )\\ & = & \displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){\left(-1\right)}^{\lambda }{Y}_{\lambda \mu }(\theta ,\phi ),\end{array}\end{eqnarray}$
therefore, only even λ is taken in equation ( $\begin{eqnarray}\begin{array}{rcl}{\hat{P}}_{z}V(r,\theta ,\phi ) & = & {\hat{P}}_{z}\displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){Y}_{\lambda \mu }(\theta ,\phi )\\ & = & \displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){Y}_{\lambda \mu }(\pi -\theta ,\phi )\\ & = & \displaystyle \sum _{\lambda \mu }{V}_{\lambda \mu }(r){\left(-1\right)}^{\lambda +\mu }{Y}_{\lambda \mu }(\theta ,\phi ),\end{array}\end{eqnarray}$
the symmetry require λ + μ should be even, thus only even μ is taken in equation (The quasiparticle wave functions Uk and Vk are expanded in terms of Dirac Woods–Saxon (DWS) basis φ(rσp) obtained from the solution of Dirac equation hD containing spherical potential S0(r) and V0(r) of Woods–Saxon shape [44]
$\begin{eqnarray}{U}_{k}({\boldsymbol{r}}\sigma p)=\displaystyle \sum _{n\kappa m}{u}_{n\kappa m}^{k}{\varphi }_{n\kappa m}({\boldsymbol{r}}\sigma p),\end{eqnarray}$
$\begin{eqnarray}{V}_{k}({\boldsymbol{r}}\sigma p)=\displaystyle \sum _{n\kappa m}{v}_{n\kappa m}^{k}{\tilde{\varphi }}_{n\kappa m}({\boldsymbol{r}}\sigma p),\end{eqnarray}$
where ${\tilde{\varphi }}_{n\kappa m}$ is the time-reversal state of φnκm. These states form a complete spherical and discrete basis in Dirac space. Then we obtain the RHB matrix in DWS basis $\begin{eqnarray}\begin{array}{l}\left(\begin{array}{cc}{H}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m} & {{\rm{\Delta }}}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}\\ -{{\rm{\Delta }}}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}^{* } & -{H}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}\end{array}\right)\left(\begin{array}{c}{u}_{n\kappa m}^{k}\\ {v}_{n\kappa m}^{k}\end{array}\right)\\ \quad ={E}_{k}\left(\begin{array}{c}{u}_{n^{\prime} \kappa ^{\prime} m^{\prime} }^{k}\\ {v}_{n^{\prime} \kappa ^{\prime} m^{\prime} }^{k}\end{array}\right),\end{array}\end{eqnarray}$
where ${H}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}$ is the matrix element of Dirac Hamiltonian hD. $\begin{eqnarray}\begin{array}{l}{H}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}=\int {\rm{d}}{\boldsymbol{r}}{\varphi }_{n\mbox{'}\kappa \mbox{'}m\mbox{'}}^{\dagger }{h}_{D}{\varphi }_{n\kappa m}\\ \quad =\,{\epsilon }_{n^{\prime} \kappa ^{\prime} m^{\prime} }{\delta }_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}\\ \quad +\int {\rm{d}}r\displaystyle \sum _{\lambda \mu }\{{G}_{n^{\prime} \kappa ^{\prime} }{G}_{n\kappa }[{V}_{\lambda \mu }(r)\\ \quad +{S}_{\lambda \mu }(r)]A(\kappa ^{\prime} m^{\prime} ,\lambda \mu ,\kappa m)\\ \quad +{F}_{n^{\prime} \kappa ^{\prime} }{F}_{n\kappa }[{V}_{\lambda \mu }(r)-{S}_{\lambda \mu }(r)]\\ \quad \times \,A(-\kappa ^{\prime} m^{\prime} ,\lambda \mu ,-\kappa m)\}\\ \quad -\int {\rm{d}}r\{{G}_{n^{\prime} \kappa ^{\prime} }{G}_{n\kappa }[{V}^{0}(r)+{S}^{0}(r)]\\ \quad +{F}_{n^{\prime} \kappa ^{\prime} }{F}_{n\kappa }[{V}^{0}(r)-{S}^{0}(r)]\}{\delta }_{\kappa ^{\prime} \kappa }{\delta }_{m^{\prime} m},\,\,\end{array}\end{eqnarray}$
with $\begin{eqnarray}\begin{array}{c}A({\kappa }^{{\rm{{\prime} }}}{m}^{{\rm{{\prime} }}},\lambda \mu ,\kappa m)=\displaystyle \sum _{\sigma }\int {\rm{d}}{\rm{\Omega }}{{ \mathcal Y }}_{\kappa \text{'}m\text{'}}^{l\text{'}\ast }(\theta ,\phi ,\sigma )\\ \,\times {Y}_{\lambda \mu }(\theta ,\phi ){{ \mathcal Y }}_{\kappa m}^{l}(\theta ,\phi ,\sigma )\end{array}\end{eqnarray}$
represents angular part of the matrix element, and ${{ \mathcal Y }}_{{jm}}^{l}(\theta ,\phi ,\sigma )$ are the spinor spherical harmonics. Here Gnκ and Fnκ represent radial wave functions of DWS basis, V0(r) and S0(r) represent the spherical potential of Woods–Saxon potential.${{\rm{\Delta }}}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}$ is the pairing matrix element. For a fixed-gap pairing force, ${{\rm{\Delta }}}_{n^{\prime} \kappa ^{\prime} m^{\prime} ,n\kappa m}$ is a diagonal matrix.
Through diagonalizing the RHB matrix, quasiparticle energy and wave functions are obtained. Then new densities and fields are calculated, which are iterated in the matrix equations until convergence is achieved. Finally, one can calculate the total energy of a nucleus by
$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{RHB}} & = & \displaystyle \sum _{k}(\lambda -{E}_{k}){V}_{k}^{2}-{E}_{\mathrm{pair}}\\ & & -\int {{\rm{d}}}^{3}r\{\frac{1}{2}{\alpha }_{S}{\rho }_{S}^{2}+\frac{1}{2}{\alpha }_{V}{\rho }_{V}^{2}+\frac{1}{2}{\alpha }_{{TV}}{\rho }_{{TV}}^{2}\\ & & +\frac{2}{3}{\beta }_{S}{\rho }_{S}^{3}+\frac{3}{4}{\gamma }_{S}{\rho }_{S}^{4}+\frac{3}{4}{\gamma }_{V}{\rho }_{V}^{4}+\frac{1}{2}{\delta }_{S}{\rho }_{S}{\rm{\Delta }}{\rho }_{S}\\ & & +\frac{1}{2}{\delta }_{V}{\rho }_{V}{\rm{\Delta }}{\rho }_{V}+\frac{1}{2}{\delta }_{{TV}}{\rho }_{{TV}}{\rm{\Delta }}{\rho }_{{TV}}+\frac{1}{2}{{eA}}_{0}{\rho }_{p}\}\\ & & +{E}_{{\rm{c}}.{\rm{m}}.},\end{array}\end{eqnarray}$
where Epair is the pairing energy and Ec.m. is the center-of-mass correction energy.Usually proper treatment of center-of-mass (c.m.) correction energy is the microscopic approach17 ) is time-consuming, the phenomenological formulas
$\begin{eqnarray}{E}_{{\rm{c}}.{\rm{m}}.}=-\displaystyle \frac{1}{2{mA}}\langle {\hat{{\boldsymbol{P}}}}_{{\rm{c}}.{\rm{m}}.}^{2}\rangle \end{eqnarray}$
with A the mass number and ${\hat{{\boldsymbol{P}}}}_{{\rm{c}}.{\rm{m}}.}$ the total momentum in the c.m. frame. As the microscopic calculation of Ec.m. in equation ( $\begin{eqnarray}{E}_{{\rm{c}}.{\rm{m}}.}=-\displaystyle \frac{3}{4}41{A}^{-1/3}\,\mathrm{MeV}\end{eqnarray}$
is often used. In this work, the phenomenological formula is adopted.The root-mean-square (rms) radius is calculated as
$\begin{eqnarray}\begin{array}{rcl}{R}_{\tau } & = & \langle {r}^{2}{\rangle }^{1/2}={\left(\displaystyle \frac{1}{{N}_{\tau }}\int {{\rm{d}}}^{3}{\boldsymbol{r}}[{r}^{2}{\rho }_{V}^{\tau }({\boldsymbol{r}})]\right)}^{1/2}\\ & = & {\left(\displaystyle \frac{\sqrt{4\pi }}{{N}_{\tau }}\int {\rm{d}}{r}[{r}^{4}{\rho }_{V}^{\tau ,\lambda =0,\mu =0}(r)]\right)}^{1/2},\end{array}\end{eqnarray}$
where τ represents the proton or neutron and Nτ represents the number of proton or neutron. The rms charge radius is calculated as $\begin{eqnarray}{R}_{{\rm{c}}}=\sqrt{{R}_{p}^{2}+0.64\,{\mathrm{fm}}^{2}}.\end{eqnarray}$
The intrinsic multipole moment Q20 and Q22 are calculated by,
$\begin{eqnarray}{Q}_{20}=\sqrt{\displaystyle \frac{5}{16\pi }}\langle 2{z}^{2}-{x}^{2}-{y}^{2}\rangle =\langle {r}^{2}{Y}_{20}\rangle ,\end{eqnarray}$
$\begin{eqnarray}{Q}_{22}=\sqrt{\displaystyle \frac{15}{32\pi }}\langle {x}^{2}-{y}^{2}\rangle =\displaystyle \frac{1}{2}\langle {r}^{2}{Y}_{22}+{r}^{2}{Y}_{2-2}\rangle .\end{eqnarray}$
Using
$\begin{eqnarray}{Q}_{20}=\displaystyle \frac{3A}{4\pi }{R}_{0}^{2}{\beta }_{20},\end{eqnarray}$
$\begin{eqnarray}{Q}_{22}=\displaystyle \frac{3A}{4\pi }{R}_{0}^{2}{\beta }_{22},\end{eqnarray}$
the quadrupole deformation parameters β and γ are obtained $\begin{eqnarray}\beta =\sqrt{{\beta }_{20}^{2}+2{\beta }_{22}^{2}},\end{eqnarray}$
$\begin{eqnarray}\gamma =\arctan [\sqrt{2}\displaystyle \frac{{\beta }_{22}}{{\beta }_{20}}],\end{eqnarray}$
with R0 = 1.2A1/3 fm.3. Numerical details
3.1. Self-consistency numerical check
In this part, we check the self-consistency of the numerical code for TDRHBWS. Quadrupole deformed nucleus has D2 symmetry, when deformation parameter γ changes 120 degrees, the ground-state properties remain the same.
In table 1, the bulk properties of 76Ge with deformation being constrained to (β = 0.3, γ = 30°) and (β = 0.3, γ = 150°) are compared, including total energy, neutron pairing energy, proton pairing energy, neutron Fermi surface, proton Fermi surface, neutron rms radius, proton rms radius, matter rms radius and charge radius. In the calculations, the density functional PC-PK1 is used. In the pp channel, the fixed gap pairing force with Δ = 12A−1/2 MeV [85] is used. It shows that the bulk properties with γ being constrained to 30° and 150° are in good agreement. The difference of total energy between these two calculations is 0.0002 MeV, and the differences of rms radii are less than 10−4 fm.
Table 1. The bulk properties of 76Ge in TDRHBWS theory with deformations being constrained to (β = 0.3, γ = 30°) and (β = 0.3, γ = 150°), respectively. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
| 76Ge | γ = 30° | γ = 150° | δ |
|---|---|---|---|
| β | 0.3000 | 0.3000 | |
| γ (deg) | 29.9997 | 149.9997 | |
| Rn (fm) | 4.2571 | 4.2571 | 0.0000 |
| Rp (fm) | 4.0488 | 4.0488 | 0.0000 |
| Rm (fm) | 4.1706 | 4.1706 | 0.0000 |
| Rc(fm) | 4.1271 | 4.1270 | 0.0001 |
| λn(MeV) | −7.8045 | −7.8045 | 0.0007 |
| λp (MeV) | −10.3645 | −10.3646 | 0.0001 |
| ${E}_{\mathrm{pair}}^{n}$ (MeV) | −5.8226 | −5.8225 | 0.0001 |
| ${E}_{\mathrm{pair}}^{p}$ (MeV) | −4.2730 | −4.2730 | 0.0000 |
| Etot (MeV) | −656.3645 | −656.3647 | 0.0002 |
3.2. Comparison with DRHBc
Taking 24Ne as an example, we compare the results of the TDRHBWS calculation with those from the DRHBc calculation. In these two calculations, the density functional PC-PK1 and the fixed gap pairing force with Δ = 12A−1/2 MeV are used, box size Rbox = 15 fm, energy cutoff for DWS basis Ecut = 120 MeV, angular momentum cutoff ${J}_{\max }=$ 19/2ℏ. The same initial deformation β0 = 0.3 is used in both calculations. In the DRHBc calculation, the truncation of Legrende expansion ${\lambda }_{\max }=$ 4. In the TDRHBWS calculation, the truncation of spherical harmonics expansion ${\lambda }_{\max }=$ 4. In table 2, the bulk properties from these calculations, including total energy, neutron pairing energy, proton pairing energy, neutron Fermi surface, proton Fermi surface, neutron rms radius, proton rms radius, matter rms radius and charge radius, are listed. It can be seen that TDRHBWS calculations and DRHBc calculations converge to the same deformation, and the bulk properties agree with each other well. The difference in total energy between these two calculations is 0.0042 MeV, and the differences in rms radii are less than 10−4 fm.
Table 2. The ground-state properties of 24Ne from TDRHBWS and DRHBc calculations, respectively. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
| 24Ne | DRHBc | TDRHBWS | δ |
|---|---|---|---|
| β | 0.2494 | 0.2494 | 0.0000 |
| γ (deg) | 0.0000 | ||
| Rn (fm) | 3.1132 | 3.1132 | −0.0000 |
| Rp (fm) | 2.8684 | 2.8684 | −0.0000 |
| Rm (fm) | 3.0136 | 3.0136 | −0.0000 |
| Rc(fm) | 2.9779 | 2.9779 | −0.0000 |
| λn (MeV) | −6.9175 | −6.9175 | 0.0000 |
| λp (MeV) | −12.6198 | −12.6199 | 0.0001 |
| ${E}_{\mathrm{pair}}^{n}$ (MeV) | −5.4425 | −5.4425 | −0.0000 |
| ${E}_{\mathrm{pair}}^{p}$ (MeV) | −4.2760 | −4.2761 | 0.0001 |
| Etot (MeV) | −188.2991 | −188.3033 | 0.0042 |
Furthermore, we compare the single-particle energy of 24Ne from TDRHBWS and DRHBc calculations in figure 1. From the figure, it is found that both neutron and proton single-particle levels from the TDRHBWS calculation agree with the results from the DRHBc calculation.
Figure 1. The single-particle levels for 24Ne in canonical basis from the TDRHBWS calculation, in comparison with results from the DRHBc calculation. |
3.3. Comparison with harmonic-oscillator basis
In this section, taking triaxially deformed nucleus 76Ge as an example, we compare the results of the TDRHBWS calculation with the ones of the relativistic Hartree–Bogoliubov theory based on the harmonic oscillator basis (RHB-HO) to check the validity of the code. The density functional PC-PK1 and the fixed gap pairing force with Δ = 12A−1/2 MeV are used in these two calculations. In the TDRHBWS calculation, the box size Rbox = 15 fm, cutoff for DWS basis Ecut = 120 MeV, angular momentum cutoff ${J}_{\max }=$ 19/2ℏ and the truncation of spherical harmonics expansion ${\lambda }_{\max }=$ 6. In the RHB-HO calculation, the major shell number Nshell = 10. The same initial deformation β0 = 0.2, γ0 = 30° is used in both calculations. The corresponding results are listed in table 3.
Table 3. The ground-state properties of 76Ge in TDRHBWS theory and RHB-HO theory. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
| 76Ge | RHB-HO | TDRHBWS | δ |
|---|---|---|---|
| β | 0.1820 | 0.1825 | 0.0005 |
| γ (deg) | 6.8800 | 5.9870 | 0.8930 |
| Rn (fm) | 4.1720 | 4.2163 | −0.0443 |
| Rp (fm) | 3.9803 | 4.0182 | −0.0380 |
| Rm (fm) | 4.0923 | 4.0971 | −0.0048 |
| Rc(fm) | 4.0599 | 4.0846 | −0.0247 |
| λn (MeV) | −8.4979 | −8.4202 | −0.0772 |
| λp (MeV) | −9.3639 | −9.3513 | 0.0126 |
| ${E}_{\mathrm{pair}}^{n}$ (MeV) | −3.5203 | −5.8356 | 2.3153 |
| ${E}_{\mathrm{pair}}^{p}$ (MeV) | −3.0485 | −4.7992 | 1.7507 |
| Etot (MeV) | −658.4101 | −657.2466 | −1.1635 |
It shows that TDRHBWS calculations and RHB-HO calculations almost converge to the same deformation. The rms radii agree with each other. The differences in rms radii are about 10−2 fm. The difference in total energy between these two calculations is 1.1635 MeV, and the relative difference is about 0.1%.
In figure 2, the single-particle levels from TDRHBWS and RHB-HO calculations are shown. It can be seen that single-particle levels from the TDRHBWS calculation are in good agreement with the results from the RHB-HO calculation.
Figure 2. The single-particle levels for 76Ge in canonical basis from the TDHBWS calculation, in comparison with results from the RHB-HO calculation. |
3.4. Convergence check
For numerical reasons, several parameters have to be introduced in the calculations, including box size Rbox for the determination of the basis wave functions by solving the spherical Dirac equation, the energy cutoff Ecut for determination of the dimension of DWS basis, the maximal λ value ${\lambda }_{\max }$ in the expansion of the deformed fields and densities, and angular momentum cutoff ${J}_{\max }$. In this section, we check the dependence of the TDRHBWS results on box size Rbox, energy cutoff Ecut, truncation of expansion ${\lambda }_{\max }$ and angular momentum cutoff ${J}_{\max }$ by taking 112Ge as an example. In the calculations, the density functional PC-PK1 and the fixed gap pairing force with Δ = 12A−1/2 MeV are used.
In figure 3, the total energy and the rms radius are plotted as functions of Ecut for a ground state of 112Ge. Here Ecut is an energy cutoff parameter for positive-energy states in the Woods–Saxon basis, and the number of negative-energy states in the Dirac sea is the same as that of positive-energy states above the Dirac gap [44]. Apparently, when we increase Ecut, these quantities converge well. When Ecut increases from 80 to 120 MeV, the relative changes in the total energy are smaller than 0.1%, and the relative changes in the radius are smaller than 0.01%. Therefore, It is reasonable to use Ecut = 80 MeV in the calculations. Figure 4 shows the same quantities as functions of the box size Rbox. The relative deviations between the total energy at Rbox = 20 fm and Rbox = 22 fm are about 0.1%. The box size Rbox = 20 fm also gives good accuracy for the rms radius.
Figure 3. The total energy and rms radius of the ground state of 112Ge as functions of the cutoff energy for DWS basis Ecut.The density functional PC-PK1 is used, the box size is Rbox = 20 fm, angular momentum cutoff is ${J}_{\max }=$ 19/2 ℏ and spherical harmonics expansion truncation ${\lambda }_{\max }=6$. |
Figure 4. The total energy and rms radius of ground state of 112Ge as functions of the box size Rbox. The density functional PC-PK1 is used, the cutoff energy for DWS basis Ecut = 80 MeV, angular momentum cutoff is ${J}_{\max }=$ 19/2ℏ and spherical harmonics expansion truncation ${\lambda }_{\max }=6$. |
In figure 5, the convergence of the TDRHBWS results with respect to the spherical harmonics expansion truncation ${\lambda }_{\max }$ is investigated. It is shown that, when the expansion truncation ${\lambda }_{\max }$ increases, the total energy and the rms radius both converge well. The relative difference of the total energy between calculations with ${\lambda }_{\max }=6$ and ${\lambda }_{\max }=10$ fm is smaller than 0.1%, and the relative changes of the rms radius are smaller than 0.01%.
Figure 5. The total energy and rms radius of ground state of 112Ge as functions of spherical harmonics expansion truncation ${\lambda }_{\max }$.The density functional PC-PK1 is used, the box size is Rbox = 20 fm, the cutoff energy for DWS basis Ecut = 80 MeV, and angular momentum cutoff is ${J}_{\max }=$ 19/2 ℏ. |
As shown in figure 6, the convergence of the TDRHBWS solutions with respect to angular momentum cutoff ${J}_{\max }$ is examined for 112Ge. In figure 6, the angular momentum cutoff ${J}_{\max }=$ 19/2ℏ leads to a relative accuracy of 0.1% for the total energy and 1% for the rms radius in comparison with ${J}_{\max }=$ 21/2ℏ.
Figure 6. The total energy and rms radius of ground state of 112Ge as functions of the angular momentum cutoff ${J}_{\max }$. The density functional PC-PK1 is used, the box size is Rbox = 20 fm, the cutoff energy for DWS basis Ecut = 80 MeV and angular momentum cutoff is ${J}_{\max }=$ 19/2ℏ. |
In conclusion, in the following calculations, we fix the box size at Rbox = 20 fm, and the cutoff energy for the Dirac Woods–Saxon basis at Ecut = 80 MeV. The spherical harmonics expansion truncation is ${\lambda }_{\max }=6$. The angular momentum cutoff is ${J}_{\max }=$ 19/2ℏ.
4. Results and discussion
In this section, we present results from the triaxially deformed RHB theory in the Woods–Saxon basis. We choose germanium (Ge) isotopes as examples and focus on the influence of triaxial deformation and continuum on ground-state properties and the neutron drip line.
One-nucleon and two-nucleon separation energies S2n provide information on whether a nucleus is stable against one or two nucleon emissions and thus define the nucleon drip lines. Figure 7 shows the two-neutron separation energy of Ge isotopes from TDRHBWS and RHB-HO calculations, in which the density functional PC-PK1 and fixed gap pairing force are used. The major shell number Nshell = 14 is used in the RHB-HO calculation. It can be seen that the two-neutron separation energies from TDRHBWS and RHB-HO calculations present a similar evolution trend, and a peak at A = 102 occurs in both calculations. The S2n from the TDRHBWS calculations are close to those from the RHB-HO calculations with A ≤ 110, and the deviation is obvious with A > 110. In the RHB-HO calculations, the S2n becomes negative after A = 116 and still decreases with the increase of neutron number, while S2n from TDRHBWS calculations and becomes negative for A = 120, i.e. the neutron drip-line nuclei of Ge isotopes are 114Ge in the RHB-HO and 118Ge in the TDRHBWS calculations, respectively. It is noted that the neutron drip-line nuclei is 114Ge in RCHB calculations [40], and the deformation would influence the position of the drip line [40, 86].
Figure 7. Two-neutron separation energy S2n from TDRHBWS calculations and RHB-HO calculations for Ge isotopes as functions of mass number. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
In figure 8, the ground-state quadrupole deformation β and γ from TDRHBWS calculations and RHB-HO calculations for Ge isotopes are shown, and the density functional PC-PK1 and fixed gap pairing force are used in the calculations. It can be seen that the ground-state quadrupole deformations from the TDRHBWS calculation agree with the ones from the RHB-HO calculation. For 88−100Ge, we observe the triaxial shapes in both calculations and 92−100Ge shows triaxial shapes with γ > 30°. In 102Ge, a spherical shape is predicted in both calculations. This results in an abrupt change of the two-neutron separation energy S2n at A = 102, as seen in figure 7. From 108Ge to 114Ge, quadrupole deformation β and γ increase with the neutron number in both calculations. With further increases of neutron number, the triaxial degree of freedom γ will still increase in the TDRHBWS calculation.
Figure 8. The ground-state quadrupole deformations β and γ from TDRHBWS calculations and RHB-HO calculations for Ge isotopes as functions of mass number. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
In figure 9, the root-mean-square radii from TDRHBWS and RHB-HO calculations for Ge isotopes are plotted as functions of mass number. We display neutron radii Rn, proton radii Rp, matter radii Rm and the r0A1/3 curve with r0 = 1.02 fm. Globally, the rms radii from TDRHBWS and RHB-HO calculations present similar evolution trends. The proton radius Rp is almost constant with a very slow increase with increasing mass number. With the increasing mass number, the neutron radii Rn and matter radii Rm increase rapidly. Due to the deformation effects, there exists a clear kink at A = 102 in both calculations. For Rm, the results from these two calculations deviate from the empirical formula. The neutron, proton and matter radii from the RHB-HO calculation are smaller than those from the TDRHBWS calculation. It is noted that, for the nuclei close to the neutron drip line, the deviation of neutron radius between these two calculations increases with the increasing neutron number.
Figure 9. The neutron, proton and matter rms radii from TDRHBWS calculations and RHB-HO calculations for Ge isotopes as functions of mass number. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. The dashed line represents empirical formula Rm = r0A1/3 (r0 = 1.02 fm). |
In figure 10, the spherical component of neutron density distribution for Ge isotopes is plotted as a function of mass number. It can be seen that it is a global trend that the neutron density distributions are extended further with neutron number. The surface expands outward rapidly, and the internal density distribution changes slightly. For the nuclei near to the drip line, the tails of the density distributions extend outward and the diffuseness increases significantly.
Figure 10. The spherical component of neutron density distribution from TDRHBWS calculations for Ge isotopes as functions of mass number. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
In figure 11, the neutron single-particle levels of 112−116Ge from the TDRHBWS calculation are shown as functions of mass number. In addition, the results from the RHB-HO calculation are given for comparison. It can be seen that for the deeply bound levels (ϵcan < −2 MeV), the TDRHBWS calculation is similar to the RHB-HO calculation. However, the TDRHBWS calculation deviates from the RHB-HO calculation for levels around the Fermi surface, especially for the continuum. In the RHB-HO calculation, the number of continuums near the threshold is small and the corresponding energy is high, while a larger number of continuums occur near the threshold and the corresponding energy is closer to the threshold in the TDRHBWS calculation. Compared to the RHB-HO calculation, the TDRHBWS calculation which treats the continuum properly by solving the RHB equation with a DWS basis gives more continuums with lower energy near the threshold, thus making the Fermi surface negative and keeping the nucleus bound.
Figure 11. The neutron single-particle energies of 112−116Ge in TDRHBWS calculation(a) and RHB-HO calculation(b). The black lines represent the levels with positive parity, and the red line represent the levels with negative parity. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
To investigate the influence of triaxial deformation, we compare the bulk properties of TDRHBWS calculation and DRHBc calculation for 114Ge and 116Ge. In table. 4, the total energy Etot, neutron pairing energy ${E}_{\mathrm{pair}}^{n}$, proton pairing energy ${E}_{\mathrm{pair}}^{n}$, neutron Fermi surface λn, proton Fermi surface λp, neutron rms radius Rn, proton rms radius Rp, matter rms radius Rm and charge radius Rc from these two calculations are listed. It can be seen that the deformation β are very close in TDRHBWS and DRHBc calculations for 114Ge and 116Ge, while the triaxial deformations γ are 9.3957° and 12.8506° in TDRHBWS calculation for 114Ge and 116Ge, respectively. For the rms radii, the results from the TDRHBWS calculation are close to the ones from the DRHBc calculation, and the difference between these two calculations is about 10−3 fm. The reason is that the rms radii only depend on the spherical part of the expansion of density distribution in both calculations. Compared to DRHBc calculations, TDRHBWS calculations give larger total energy and pairing energy after considering the triaxial deformation. That is because considering γ deformation would bring more correlations and influence the single-particle levels.
Table 4. The bulk properties of 114Ge and 116Ge in TRDHBWS calculation and DRHBc calculation, including total energy Etot, neutron pairing energy ${E}_{\mathrm{pair}}^{n}$, proton pairing energy ${E}_{\mathrm{pair}}^{n}$, neutron Fermi surface λn, proton Fermi surface λp, neutron rms radius Rn, proton rms radius Rp, matter rms radius Rm and charge radius Rc. The density functional PC-PK1 and fixed gap pairing force are used in the calculations. |
| 114Ge | 116Ge | |||||
|---|---|---|---|---|---|---|
| DRHBc | TDRHBWS | δ | DRHBc | TDRHBWS | δ | |
| β | 0.2665 | 0.2709 | −0.0044 | 0.2877 | 0.2798 | 0.0079 |
| γ (deg) | 9.3957 | 12.8506 | ||||
| Rn (fm) | 5.3167 | 5.3199 | −0.0032 | 5.3859 | 5.3820 | 0.0040 |
| Rp (fm) | 4.3948 | 4.3954 | −0.0006 | 4.4112 | 4.4125 | −0.0013 |
| Rm (fm) | 5.0748 | 5.0774 | −0.0026 | 5.1356 | 5.1328 | 0.0028 |
| Rc(fm) | 4.4670 | 4.4676 | −0.0006 | 4.4832 | 4.4844 | −0.0012 |
| λn (MeV) | −0.2327 | −0.2769 | 0.0442 | −0.1060 | −0.0961 | −0.0099 |
| λp (MeV) | −23.2615 | −23.2113 | −0.0502 | −23.6612 | −23.6844 | 0.0232 |
| ${E}_{\mathrm{pair}}^{n}$ (MeV) | −5.4860 | −6.7704 | 1.2844 | −5.5614 | −6.8783 | 1.3169 |
| ${E}_{\mathrm{pair}}^{p}$ (MeV) | −1.9080 | −2.4588 | 0.5508 | −1.9085 | −2.3700 | 0.4615 |
| Etot (MeV) | −766.0057 | −768.2732 | 2.2675 | −766.3153 | −768.4260 | 2.1107 |
In figure 12, we display the self-consistent TDRHBWS triaxial quadrupole binding energy maps of 110−118Ge in the β–γ plane (0 ≤ γ ≤ 60°), obtained by imposing constraints on the expectation values of the quadrupole moments. All energies are normalized with respect to the binding energy of the absolute minimum. It can be seen that, the ground-state deformation β increases from 0.20 to 0.28 with mass number increasing, and the ground-state deformation γ increases from 4.5° to 15.2° at the same time. Furthermore, the softness in γ direction of the potential surface gradually grows with the increasing mass number.
Figure 12. The potential energy surfaces calculated by TDRHBWS for 110−118Ge. All energies are normalized with respect to the binding energy of the absolute minimum (large black dot). The energy difference between the neighboring contour lines is 0.2 MeV. |
5. Summary
In summary, the triaxially deformed relativistic Hartree–Bogoliubov theory in the Woods–Saxon basis is developed which could treat the triaxial deformation, pairing correlations and continuum. In order to consider triaxial deformation, the deformed potentials are expanded in terms of spherical harmonic functions in the coordinate space. The triaxially deformed RHB equations are solved in a Dirac Woods–Saxon basis where the radial wave functions have a proper asymptotic behavior at a large distance from the nuclear center, so as to treat pairing correlations and continuum properly. The formalism of the triaxially deformed RHB theory is presented. Taking the axially deformed nuclei 24Ne as an example, the comparison between the TDRHBWS calculation and DRHBc calculation is performed. Taking 76Ge as an example, the self-consistency of the code is confirmed by comparing the constraint calculation with γ = 30° to the constraint calculation with γ = 150°. Taking the triaxially deformed nuclei 76Ge as an example, the comparison between the TDRHBWS calculation and the triaxially deformed relativistic Hartree–Bogoliubov on the RHB-HO calculation is performed. The converge checks for the size of the Woods–Saxon basis, the box size, the spherical harmonics expansion truncation, and the angular momentum cutoff are performed.
Using TDRHBWS theory with PC-PK1 and fixed gap Δ = 12A−1/2 MeV, the nuclei far from β stability line 88−118Ge are investigated. The evolutions of the two-neutron separation energy, deformation and root-mean-square (rms) radii with mass number are analyzed. Compared to the RHB-HO theory, the TDRHBWS theory extends the neutron drip line from 114Ge to 118Ge. The canonical single-particle levels of 112−116Ge from TDRHBWS calculations are analyzed. For the single-particle levels with energy lower than −2 MeV, the results of the TDRHBWS theory are very similar to the RHB-HO theory. However, for the single-particle levels with energy higher than −2 MeV, the results are quite different, especially for the continuum. In the TDRHBWS theory, there are lots of continuums around or above the Fermi surface, and the energies of the continuum are very closed zero, while in the RHB-HO theory, the continuum around the Fermi surface are quite few and have higher energies. Taking 114Ge and 116Ge as examples, the total energy and radii from DRHBc and TDRHBWS calculations are compared. It is found that the total energy from TDRHBWS calculations is larger than the results from DRHBc, due to considering triaxial deformation leads to more correlations.