Anatomy of octupole correlations in $^{96}$Zr with a symmetry-restored multidimensionally-constrained covariant density functional theory

A recent analysis based on the STAR measurement in relativistic heavy-ion collision experiments provides evidence of octupole correlation in the ground state of $^{96}$Zr, the description of which presents a challenge to nuclear structure models. In this work, we perform projection-after-variation calculations for $^{96}$Zr based on a multidimensionally-constrained relativistic Hartree-Bogoliubov model. Our results show that an octupole deformed shape is favored in energy after symmetry restoration, and this phenomenon cannot be reproduced in the pure mean-field calculations. These complex structures originate from the competition among various shell structures in this mass region. Our results suggest that the allowance of symmetry breaking at the mean-field level and the restoration of broken symmetries are essential elements for understanding the structure of $^{96}$Zr in density functional theories.


INTRODUCTION
In atomic nuclei, the occurrence of spontaneous symmetry breaking results in various shapes. The most common nuclear shape is the ellipsoid. For some specific proton or neutron numbers, we can also find pronounced pear-like shapes [1][2][3][4][5][6]. For a long time, the only reliable method for studying these reflection-asymmetric shapes has been to search for the characteristic patterns of parity and angular momentum in nuclear spectroscopy. For example, the enhanced reduced electric-octupole transition probabilities, B(E3), in 224 Ra [7], 144 Ba [8], 146 Ba [9], and 228 Th [10] can be viewed as direct experimental evidences of strong octupole deformations.
Recently, it has been proposed that nuclear deformations can also be extracted from relativistic heavyion collisions (RHIC) by analysing the hydrodynamic collective flow of the final-state particles [11][12][13][14][15]. Later this method was applied to study realistic experimental ativistic Hartree-Bogoliubov (RHB) theory [5], the Hartree-Fock-Bogoliubov (HFB) theories with Gogny interaction [4] and Skyrme interactions [6], do not find any octupole deformation for 96 Zr. Instead, all these theoretical investigations except for the MM model find octupole deformations near N = Z = 40, which is incompatible with both the experiments and the argument of octupole magic numbers. To reconcile this inconsistency between the theories and experiments, we need to consider the beyond-mean-field effects not included in the systematic calculations.
The nuclear density functional theories are very successful in describing the nuclear deformations microscopically. In these calculations, the ground states are obtained by applying the variational principle with a Slater determinant or a Bogoliubov vacuum. The resulting states usually violate fundamental symmetries and conservation laws, such as rotational symmetry, translational invariance and particle-number conservation [31].
Besides the axial deformations, recent PAV calculations based on Gogny interaction predicted a tetrahedral shape associated with a pure non-axial β 32 deformation in 96 Zr [43]. Nevertheless, so far there is no systematic investigations for the evolution of octupole deformations in the A ∼ 100 region based on full angular momentum and parity projections. Such calculations will provide necessary knowledge that can help understand the recent RHIC experiments as well as the spectroscopic data.
For studying the various nuclear shapes in a unified framework, the multidimensionally-constrained relativistic Hartree-Bogoliubov (MDCRHB) model has been developed, in which both the axial and spatial reflection symmetries are broken but the simplex symmetry is preserved [52,53]. Recently we developed a projected multidimensionally-constrained relativistic Hartree-Bogoliubov (p-MDCRHB) model by incorporating the parity and angular momentum projections into the MDCRHB model to study the low-lying states related to exotic nuclear shapes, e.g., the triangular shape associated with the three-α configuration in 12 C [54].
In this paper, we investigate the octupole correlations in 96 Zr by employing the p-MDCRHB model. We present a systematic beyond-mean-field study with both triaxial and octupole shapes included simultaneously. In Sec. 2, the theoretical framework of this model is introduced. In Sec. 3, we use the p-MDCRHB model to calculate the projected PESs in 96 Zr and extract the deformations and excitation energies of the low-lying states. The results are discussed in comparison with spectroscopic measurements and relativistic collision experiments. The microscopic mechanism of the octupole correlations are analysed based on the evolution of single-particle shell structures. Finally, a summary is given in Sec. 4.

