Relativistic Hartree-Fock model for axial-symmetric nuclei with quadruple and octupole deformations

Over the past decades, octupole deformation in nuclei has been an active research field in nuclear physics [1, 2]. Researchers widely acknowledge that the shape of the ground state for most deformed nuclei is symmetric under space inversion, and therefore, the dominant intrinsic deformation has a quadruple character. However, some nuclear systems exist in which the reflection symmetry is not preserved, leading to an octupole deformation or pear-like shape. Microscopically, the occurrence of the octupole deformation is attributed to the coupling between two single-particle orbits near the Fermi surface, which differ by $ \Delta l $ = 3 and $ \Delta j $= 3. Such nuclei are located close to proton and neutron numbers of 34, 56, and 88 and the neutron number of 134 [1].

In recent years, the occurrences of octupole deformation have been confirmed in various regions of the nuclear chart, such as $ ^{224} {\rm{Ra}}$ [3], $ ^{144} {\rm{Ba}}$ [4], $ ^{146} {\rm{Ba}}$ [5], $ ^{228} {\rm{Th}}$ [6], and $ ^{96} {\rm{Zr}}$ [7]. More interestingly, several nuclear novel phenomena related to the octupole deformation have also been discovered, such as the alternating-parity rotational bands, low-lying $ 1^- $ and $ 3^- $ states, and $ E3 $ transitions [1]. In parallel, considerable theoretical efforts have been devoted to understanding these colorful phenomena, including macroscopic-microscopic models [8−10], self-consistent mean-field models based on nuclear density functional theory [11−14], the interacting boson model [15−19], geometrical collective models [20, 21], and cluster models [22, 23].

As a representative, the relativistic mean field (RMF) theory based on the meson exchange diagram of nuclear force [24], also referred as the covariant density functional theory (CDFT), has achieved significant success in exploring the structural properties of nuclei that spread over almost the entire nuclear chart [25−31]. In particular, owing to the covariant representation, the RMF theory can provide a natural interpretation on the strong spin-orbit couplings [32, 33] and the origin of the pseudo-spin symmetry [34, 35]. However, in the RMF theory, the Fock terms, the inseparable parts of the meson exchange diagram of nuclear force, are excluded for simplicity. Thus, the important degrees of freedom associated with the π- and ρ-tensor couplings, which work only through the Fock terms, are missing within the RMF scheme.

Over the past decades, the relativistic Hartree-Fock (RHF) theory and its extension, which incorporate the density-dependent meson-nucleon coupling strengths [36−38], have achieved comparable accuracy as popular mean-field models in describing various nuclear structure properties, by the RHF Lagrangians PKOi(i = 1,2,3) [36, 39] and PKA1 [37]. Because of the Fock terms, notable improvements have been achieved in the self-consistent description of nuclear structure properties, such as shell evolutions [39−41], symmetry energy [42, 43], new magicity [44, 45], pseudo-spin symmetry [46, 47], and spin and isospin excitations [48, 49]. In particular, the tensor force, an important ingredient of nuclear force, has been considered naturally via the Fock terms [50, 51], e.g., using the π- and ρ-tensor couplings. Without being limited to the RHF scheme, studies have revealed that the tensor force can play an important role in nuclear shell evolution [39, 41], symmetry energy [50, 52], and nuclear excitations [53, 54].

Recently, using the Dirac Woods-Saxon (DWS) basis [55], both the RHF and relativistic Hartree-Fock-Bogoliubov (RHFB) theories [36−38] were extended for quadruple-deformed nuclei with axial symmetry, giving the D-RHF and D-RHFB models [56, 57], respectively. Note that the Dirac equations become integro partial-differential ones, which are solved by expanding the Dirac spinors on the spherical DWS basis [56, 57]. The DWS basis [55] provides reasonable asymptotic behaviors for the wave functions of unstable nuclei and microscopic insights into deformed nuclei. Taking $ ^{20} {\rm{Ne}}$ as an example, the tensor force carried by the π-coupling was demonstrated to have an essential role in determining the deformed single-particle structures [56]. Furthermore, the D-RHFB model with PKA1 reproduces both the even-parity ground state and halo structure of $ ^{11} {\rm{Be}}$ well, from which a microscopic picture of deformed halo was indicated in terms of the DWS basis [58].

For a long time, the attempt to understand the occurrence of stable octupole deformation in the nuclear chart, in particular for the role of nuclear tensor force, has been an interesting topic. As encouraged by the successes of both D-RHF and D-RHFB models, we are motivated to extend the RHF model by incorporating both quadruple and octupole deformations, resulting in the O-RHF model in this study. Similarly, the Dirac equation will be solved by expanding the Dirac spinor on the spherical DWS basis because of its advantages. Owing to the complicated Fock terms, the axial symmetry is still imposed for octupole-deformed nuclei. Meanwhile, the pairing correlations are treated within the BCS scheme by using the central part of the finite-range Gogny force D1S [59] as the pairing force.

The remainder of this paper is organized as follows. In Sec. II, the general formalism of the O-RHF model is presented using the spherical DWS basis. Section III presents the results and discussions, including the space truncation, convergence check, and description of the nucleus $ ^{144} {\rm{Ba}}$, where the effects of the Fock terms and tensor force are emphasized. Finally, a summary is given in Sec. IV.