Theoretical framework
We start with the RHB theory with effective point coupling interactions [55,56]. The RHB equation in coordinate space reads [57,58] where ∆ is the pairing field, λ is the Fermi energy, E k is the quasi-particle energy and (U k (r), V k (r)) T is the quasi-particle wave function. h is the single-particle Hamiltonian, where M is the nucleon mass, S (r) is the scalar potential, and V(r) is the vector potential. In this work, we use the PC-PK1 parameter set [59] and solve the RHB equation by expanding the wave functions on an axially-deformed harmonic oscillator basis. For pairing interaction, we adopt a separable pairing force of finite range in the spin-singlet channel [60,61], where G is the pairing strength, R = (r 1 + r 2 )/2 and r = r 1 − r 2 are the center of mass and relative coordinates, respectively. The symbols with and without the prime denote the quantities before and after the interaction, respectively. P σ is the spin-exchange operator. P(r) = (4πa 2 ) −3/2 e −r 2 /4a 2 is the Gaussian function. a = 0.644 fm is the effective range of the pairing force. In this work we use different pairing strength for neutron and proton, G n = 728.00 MeV·fm 3 and G p = 815.36 MeV·fm 3 . These parameters were adjusted to reproduce the available empirical pairing gaps of 102,104 Zr [62]. Nuclear shapes are characterized by the deformation parameters β λµ which are defined as where R = 1.2A 1/3 fm with A the mass number, Q λµ is the multipole moment of the intrinsic densities. With symmetry imposed in this model, all β λµ with even µ can be considered simultaneously.
In the p-MDCRHB model [54], the nuclear wave function |Φ(q) is projected onto the parity π and angular momentum J where f JKπ α is the weight function, q represents a collection of the deformation parameters, is an operator which projects out the component with angular momentum J and its projection M from the deformed mean-field wave function. K is angular momentum projection onto z-axis in the intrinsic frame. D J MK (Ω) is the Wigner function with a Euler angle Ω ≡ {φ, θ, ψ}.R represents spatial rotation. The parity projection operator writeŝ whereP is the spatial reflection operation. Considering that the projection calculation with both rotational and parity symmetry breaking is already costly, in this first step towards full projection, we do not include the exact particle number projection. Instead, we add two correction terms to the Hamiltonian kernel as in Ref. [63] to approximately restore the average proton and neutron numbers, i.e., where λ p and λ n are the Fermi energies for protons and neutrons of the mean-field state |Φ(q) , Z(r; q, q ; Ω) and N(r; q, q ; Ω) are the mixed vector densities for protons and neutrons, respectively. Z 0 and N 0 are the target proton and neutron numbers, respectively. After calculating the mixed energy density, the weight function f JKπ α and the eigenvalue E Jπ α are obtained by solving the generalized eigenvalue equation in the standard way [64,65]: (9) where N is the norm kernel. To check the validity of the ansatz of introducing the correction terms in Eq. (8), we also carry out a beyond relativistic meanfield calculation with exact particle-number projections by employing the approach developed in Ref. [39]. We note that the restoration of broken symmetries with modern nuclear energy density functionals usually meet the problems of spurious divergencies [66,67] and finite steps [68][69][70][71]. The solution to this problem in a proper way is still an open question. Of course, these issues can be avoided in Hamiltonian-based frameworks if all the terms in the projected energy are taken into account in a consistent way, in particular Fock terms and contributions of the Coulomb and spin-orbit potentials to pairing, etc. [66,72,73]. In the current covariant density functional theories, these problems remain in one way or another. They are not serious in this work as the symmetry restoration is carried out after variation. Besides, to avoid the numerical ambiguity from the division of zero by zero, an odd number of mesh points in the integration over the gauge angles for particle number projection is adopted. More arguments have already been made in Ref. [65].

Results and discussions
In MDCRHB model and its extensions, thanks to the usage of an unified axial symmetric basis, we can alleviate the computational burden by simply disregarding the matrix elements that violate the specific symmetries. We break a symmetry if and only if it is an essential element in describing certain processes we are focusing on. Therefore, four types of spatial symmetries can be imposed and associated with some non-zero quadrupole-octupole deformations: (a) axialreflection symmetry with only β 20 ; (b) non-axial but reflection symmetry with (β 20 , β 22 ); (c) axial symmetry but reflection asymmetry with (β 20 , β 30 ); (d) non-axial and reflection asymmetry with (β 20 , β 22 , β 30 , β 32 ), which are labeled as Kπ, / Kπ, K/ π, and / K/ π, respectively. Fig. 1(a-e) shows the total energy of 96 Zr as a function of β 20 , β 22 , β 30 , β 32 and β 40 , respectively, where all other deformation parameters are constrained to be zero. The dash-dotted and solid lines correspond to the meanfield and projected 0 + PESs, respectively. It is shown in Fig. 1(a) that the mean-field energy surface is very soft around the spherical shape. Any slight changes of the pairing strength may shift the location of the energy minimum. Specifically, with G n = 728.00 MeV·fm 3 and G p = 815.36 MeV·fm 3 , we obtain two energy minima located at prolate (β 20 = 0.21) and oblate (β 20 = −0.19) respectively with energy difference 0.32 MeV. If we enhance (quench) the neutron (proton) pairing strength by 20%(10%), the PES is still soft, but the energy minimum is shifted to the spherical shape. A proper description of the ground state needs to go beyond the mean-field approximation.
For E(β 22 ) shown in Fig. 1(b), the mean-field PES is very soft. The projected 0 + energy minimum is located at β 22 = 0.17 and the corresponding binding energy is slightly overestimated. The PESs E(β 30 ) and E(β 32 ) with octupole deformations are shown in Fig. 1(c) and (d). One can see that the energy minimum in the mean-field PES is shifted from spherical to octupole deformed state with β 30 = 0.23 or β 32 = 0.16 respectively after symmetry restorations. The energies gained from projection are 9.33 MeV and 8.65 MeV, respectively, which are consistent with the results of HFB model with Gogny interaction [43]. For E(β 40 ) in Fig. 1(e), The energy minimum changes from spherical to hexadecapole shape after projections and we gain 3.9 MeV extra binding energy by symmetry restorations. We define the projection energy E proj as where E(MF) min and E(0 + ) min are the global energy minima on the mean-field and 0 + PESs under the same symmetry restriction respectively, and show them in Fig. 1(f). One can find E proj 's for octupole deformations are the dominant ones, which are larger than those for quadrupole and hexadecapole deformations. Next we turn to the two-dimensional PESs. The results are shown in Fig. 2. We also list the deformation β λµ , projection energy E proj and excitation energy E x corresponding to the global energy minima E min on each PES in Table 1. We define E x as where we compare the energy minima on different PESs which usually take different shapes. Note that this approximation only works for states with definite shapes. For soft PESs we also need the GCM to take into account the shape fluctuations or shape coexistence. In / Kπ calculations, thanks to the six-fold symmetry with pure quadrupole deformations, we constrain the Hill-Wheeler coordinates β 2 = β 2 20 + 2β 2 22 , γ = arctan( √ 2β 22 /β 20 ) with γ ∈ [0, 60 • ] instead of (β 20 , β 22 ). Note that the octupole deformations will spoil the elegant six-fold symmetry and we constrain β λµ 's directly in other cases. The mean-field PES is soft along the γ direction and the two energy minima correspond to the prolate and oblate shape isomers in Fig. 1(a), respectively. After symmetry projections, there is only one energy minimum with (β 20 , β 22 ) = (0.23, 0.07), i.e., (β 2 , γ) = (0.25, 23 • ) in Hill-Wheeler coordinates. Thus the ground state of 96 Zr is triaxially deformed, consistent with the Monte Carlo shell-model (MCSM) calculations [76]. The triaxial degree of freedom also appears in the energy minima of the 2 + and 4 + PESs. These energy minima assume similar deformations, indicating that they belong to the same collective band. The excitation energies E x (2 + ) = 0.91 MeV and E x (4 + ) = 2.31 MeV are lower than their corresponding experimental values 1.75 MeV and 2.75 MeV [75], respectively.
In Fig. 2(b) we present the mean-field and projected PESs E(β 20 , β 30 ) with the K/ π symmetries. On the K/ π symmetry imposed. The (c) K/ π+PNP(δ) case calculated by using three-dimensional harmonic oscillator expansion method and densitydependent δ pairing force developed in Ref. [39] are shown for comparison. The contours join points with the same energy, and the separation between adjacent contours is 0.5 MeV. On each PES the red star marks the corresponding global energy minimum. The values marked in (d) are the constrained (β 20 , β 30 ) in this column, which are found in the global minima of (b). In each column the energies are relative to the mean-field ground state. Table 1: The deformations β λµ and excitation (projection) energies E x (E proj ) of various J π (MF) states for 96 Zr. E(0 + ) min denote the energy of the projected 0 + states. All excitation energies are relative to the 0 + states with the same symmetry. For / K/ π symmetry, (β 20 , β 30 ) deformations are the same as those in the K/ π global energy minimum. In the last line, the deformations inferred from the RHIC experiments [20] are listed for comparison. The experimental excitation energies are taken from Ref. [75]. All energies are in MeV.  [20] 0.062 0.202 mean-field PES, the two energy minima are located in (β 20 , β 30 ) = (0.21, 0.00) and (−0.19, 0.00), respectively. After symmetry projections, we find a single energy minimum at (β 20 , β 30 ) = (0.01, 0.23) on the 0 + PES. The symmetry restorations give a large energy correction E proj = 8.53 MeV. The projected 2 + and 4 + PESs have similar structures: they both exhibit energy minimum at (β 20 , β 30 ) ≈ (0.2, 0.1). One may conjecture that these states together with another 0 + state with similar deformations belong to a typical rotational band. To examine the existence of such 0 + state, we also need further GCM calculations. On the 3 − PES in Fig. 2(b), we find a single energy minimum at (β 20 , β 30 ) = (0.03, 0.24) with E x (3 − ) = 4.90 MeV. In contrary to the 2 + and 4 + states, its shape is close to that of the 0 + state. We also present the 1 − PES. We find that the 1 − and 3 − minima have very different shapes. The very small experimental E x (3 − ) = 1.90 MeV and the deformations of the states point out the 3 − state may be a octupole vibrational collectivity [29] and the 0 + is the band head of octupole vibration instead of that of a nearly spherical band predicted by the MCSM calculations [76].
We also calculate E(β 20 , β 30 ) in K/ π case with full particle-number projection (PNP) using the method developed in Ref. [65]. The results are shown in Fig. 2(c). Note that in these calculations we used different parameter settings compared to the other columns, such as the dimension of the basis space, the pairing force (in (c) we used a δ-pairing force, while in the other columns we used a separable pairing force) and the particle-number corrections (in (c) we performed an exact particle-number projection, while in the other columns we used the approximate ansatz Eq. (8)). For all J π states with the same symmetries imposed, the PESs are qualitatively similar. The only noticeable difference is that the approximate scheme in (b) gives prolate minima for mean-field and 2 + projected PESs, while the exact PNP method in (c) finds oblate minima on these PESs. The reason is that the PES of this nucleus is γ-soft and both axial-symmetric minima are nearly degenerate. A slight difference in the parameters can reverse the order of them. Thus, for more precise calculations the GCM must be included to mix the shapes.
In Fig. 2(d) we show the energy E(β 22 , β 32 ) of meanfield state and the states with projections onto different spin-parity as a function of (β 22 , β 32 ), while (β 20 , β 30 ) are fixed to the values corresponding to energy minima in Fig. 2(b). For mean-field PES, the energy minimum is located at (β 22 , β 32 ) = (0.0, 0.0), indicating that the nonaxial deformations are not favored, consistent with the prediction of the RHB theory [5], the HFB theories with Gogny interaction [46] and Skyrme interactions [6]. In contrast, after projection onto spin-parity J π = 0 + , the energy minimum shifts to the deformations (β 20 , β 22 , β 30 , β 32 ) = (0.01, 0.06, 0.23, 0.00). However, the projected 0 + PES is soft against (β 22 , β 32 ) deformations. One may expect large shape fluctuations, the consideration of which is beyond the scope of this work. Nevertheless, our results show nontrivial contribution from both triaxial β 22 and tetrahedral β 32 deformations, which is consistent with the finding in Refs. [43,76] and the RHIC experiment [77]. Of course, this argument is based on the assumption that the location of the energy minimum in (β 20 , β 30 ) plane does not move when varying (β 22 , β 32 ). A more solid conclusion should be made based on the result of calculation for the states within the whole (β 20 , β 22 , β 30 , β 32 ) deformation space, which is too expensive to be carried out in this work. Similarly, the above argument would also apply to all excited levels.
Next we examine the calculated energy levels and transition amplitudes. In Fig. 3, the calculated lowlying excitation energies, defined in Eq. (11), are shown and compared with the experimental data [75].
The E x (2 + ) with / Kπ symmetry is lower than the experimental value but the situation is opposite with K/ π case. The excitation energies extracted from full PNP calculations are larger than the experimental values by about 2 MeV, lower than the case with particlenumber restriction for about 1 MeV. By breaking the axial symmetry, the energy levels 2 + , 4 + and 1 − in / K/ π calculations are lowered compared with the results in the K/ π calculations. Conversely, the 3 − energy is almost unaffected by the triaxial deformations. These level schemes are qualitatively similar to that from MCSM calculations [29].
In the current framework, the electric multipole transitions are calculated between states that are projected out with the same deformations. In other words, the shift of the nuclear shapes with spin is neglected. Under this approximation, we calculate the B(E3; 3 − → 0 + ) using the symmetry-projected wave functions with intrinsic deformations (β 20 , β 30 ) = (0.01, 0.23), corresponding to the energy minimum of 0 + state labelled with K/ π. The result is 20.68 × 10 3 e 2 fm 6 , which is comparable to the recent experimental value of 23(2) × 10 3 e 2 fm 6 [29]. Considering the fact that the energy minimum of the 3 − state coincidents with that of the 0 + state, as shown in Fig. 2(b), it is reasonable to compare this calculated value with the data. From the energy spectrum, the lowest 3 − state probably belongs to an octupole vibration band, the proper description of which requires the inclusion of the shape mixing. It is beyond the current framework. However, it is interesting to note that the B(E3) is in surprisingly good agreement with experimental data.
We also calculate the mean-field and projected PESs with DD-PC1 [78] and PC-F1 [56] parameter sets. The results are qualitatively the same as those we present here. Microscopically, the emergence of octupole collectivity stems from the specific shell structures of singleparticle levels around the Fermi level [1,79]. In Fig. 4 we plot the mean-field single-neutron levels of 96 Zr as functions of β 20 (β 30 ) deformation. The single-proton levels below and around Z = 40 are qualitatively similar and not shown here. For the spherical shape β 20,30 = 0, the neutron Fermi level lies between the 2d 5/2 and 1g 7/2 shells. In the left half of Fig. 4, we present the dependence of the single-neutron levels on pure β 30 deformation. We see that the energy splittings caused by non-zero β 30 are small, and the spherical gaps at N = 40, 50, and 70 are still remarkable. In the right half of Fig. 4, we show the splitting of the neutron energy levels due to a pure β 20 deformation. Here the 1h 11/2 state lying above the Fermi level splits into six levels. Thus, neutrons have the possibilities to occupy the 1h 11/2 orbital with considerable amplitude, leading to strong octupole coupling with the orbital 2d 5/2 . As for single-proton energy levels, the 2p 3/2 and 1 f 5/2 are the two states coupled via pure β 20 deformation while the 1g 9/2 and 2p 3/2 can be coupled via pure β 30 correlation. The energy splitting for non-zero β 30 is again small compared with that for non-zero β 20 , resulting in a large shell gap around the proton Fermi surface. For non-zero β 20 , some proton levels from the 1g 9/2 and 2p 3/2 states are close to each other near the Fermi surface. Thus, the octupole deformation in 96 Zr is dominated by the proton levels 2p 3/2 → 1g 9/2 and neutron levels 2d 5/2 → 1h 11 Fig. 2(b), respectively. The principal, orbital, and angular momentum quantum numbers nl j for spherical symmetry are presented. The dashed line denotes the Fermi level.

Summary
The competitions of various shell structures near A 100 create complex structures on the potential energy surfaces (PESs) of Zr isotopes. In particular, it was found that both triaxial and octupole deformations play important roles in 96 Zr, which presents a great challenge for calculations based on microscopic theories. Recently we implement a unified framework called p-MDCRHB model to calculate the potential energy surfaces with various symmetries. In this model both axial and reflection symmetries can be broken in the mean-field level and restored by projection-aftervariation methods. Various symmetries can be imposed to speed up the calculations. With this model, we show that the PESs of 96 Zr have strong dependence on the angular momentum and parity. Most mean-field and projected PESs are soft for triaxial deformation β 22 and tetrahedral deformation β 32 . We also extract the energy levels at the global minima on each projected PES and compare them with the experiments. Our results suggest that both triaxial and octupole degrees of freedom are essential elements for precisely describing the structure of the 96 Zr nucleus. A more solid conclusion cannot be made before the full symmetryrestored GCM calculation with the mixing of all types of deformed configurations is carried out. Work along this direction is in progress.