This section provides a brief introduction of the general formalism of the RHF theory. To provide readers with a comprehensive understanding, we also present some details of the RHF energy functional and Dirac equations for the nuclei with the reflection asymmetry and axial symmetry using the spherical DWS basis.

Under the meson-exchange diagram of nuclear force, the Lagrangian of nuclear systems, which is the starting point of the RHF theory, can be obtained by considering the degrees of freedom associated with the nucleon (ψ), isoscalar σ- and ω-mesons, isovector ρ- and π-mesons, and photon field ($ A^\mu $). Starting from the general Lagrangian as detailed in Ref. [56], the Hamiltonian can be derived via the Legendre transformation as

where the kinetic energy (T) and potential energy ($ V_\phi $) terms are expressed as

Here, x and $ {\boldsymbol{x}}' $ denote the space coordinate vectors, and ψ is the Dirac spinor of a nucleon. In the potential energy terms $ V_\phi $, ϕ represents the two-body interaction channels, including the scalar (σ-S), vector (ω-V, ρ-V and A-V), tensor (ρ-T), vector-tensor (ρ-VT), and pseudo-vector (π-PV) couplings. Accordingly, the interaction vertices $ \Gamma_{\phi}(x,x') $ have the following forms [56, 60]:

In the above expressions, the symbol $ \vec{\tau} $ represents the isospin vector with $ \tau_3 $ for the projection, $ x = (t, {\boldsymbol{x}}) $, with the bold type representing space vectors, and M and $ m_\phi $ are the masses of the nucleon and mesons, respectively. After neglecting the retardation effects, that is, ignoring the time component of the four-momentum carried by the mesons and photon, the propagators $ D_{\phi}({\boldsymbol{x}}, {\boldsymbol{x}}') $ in the potential terms $ V_\phi $ are expressed as

where ϕ represents the σ-, ω-, ρ-, and π-meson fields.

Restricted on the level of the mean field approach, the no-sea approximation is introduced as usual, which is amount when the contributions from the negative energy states are neglected [60]. Thus, the nucleon field operator ψ in the Hamiltonian can be quantized as

where $ x = (t,{\boldsymbol{x}}) $, the index α denotes the positive-energy solutions of Dirac equation, and $ \psi_\alpha({\boldsymbol{x}}) $, $ \varepsilon_\alpha $, and $ c_\alpha $ ($ c_\alpha^\dagger $) represent the single-particle wave functions, energies, and related annihilation (creation) operators, respectively. Subsequently, the nuclear energy functional E can be obtained from the expectation of the Hamiltonian with respect to the Hartree-Fock ground state $ \left|{\Phi_0}\right> $ [56, 60]:

where $ \left|{-}\right> $represents vacuum, and A is the nuclear mass number. For the two-body interaction $ V_\phi $, the above expectation leads to two types of contributions, namely, the direct Hartree and exchange Fock terms. If only the Hartree terms are considered, this leads to the so-called RMF theory. If both Hartree and Fock terms are considered, it yields the RHF theory.

Taking the variational of the energy functional E, an integro-differential Dirac equation can be obtained as

where the single-particle Hamiltonian h consists of three parts, i.e., the kinetic term $ h^{\text{kin}} $, local potential term $ h^D $, and non-local terms $ h^E $ contributed by the Fock terms. The details are available in Ref. [56]. Note that the integral Dirac equation is difficult to solve directly in coordinate space, particularly when considering nuclear deformations.

In the popular RHF Lagrangians PKOi ($ i = 1,2,3) $ and PKA1, the density dependencies are introduced for the meson-nucleon coupling strengths $ g_\phi $ ($ \phi = \sigma, \omega^\mu, \vec{\rho}^\mu $) and $ f_{\phi'} $ ($ \phi' = \vec{\rho}^\mu, \vec{\pi} $) in Eq. (4), which are taken as functions of nucleon density $ \rho_b = \bar\psi\gamma^0\psi $ [36, 37]. Note that the density dependencies of the coupling strengths lead to an additional contribution to the self-energy $ \Sigma_0 $, i.e., the rearrangement terms $ \Sigma_R $ [36−38], which should be considered to preserve the energy-momentum conservation [61]. Together with the Fock terms, the density dependencies of the coupling strengths remarkably increase the numerical complexity when the deformation degree of freedom is involved [56, 57].

As limited by complicated Fock terms, we restrict ourselves under the axial symmetry for the deformed nuclei with quadruple and octupole deformations. Owing to the reflection asymmetry, the parity does not remain as a good quantum number, and the projection m of the total angular momentum remains a good one. To label the single-particle states, we introduce the symbol $ \alpha = (\nu m) $ by using the index ν to denote the orbits with different energies in a given m-block. In the following, we will derive the RHF energy functional (7) and solve the integro-differential Eq. (8) by expanding the single-particle wave functions on the spherical DWS basis [55].

In expanding the wave functions, both positive and negative energy states in the spherical DWS basis are considered. This does not conflict with the no-sea approximation, which is considered in quantizing the nucleon field operator ψ [see Eq. (6)] by neglecting the contributions from the Dirac sea, e.g., to nucleon densities. Conversely, considering the negative energy states of the spherical DWS basis is demanded by the mathematical completeness. The DWS states with negative energies were proved to be essential for correctly describing nuclear matter [62] and nuclear structures [57]. In terms of the DWS basis, the expansion of the wave function of the orbit $ \alpha = (\nu m) $ is expressed as

where $ \psi_{am} $ represents the wave functions of the spherical DWS basis states denoted by $ a = (n\kappa) $, where n represents the principle number and $ \kappa = \pm (j+1/2) $ with $ j = l\mp1/2 $. For convenience, the expansion coefficients $ C_{a,\alpha} $ are set as real numbers. To provide compact expressions, the sum over the principal number n in above equation is absorbed by $ \psi_{ \nu \kappa m} $ as

where $ {{\cal{G}}}_{\alpha\kappa} = \sum\nolimits_n C_{n\kappa,\alpha} G_{n\kappa} $, $ {{\cal{F}}}_{\alpha\kappa} = \sum\nolimits_n C_{n\kappa,\alpha} F_{n\kappa} $, and $ \Omega_{\kappa m} $ (also referred as $ \Omega_{j m}^l $) is the spherical spinor [63]. Note that the sum in the expansion (9) contains the DWS state $ \left|{n\kappa m}\right> $ with both even- and odd-parity owing to the reflection asymmetry.

For the propagators (5), the following expansion in the spherical coordinate space ($ r,\vartheta,\varphi $) [56, 60, 63] is used to derive the energy functional:

where $ {\bf{\Omega}} = (\vartheta,\varphi) $, and the index $ \lambda_y $ denotes the expansion terms of the Yukawa propagator, with $ R_{\lambda_y\lambda_y} $ representing the radial part. For details, please refer to Ref. [56].

Owing to the density dependence, the coupling strengths in the interaction vertex (4) are expanded in a series of spherical harmonic functions $ Y_{\lambda_p0} $:

where $ g_\phi $ represents the coupling strengths $ g_\sigma $, $ g_\omega $, $ g_\rho $, and $ f_\pi $ in Eq. (4), and the index $ \lambda_p $ denotes the expansion terms of the density-dependent parameters in the effective Lagrangians. In contrast to the D-RHF and D-RHFB models [56, 57], which preserve the reflection symmetry, here $ \lambda_p $ does not remain as only even integers owing to the octupole deformation.

In this study, we focus on the RHF Lagrangians PKOi ($ i = 1,2,3 $) [36, 39]. PKO2, which shares the same degrees of freedom as the popular RMF models, contains the σ-S, ω-V, ρ-V, and A-V couplings. In addition, PKO1 and PKO3 consider the π-PV coupling. Thus, the energy function E can be expressed as

where the kinetic energy $ E^{\text{kin}} = \left<{\Phi_0}\right| T\left|{\Phi_0}\right> $, and $ E_\phi^D $ and $ E_\phi^E $ are the Hartree and Fock terms of the potential energies, respectively, namely $ E_\phi = \left<{\Phi_0}\right|V_\phi\left|{\Phi_0}\right> $ with $ \phi = \sigma $-S, ω-V, ρ-V, π-PV, and A-V.

With the expansions (9), (11) and (12), we can derive the energy functional terms $ E^{\text{kin}} $, $ E_{\phi}^{D} $, and $ E_{\phi}^E $ of Eq. (13) in detail. The detailed expressions for these terms are not provided in this paper, as their forms are identical to those given in Eqs. (37), (40), and (41) of Ref. [56], with the exception of the following points. As a consequence of the reflection asymmetry, the sums with respect to $ \lambda_d $ in Eq. (40) of Ref. [56] will be released with both even and odd numbers. For the Fock terms, the restrictions on the parity index, π and $ \pi' $, in Eq. (41) and Appendix 3 of Ref. [56], are not valid for octupole-deformed nuclei. To avoid misunderstanding, the symbols $ \bar{{\mathscr D}}_{\kappa_1m_1;\kappa_2m_2}^{\lambda \bar\mu} $ and $ \bar{{\cal{Q}}}_{\kappa_1m_1;\kappa_2m_2}^{\lambda\bar\mu\sigma} $ appearing in the non-local self-energies, namely Eqs. (A5) and (A7) in Ref. [56], are redefined for octupole-deformed nuclei as

where $ \pi_1 = 0 $ and $ 1 $ for $ \kappa_1 $ represent the blocks with even- and odd-parity, respectively.

For the open-shell nuclei, the pairing correlations are considered within the BCS scheme, and the central part of the finite-range Gogny force D1S [59] is adopted as the pairing force, similar to what we did for the D-RHF model [56]. However, in contrast to the D-RHF model, where only the contribution of the main component $ J = 0 $ is included [56], the full contributions from all the J-components are considered in this study.

Because the Dirac spinors are expanded on the spherical DWS basis, the variation of the RHF energy functional (13) with respect to the expansion coefficient $ C_{a,\alpha} $ leads to a series of eigenvalue equations as

where the symbol $ \widehat C_\alpha $ is a column matrix composed of the coefficients $ C_{a,\alpha} $. By diagonalizing the Hamiltonian $ H_{m} $ for a given m-block, the eigenvalue, i.e., the single-particle energy $ \varepsilon_\alpha $, can be obtained, as well as the eigenvector $ \widehat C_\alpha $ for deformed single-particle orbit $ \alpha = (\nu m) $.

Similar to the single-particle Hamiltonian in Eq. (8), the matrix $ H_{m} $ in Eq. (16) consists of three parts, i.e., the kinetic $ H^{\text{kin}} $, local $ H^D $, and non-local $ H^E $ terms:

where the subscript $ (m) $ is omitted. These terms have the same form as those in Eqs. (54), (55), and (56) of Ref. [56]. Consistently, $ \lambda_d $ in Eq. (55) of Ref. [56] is also released with both the even and odd numbers, and the restrictions on the parity index π in Eq. (56) of Ref. [56] do not apply to octupole-deformed nuclei.

To provide some standard references for future applications of the O-RHF model, we first test the space truncations of the spherical DWS basis and the cutoff on $ \lambda_p $ in Eq. (12) using $ ^{144} {\rm{Ba}}$ as an example, which is known to exhibit octupole deformation [4]. In determining the spherical DWS basis, setting the spherical box size as $ 20 $ fm with a radial mesh step of $ 0.1 $ fm is sufficiently precise, and the continuum states with both positive and negative energies are discretized within such a spherical box. Furthermore, with the determined space truncations and the cutoff of $ \lambda_p $, we analyze the structural properties of $ ^{144} {\rm{Ba}}$ using the RHF Lagrangians PKOi ($ i = 1,2,3 $) [36, 39] and the RMF Lagrangian DD-ME2 [64]. Particular effort is devoted to the quadruple and octupole deformation effects and the role of the tensor force component carried by the π-PV coupling.

In the O-RHF calculations, two independent truncations are subjected to rigorous examination. These are the cutoffs on the quantities $ (n\kappa) $ of the spherical DWS basis [Eq. (9)] and the expansion term $ \lambda_p $ of the density-dependent coupling strengths (12). Note that the expansion terms $ \lambda_y $ in the propagators (11) are truncated automatically by the selected DWS basis states and expansion terms $ \lambda_p $ of the coupling strengths, as detailed in Eq. (A2) of Ref. [56]. Owing to the reflection asymmetry, both odd and even $ \lambda_p $ terms must be incorporated in expanding the coupling strengths (12). This differs from the D-RHF and D-RHFB models [56, 57], in which only even $ \lambda_p $ terms are required for quadruple deformed nuclei with the reflection symmetry. Generally, it is sufficient to set $ \lambda_p = 0, 1, 2, \cdots, 8 $ for the majority of nuclei.

Practically, the maximum value of m, designated as $ m_{\max} $, depends on the nucleus under consideration. Let us set $ K_m $ to denote the number of κ-blocks included in expanding the Dirac spinors with $ m_{\max} $. Consistently, $ m_{\max} $ and $ K_m $ together determine the maximum value of $ |\kappa| $, i.e., $ \kappa_{\max} = m_{\max} + K_{m} - 1/2 $. Thus, for an arbitrary Dirac spinor $ \psi_{\nu m} $, the κ-quantities in the expansion (9) are expressed as $ \kappa = \pm (m + 1/2), \pm(m + 3/2), \cdots, \pm\kappa_{\max} $, including the DWS states with both even- and odd-parity owing to the reflection asymmetry. In principle, the $ m_{\text{max}} $ value can be determined by referring to the conventional shell-model picture. For $ ^{144}_{\ 88} {\rm{Ba}}$$ _{56} $, the $ m_{\max} $ value is determined as $ 15/2 $ for both neutron and proton orbits according to the magic number $ 184 $ as predicted in the conventional shell-model picture. For such a large $ m_{\max} $ value, setting $ K_{m} = 4 $ is sufficiently accurate. As a comment, for cases with fairly large deformation, some high-j orbits may penetrate the well-known major shells, and a rigorous test shall be conducted, e.g., by increasing the $ m_{\max} $ value.

Analogous to the D-RHFB model [57], the maximum values of the principal number n for each κ-block in the expansion (9) are determined by the energy cutoff $ E_{\pm}^C $. This cutoff is defined by the single-particle energy ε of the spherical DWS state, with the sign $ + $ or $ - $indicating the positive or negative energy states, respectively. Specifically, the states with positive/negative energies, that is, $ E_{+}^C+ M>\varepsilon> E_{-}^C-M $, are considered in the expansion (9). Caution must be taken when testing the energy cutoff $ E_\pm^C $.

As mentioned earlier, the Fock terms significantly increase the numerical complexity. Thus, the O-RHF calculations become extremely time consuming. Similar to developing the D-RHF and D-RHFB models, we apply a GPU parallel speedup technology to calculate the nonlocal Fock mean fields and pairing matrix elements. This notably decreases the computing time, which makes the extension to superheavy region expectable. However, for a single calculation of the nuclide $ ^{144} {\rm{Ba}}$ with $ m_{\max} = 15/2 $, approximately 19.1 and 51.9 h are still required for PKO2 and PKO3, respectively, when fully using eight GPU units (Tesla A100-40G) and 80 CPU cores. When the $ m_{\max} $ value increases, the time cost increases further. Fortunately, it does not increase significantly when the $ E_\pm^C $ values are increased.

Figures 1(a), 1(b) and 1(c), 1(d) present the tests of the energy cutoff $ E_\pm^C $ for the total energy E (MeV) and deformation $ (\beta_2,\beta_3) $ of $ ^{144} {\rm{Ba}}$, respectively. The results were obtained by employing the RHF Lagrangians PKO2 and PKO3 and the RMF Lagrangian DD-ME2. Figures 1(a) and 1(c) depict the convergence with respect to $ E_{+}^C $, with $ E_{-}^C $ set as $ 0 $ MeV, and Figs. 1(b) and 1(d) demonstrate the convergence with respect to $ E_{-}^C $, in which $ E_{+}^C $ is fixed as $ 350 $ MeV. Both the total energy E [Fig. 1(a)] and deformation $ (\beta_2, \beta_3) $ [Fig. 1(c)] converge when $ E_{+}^C > 200 $ MeV for $ ^{144} {\rm{Ba}}$. When the $ E_-^C $ value is changed from zero to $ -50 $ MeV, namely, considering the continuum negative energy states of the DWS basis in the expansion (9), the total energy E remains almost unchanged for DD-ME2 and PKO3, whereas PKO2 exhibits a slight but visual change, as shown in Fig. 1(b).

Figure 1.  (color online) Total energy E (MeV) [plots a and b] and deformations $(\beta_2, \beta_3)$ [plots c and d] for $^{144}$Ba with respect to the positive ($+$) and negative ($-$) energy cutoffs $E_{\pm}^C$ (MeV) in expanding the Dirac spinors $\psi_{\nu m}$. The results are calculated using DD-ME2, PKO2, and PKO3.

In contrast to the total energy E, the deformations $ (\beta_2, \beta_3) $ appear to be more sensitive to the negative energy cutoff $ E_-^C $. As shown in Figs. 1(c) to 1(d), an extended truncation of the negative energy states of the DWS basis is required to give converged deformations, particularly for the octupole deformation $ \beta_3 $ given by the RHF Lagrangians PKO2 and PKO3. Despite these evident alterations, both deformations $ \beta_2 $ and $ \beta_3 $ converge rapidly with respect to the negative $ E_{-}^C $ values, namely $ E_-^C<-100 $ MeV.

To understand the convergence of the octupole deformation $ \beta_3 $, as shown in Figs. 1(c) and 1(d), Fig. 2 shows the proton single-particle spectra of $ ^{144} {\rm{Ba}}$ calculated using DD-ME2 and PKO2 with $ E_-^C = 0 $ and $ -50 $ MeV, where the Fermi levels are denoted by $ E_{\text{F}} $ and $ m_\nu $ for the deformed single-particle orbits. The results for the RHF Lagrangian PKO3 are not shown becuse they have a similar description as those of PKO2.

Figure 2.  (color online) Proton single-particle spectra of $^{144}$Ba given by DD-ME2 (plot a) and PKO2 (plot b) with negative energy cutoff $E_{-}^C = 0$ MeV and $-50$ MeV, in which the positive energy cutoff is fixed as $E_{+}^C = 350$ MeV. The ultra-thick bars represent the occupation probabilities of the orbits $m_\nu$, with the index ν representing the νth state in m-block and $E_F$ denoting the Fermi levels.

Consistent with relatively small alterations from Figs. 1(c) to 1(d), the proton single-particle spectra given by DD-ME2 remain almost unchanged when the $ E_-^C $ value varies from zero to $ -50 $ MeV [Fig. 2(a)]. Conversely, as shown in Fig. 2(b), the PKO2 results exhibit remarkable alterations. Similar systematics were observed for the neutron single-particle spectra, which are not shown. In other words, the results given by DD-ME2 remain almost unchanged, whereas the PKO2 results manifest rather distinct alterations when the $ E_-^C $ value changes from zero to $ -50 $ MeV.

In this study, the deformed single-particle orbits are expanded on the spherical DWS basis. This provides us an insight into the microscopic properties of the octupole nucleus $ ^{144} {\rm{Ba}}$, as well as the noticeable alterations in Fig. 2(b). Table 1 lists the proportions (%) of the main DWS waves for the proton valence orbits $ 3/2_8 $ and $ 1/2_{13} $ of $ ^{144} {\rm{Ba}}$. In contrast to the DD-ME2 results, the proportions of the main spherical DWS waves given by PKO2 undergo a notable redistribution from $ E_-^C = 0 $ to $ E_-^C = -50 $ MeV, particularly for the $ 1h_{11/2} $ and $ 2d_{5/2} $ waves, which are essential for the occurrence of the octupole deformation. This indicates that the completeness regarding the DWS basis states with negative energies can be essential for correctly describing the octupole deformation.

orbits	$E_{-}^C$	$2d_{5/2}$	$1g_{7/2}$	$1g_{9/2}$	$1f_{5/2}$	$1h_{11/2}$
$3/2_{8}$	0	5.9	65.5	5.4	2.4	11.5
	−50	6.6	59.5	2.6	6.9	12.6
$1/2_{13}$	0	21.8	24.4	7.7	0.7	27.1
	−50	25.8	20.1	2.6	2.1	30.1
$3/2_{8}$	0	2.1	43.8	25.3	4.0	2.4
	−50	6.5	53.8	6.2	3.5	20.2
$1/2_{13}$	0	2.9	23.5	37.7	3.3	5.0
	−50	24.7	2.8	8.8	1.1	30.1

Table 1.  Proportions (%) of the main expansion components of proton orbits $3/2_8$ and $1/2_{13}$. These results are calculated with DD-ME2 (upper panel) and PKO2 (lower panel) by selecting negative energy cutoffs $E_{-}^C$ of 0 and –50 MeV, respectively, and a positive energy cutoff of $E_{+}^C = 350$ MeV.

To ensure the completeness of the expansion (9), we must consider the spherical DWS states with negative energies, although the relevant contributions are rather small. As emphasized earlier, the mixing of the high-j waves such as $ h_{11/2} $is essential for obtaining the octupole deformation. In contrast to the low-j waves, fewer high-j waves with negative energies are bound in the DWS basis, which corresponds to $ E_-^C = 0 $ MeV. Thus, when the high-j waves become essential, considering only bound negative-energy states ($ E_-^C = 0 $) to guarantee the completeness of the expansion (9) is insufficient. Therefore, the notable alterations of the deformation ($ \beta_2, \beta_3 $) from Figs. 1(c) to 1(d) are an indication of the improved completeness. In particular, the Fock terms involve more two-body correlations, in contrast to the Hartree terms [57]. Consequently, the RHF calculations with PKO2 and PKO3 appear to be more sensitive to the negative energy cutoff $ E_-^C $ than the RMF one, as illustrated in Figs. 1 and 2.

With the given space truncations, namely $ E_+^C = 350 $ MeV, $ E_-^C = -150 $ MeV, $ m_{\max} = 15/2 $ and $ K_{m} = 4 $, we perform the O-RHF calculations for the octupole nucleus $ ^{144} {\rm{Ba}}$, using the RHF Lagrangians PKOi ($ i = 1,2,3 $) [36, 39] and the RMF one DD-ME2 [64]. Table 2 lists the total energy E (MeV), deformations $ (\beta_2, \beta_3) $ , charge radius $ r_c $ (fm) for the ground state (g.s.), quadruple local minimum (q.m.), and scenario with the spherical shape (sph.), respectively.

	Cases	E	$(\beta_2, \beta_3)$	$r_c$	$\Delta E$
Exp.		$-1190.22$	(0.18, 0.17)	4.9236
DD-ME2	g.s.	−1187.95	(0.23, 0.09)	5.0340	−1.51
	q.m.	−1186.44	(0.20, 0.00)	5.0013	−5.77
	sph.	−1180.87	(0.00, 0.00)	4.9013
PKO2	g.s.	−1189.47	(0.22, 0.14)	5.0234	−3.69
	q.m.	−1185.78	(0.19, 0.00)	4.9971	−3.28
	sph.	−1182.50	(0.00, 0.00)	4.9828
PKO1	g.s.	−1192.54	(0.24, 0.14)	5.0525	−3.15
	q.m.	−1189.39	(0.21, 0.00)	5.0224	−4.27
	sph.	−1185.12	(0.00, 0.00)	5.0047
PKO3	g.s.	−1191.77	(0.23, 0.13)	5.0365	−2.69
	q.m.	−1189.08	(0.21, 0.00)	5.0066	−4.61
	sph.	−1184.47	(0.00, 0.00)	4.9878

Table 2.  Total energy E (MeV), quadruple and octupole deformations $(\beta_2, \beta_3)$, and charge radius $r_c$ (fm) of $^{144}$Ba, calculated using PKOi ($i = 1,2,3$) and DD-ME2 for the octupole g.s., quadruple local minima (q.m.), and spherical shape (sph.). The experimental data from Refs. [4, 65, 66] are shown as references. The last column $\Delta E$ (MeV) shows the energy differences between the g.s. and q.m. cases and those between the q.m. and sph. cases, respectively.

As shown in Table 2, $ ^{144} {\rm{Ba}}$ is predicted to be deeper bound from the sph. to the q.m. cases, and the total energy E becomes more closely aligned with the experimental value [65]. Further considering the octupole deformation, all the selected Lagrangians yield the octupole g.s. for $ ^{144} {\rm{Ba}}$, and PKO2 exhibits the best agreement with the experimental data, including the total energy E, deformations $ (\beta_2, \beta_3) $, and charge radius $ r_c $. Note that a value of $ \beta_3 = 0.17(_{-6}^{+4}) $ was derived experimentally with $ \beta_2 = 0.18 $ [4]. As shown in Table 2, the $ \beta_3 $ value given by the RMF model DD-ME2 is even less than the lower limit $ 0.11 $, similar to previous mean-field calculations [8, 67−69]. In contrast, the values of $ \beta_3 $ given by the RHF models PKOi $ (i = 1,2,3) $ are all within the uncertainty of the experimental value, which deserves to be explored with further detail.

As shown in Table 2, the quadruple deformations $ \beta_2 $ for the g.s. slightly increase from those in the q.m. cases. Thus, the energy differences between the g.s. and q.m. cases and between the q.m. and sph. cases, namely, the $ \Delta E $ values (MeV) in the last column of Table 2, can be approximately considered as the effects of the octupole and quadruple deformations, respectively. We observe that the RHF Lagrangians PKOi ($ i = 1,2,3 $) produce more gains in binding caused by the octupole deformation than the RMF Lagrangian DD-ME2. Generally, this can be attributed to the effects of the Fock terms, which consider more two-body correlations than the Hartree terms. Furthermore, as compared with the value of $ \sim 0.5 $ MeV given by the point coupling RMF model PC-PK1 [13] and $ \sim1.8 $ MeV given by the projected HFB calculations [12], the selected RHF models also predict more pronounced effects of octupole deformation.

Moreover, it is somewhat beyond our expectation that PKO1 and PKO3, which contain the degree of freedom associated with the π-PV coupling, yield weaker octupole enhancements than PKO2. We must mention that PKO2 does not contain the π-PV coupling, and PKO3 carries stronger π-PV coupling than PKO1. We observe that the $ \Delta E $ values between the g.s. and the q.m. cases decrease in sequence from PKO2 to PKO1 and further to PKO3, and vice versa for the values of $ \Delta E $ between the q.m. and sph. cases. This suggests that the π-coupling does not favor the octupole deformation but rather enhances the quadruple deformation. Indeed, PKO1 and PKO3 have been shown to produce an enhanced quadruple deformation effect in comparison with PKO2. For instance, in the well-deformed nucleus $ ^{20} {\rm{Ne}}$, the tensor force component carried by the π-PV coupling has been demonstrated to enhance the deformation effects [56]. Moreover, the ρ-tensor coupling, owing to the nature of Lorentz tensor coupling [57], can enhance the deformation effects, as evidenced in $ ^{24} {\rm{Mg}}$ [57] and $ ^{11} {\rm{Be}}$ [58, 70]. This may be applicable for nuclei with octupole deformation, which calls for the future implementation of the ρ-tensor coupling in the O-RHF model. A qualitative understanding of the systematics of octupole and quadruple enhancements can be obtained by combining these with the single-particle spectra of $ ^{144} {\rm{Ba}}$.

Figures 3(a) and 3(b) show the neutron and proton spectra of $ ^{144} {\rm{Ba}}$ given by PKO2, respectively. Note that the single-particle spectra given by PKO2 are not significantly different from the ones given by PKO1 and PKO3, which are not shown here. As shown in Fig. 3, from q.m. to the octupole g.s., both neutron and proton valence orbits (marked in blue color) become much deeper bound, and the induced large shell gaps (denoted by double arrows) stabilize the octupole deformation for the g.s. of $ ^{144} {\rm{Ba}}$. To further understand the effects of octupole deformation, Table 3 shows the proportions (%) of the main DWS waves for the neutron and proton valence orbits. The intrusions of both high-lying neutron $ 1i_{13/2} $ and proton $ 1h_{11/2} $ waves play a key role in producing the octupole deformation, which is commonly supported by the selected models in this study.

Figure 3.  (color online) Neutron (left) and proton (right) spectra of $^{144}$Ba for the sph., q.m., and octupole g.s. given by PKO2. The ultrathick bars represent the occupation probabilities of the orbits, and $m_\nu$ and $E_F$ denote the Fermi levels.

Neutron		$1i_{13/2}$	$3p_{3/2}$	$2f_{5/2}$	$2f_{7/2}$	$1h_{9/2}$	$1h_{11/2}$
$3/2_{12}$	q.m.	0.0	2.3	16.5	11.7	62.4	4.4
	g.s.	24.5	2.9	12.1	22.5	11.5	14.1
$1/2_{18}$	q.m.	0.0	23.9	0.8	42.7	18.8	8.3
	g.s.	23.9	2.4	15.7	12.4	22.0	3.1
$1/2_{17}$	q.m.	0.0	6.5	27.4	4.0	48.2	3.1
	g.s.	3.5	15.6	15.7	5.6	11.3	14.6
Proton		$1h_{11/2}$	$3s_{1/2}$	$2d_{3/2}$	$2d_{5/2}$	$1g_{7/2}$	$1g_{9/2}$
$3/2_{8}$	q.m.	0.0	−	4.1	3.0	86.9	3.8
	g.s.	20.2	−	1.9	6.6	53.5	5.9
$1/2_{13}$	q.m.	0.0	14.1	4.4	55.7	14.2	9.1
	g.s.	29.8	7.1	5.7	25.3	2.3	8.6
$1/2_{12}$	q.m.	0.0	2.9	18.8	2.2	70.2	2.6
	g.s.	3.4	3.3	21.4	0.0	49.6	3.3

Table 3.  Proportions (%) of the main DWS waves in the neutron orbits $3/2_{12}$, $1/2_{18}$, and $1/2_{17}$ (upper panel) and proton orbits $3/2_8$, $1/2_{13}$, and $1/2_{12}$(lower panel) of $^{144}$Ba. These results are calculated using the RHF Lagrangian PKO2.

Taking proton as an example, the couplings between the $ 2d_{5/2} $ and intruding $ 1h_{11/2} $ waves, which fulfill the conditions $ \Delta l = 3 $ and $ \Delta j = 3 $, play an essential role in giving a stable octupole deformation. At the q.m., the orbits $ (3/2_8, 1/2_{12}) $ and $ 1/2_{13} $ are dominated by the $ 1g_{7/2} $ and $ 2d_{5/2} $ waves (in bold types), respectively, as shown in the lower panel of Table 3. Consistently, when the octupole deformation is involved, more evident couplings between the $ 1h_{11/2} $ and $ 2d_{5/2} $ waves are observed for the orbit $ 1/2_{13} $ than the others $ (3/2_8, 1/2_{12}) $. This partly explains the octupole enhancements, particularly for the proton orbit $ 1/2_{13} $, which becomes significantly more deeper bound from the q.m. to the octupole g.s., as depicted in Fig. 3(b). For a neutron, the couplings between the $ 2f_{7/2} $ and intruding $ 1i_{13/2} $ waves are significant to produce the octupole deformation, as shown in the upper panel of Table 3. In particular, the proportions of both neutron and proton pseudo-spin partners, ($ 2f_{7/2}, 1h_{9/2} $) and ($ 2d_{5/2}, 1g_{7/2} $), respectively, undergo significant redistributions with the penetrations of the high-j waves $ 1i_{13/2} $ and $ 1h_{11/2} $. This can favor the emergence of octupole deformation in $ ^{144} {\rm{Ba}}$.

As demonstrated by the $ \Delta E $ values given by PKOi $ (i = 1,2,3) $ in Table 2, the π-PV coupling, which contributes only via the Fock terms, appears to have an effect against the octupole deformation but enhance the effect of the quadruple deformation for $ ^{144} {\rm{Ba}}$. Combining Fig. 3 and Table 3, we can easily understand the role of the π-PV coupling, in which the tensor force components lead to repulsive couplings between the $ j_> = l+1/2 $ ($ j_< = l-1/2 $) and $ j_> $ ($ j_< $) waves but attractive ones between the $ j_< $ and $ j_> $ waves [71]. At the q.m., as shown in Fig. 3 and Table 3, the valence neutron ($ 3/2_{12}, 1/2_{17} $) and proton $ (3/2_8, 1/2_{12}) $ orbits are dominated by the $ j_< $ waves, i.e., the $ 1h_{9/2} $ and $ 1g_{7/2} $ waves, respectively. In contrast, as the results of strong spin-orbit couplings, the core of $ ^{144} {\rm{Ba}}$ corresponds to $ N = 82 $ and $ Z = 50 $, both representing the nature of $ j_> $ on average. Thus, owing to the attractive tensor force effects carried by the π-PV coupling, PKO1 and PKO3 present a stronger quadruple enhancement than PKO2, similar to that revealed in Ref. [56]. Notice that the RMF Lagrangian DD-ME2, which do not contain either the Fock terms or the π-PV coupling, present a stronger quadruple enhancement than the RHF models PKOi ($ i = 1,2,3 $). This can be attributed to the fact that the modeling of nuclear binding, primarily the interplay between strong attractive σ-S and repulsive ω-V couplings, is changed by the presence of the Fock terms, particularly by the much enhanced isovector ρ-V coupling, compared with that in RMF models like DD-ME2 [36, 47].

As mentioned earlier, the intrusion of the $ j_> $ waves, namely neutron $ 1i_{13/2} $ and proton $ 1h_{11/2} $ waves, is essential for the emergence of the octupole deformation in $ ^{144} {\rm{Ba}}$. However, owing to the nature of tensor force, this results in repulsive tensor couplings between the core of $ ^{144} {\rm{Ba}}$ and the intruding waves, because both represent the nature of $ j_> $. Meanwhile, as shown in Table 3, the $ j_< $ waves in both neutron and proton valence orbits are averagely reduced from the q.m. to the octupole g.s., which further weakens the attractive tensor couplings between the core and $ j_< $ waves of valence nucleons. Thus, we can easily understand that the π-PV coupling, primarily owing to its tensor force components, manifests an effect against the octuple deformation for $ ^{144} {\rm{Ba}}$, as revealed by the $ \Delta E $ values in Table 2. This may be applicable for the regions with $ N/Z = 34 $, $ 56 $, and $ 88 $ and $ N = 136 $, in which the couplings $ (1g_{9/2}, 2p_{3/2}) $, $ (1h_{11/2}, 2d_{5/2}) $, $ (1i_{13/2}, 2f_{7/2}) $, and $ (1j_{15/2}, 2g_{9/2}) $ play an essential role in producing the octupole deformation, respectively. This is because both intruding waves $ (1g_{9/2}, 1h_{11/2}, 1i_{13/2}, 1j_{15/2}) $ and cores ($ N = 34, 56, 88,136 $ and/or $ Z = 34, 56, 88 $) carry the nature of $ j_> $, between which the tensor force couplings are repulsive. This also deserves a special focus in future applications of the O-RHF model.

In this study, the relativistic Hartree-Fock (RHF) theory is extended to accommodate the octupole-deformed nuclei, giving the O-RHF models. In terms of the spherical Dirac Woods-Saxon (DWS) basis, the general formalism of the O-RHF model is introduced, and the space truncations of the DWS basis are tested. Taking $ ^{144} {\rm{Ba}}$ as an example, the mechanism related to the occurrence of the octupole deformation is analyzed with a special focus on the Fock terms and the tensor force effects carried by the π-PV coupling.

We find that the selected RHF models can reproduce the octupole deformation of $ ^{144} {\rm{Ba}}$ within the uncertainty of the experimental measurements, in contrast to the other popular mean field models. Moreover, the intrusion of the high-lying neutron $ 1i_{13/2} $ and proton $ 1h_{11/2} $ waves, which results in the octupole deformation in $ ^{144} {\rm{Ba}}$, is notably enhanced by the Fock terms. This leads to an enhanced effect of octupole deformation. In particular, owing to the repulsive tensor couplings between the core of $ ^{144} {\rm{Ba}}$ and the intruding $ j_> $ waves ($ 1i_{13/2} $ and $ 1h_{11/2} $), the tensor force components carried by the π-PV coupling present an effect against the octupole deformation for $ ^{144} {\rm{Ba}}$, which deserves special attention in future research on octupole nuclei.