Nuclear masses near $N=z$ from Nilsson-Strutinsky calculations with pairing corrections beyond BCS from an isospin-conserving pairing force

Since the late 1990s the masses of nuclei near the line $N=Z$ in the chart of nuclides, where $N$ and $Z$ are the numbers of neutrons and protons, have attracted much interest from the nuclear physics community. In particular, the origin of the so-called Wigner energy Myers and Swiatecki (1966), a depression of the mass at $N=Z$ relative to a trend described by a symmetry energy quadratic in $N-Z$, has been a matter of debate. For a review of this discussion, see Neergård Neergård (2009). In 1995, Duflo and Zuker Duflo and Zuker (1995) published a semiempirical mass formula with a symmetry energy proportional to $T(T+1)$, where $T$ is the isospin quantum number, in the ground states of doubly even nuclei equal to $|T_z|$, where $T_z=(N-Z)/2$. This expression includes a Wigner energy in a natural manner. It was observed by Frauendorf and Sheikh Frauendorf and Sheikh () that a symmetry energy with the factor $T(T+1)$ resembles the spectrum of a quantal, axially symmetric rotor. These authors identified the nuclear superfluidity as the deformation in isospace that could give rise to collective rotation in this space. The Bardeen-Cooper-Schrieffer (BCS) pair gaps ${\Delta }_{n,p}$ are indeed components of an isovector perpendicular to the isospin. Their magnitude is a measure of the collectivity of the isorotation. Obviously the dependence of the energy on $T$ will become more regular as their magnitude increases as a result of the progressively more gradual change of the occupation numbers around the Fermi surface, which will wash out fluctuations of the level density. On the other hand the isorotational moment of inertia, which is determined by the average level density, will not change much. The results of calculations in Ref. Satuła and Wyss (2001) exemplify these generic features.

Neergård Neergård (2002, 2003, 2009) set up a microscopic theory of such a superfluid isorotation based on the Hartree-Bogolyubov (HB) plus random phase approximation (RPA). His Hamiltonian involves independent nucleons in a charge-independent potential well, an isovector pairing force, and an interaction of the nucleonic isospins, which he calls the symmetry force. The latter is shown to contribute merely a term in the total energy proportional to $T(T+1)$. In the idealized case of equidistant single-nucleon levels, the total symmetry energy is found to be proportional to $T(T+1)$ provided the pairing force is sufficiently strong to produce an HB energy minimum with nonzero BCS gaps. If this condition fails to be satisfied, the model gives in this case a total symmetry energy proportional to $T(T+X)$ with $X<1$ Neergård (2003). In realistic cases with nonuniform single-nucleon spectra, major modifications of these simple expressions arise from shell effects Neergård (2009).

Recently, Bentley and Frauendorf Bentley and Frauendorf (2013) calculated exactly the lowest eigenvalue of Neergård’s Hamiltonian employing small valence spaces of dimensions 24 or 28. They demonstrated (see their Fig. 7) that for a sufficiently strong pairing interaction the $T(T+1)$ limit of rigid isorotation is approached for various kinds of bunched single-particle level distributions. Unlike Neergård, who keeps the Hamiltonian constant along each isobaric chain, these authors take into account the variation of the nuclear shape with the isospin. In a survey of the range of mass numbers $A=N+Z=$ 24–100, they find that their model accounts generally for the fluctuations due to the shell structure observed in several combinations of the masses near $N=Z$ taken as functions of $A$: (i) the mass difference of doubly odd and doubly even $N=Z$ nuclei, (ii) the difference in excitation energy of the lowest $T=1$ and $T=0$ states in the the doubly odd $N=Z$ nuclei, and (iii) the constants $\theta$ and $X$ in an expansion of the symmetry energy of the form $T(T+X)/2\theta$ extracted from doubly even masses close to $T=0$. However, the model underestimates $\theta$ when the symmetry force constant is fit to the difference of $T=1$ and $T=0$ excitation energies.

The exact results of Bentley and Frauendorf provide a background on which the accuracy of the HB + RPA may be tested. We show in the present study that the HB + RPA gives a very good approximation to the exact lowest eigenvalue of the Hamiltonian except near the values of the pairing force coupling constant $G$ where the HB solution shifts from a zero to a non-zero pair gap. We show, as well, that the HB + RPA reproduces the exact eigenvalue asymptotically in the limit $G\to \infty$. We then devise a recipe for interpolating the HB + RPA energy across the critical region of $G$ leading to an algorithm which accurately approximates the exact eigenvalue in the entire range of $G$ from zero to infinity.

This algorithm is simple enough to allow calculations with valence spaces of dimension $2A$. More specifically, we include in our present calculations all single-nucleon states below the $T=0$ Fermi level and equally many states above this level. To allow actual nuclear masses to be calculated from this schematic model, we add a Strutinsky renormalization. As a side effect, we can then dispense with the symmetry force. Its contribution to the total energy proportional to $T(T+1)$ may thus be considered a part of the macroscopic liquid drop symmetry energy. Our microscopic Hamiltonian thus consists merely of a charge-independent independent-nucleon term and the isovector pairing force.

Using this scheme we calculate once more the combinations of masses near $N=Z$ previously considered by Bentley and Frauendorf and show that the issue with the underestimation of the constant $\theta$ is resolved by the enlargement of the valence space. We also show that the present scheme allows to calculate the masses themselves with a small root-mean-square deviation from the measured ones.

The following Secs. II–IV and VI contain the formal presentation of our model. In Sec. V we compare the exact lowest eigenstate of our microscopic Hamiltonian with the HB + RPA and explain our scheme of interpolation. Then, in Sec. VII, we present and analyze our results of calculation. Sec. VIII provides arguments detailing why we disregard isoscalar pair correlations. Finally, in Sec. IX, the article is summarized and some perspectives drawn.

The five parameters in Eq. (2) are determined by fitting Eq. (1) to the 112 measured binding energies of doubly even nuclei considered in the present study according to the 2012 Atomic Mass Evaluation Audi et al. (2012) with $E_{\text{s.p.}}$, ${\tilde{E}}_{\text{s.p.}}$, $P$, and $\tilde{P}$ calculated as described in Secs. III–VI. The result is $a_v=15.19$ MeV, $a_{vt}=110.8$ MeV, $a_s=16.35$ MeV, $a_{st}=135.3$ MeV, and $a_c=0.6615$ MeV. These parameters are similar to those obtained by Mendoza-Temis et al. in a global fit of nuclear masses with minus the binding energy given by Eq. (2) plus a phenomenological, negative definite, pairing energy and $B_S=B_C=1$ Mendoza-Temis et al. (2010). The parameters of Mendoza-Temis et al. cannot be used in the present context because unlike their negative definite pairing energy the sum of the liquid drop deformation energy and the shell and pairing correction terms in Eq. (1) average approximately to zero. Moreover, the semiempirical formula of Mendoza-Temis et al. deviates quite a lot from the empirical binding energies locally in the present region of nuclei. This deviation increases with increasing $A$ along the $N=Z$ line and amounts to about 20 MeV for $A\approx 100$. Since our liquid drop parameters are optimized for the present region of $N$ and $Z$ they may not reproduce accurately binding energies in other regions.

The sum $E=E_{\text{s.p.}}+P$ in Eq. (1) is calculated as the lowest eigenvalue of the Hamiltonian

where $k$ labels orthogonal quartets of a single-proton and a single-neutron state and their time reversed. The annihilator of a nucleon in one of these states is denoted by ${\hat{p}}_k$, ${\hat{n}}_k$, ${\hat{p}}_{\overline{k}}$, or ${\hat{n}}_{\overline{k}}$, and

The single-nucleon energies ${ϵ}_k$ are derived from a calculation with the Nilsson potential employing the parameters of Bengtsson and Ragnarsson Bengtsson and Ragnarsson (1985). To conserve isospin we take the average of the neutron and proton energies with given ordinal number counted from the bottom of the spectrum. The resulting quartets are labeled by $k$ in the order of increasing ${ϵ}_k$, and the first $\Omega =A/2$ of them in this order included in the calculation.

The deformation parameters ${ϵ}_2$ and ${ϵ}_4$ are taken from a recent survey of deformations based on the Nilsson-Strutinsky plus BCS theory. The equilibrium deformation of a given nucleus was calculated in this survey by minimizing with respect to Larsson’s triaxial deformation parameters $ϵ$ and $\gamma$ Larsson (1973) as well as ${ϵ}_4$ the Nilsson-Strutinsky plus BCS energy calculated with the tac code Frauendorf (1993). The latter employs an expression for the liquid drop energy similar to Eq. (2) but with symmetry energy terms quadratic in $T_z$, the Nilsson potential with the parameters of Ref. Bengtsson and Ragnarsson (1985), and a pairing correction without a smooth counterterm calculated from BCS gaps ${\Delta }_n={\Delta }_p=12A^{-1/2}$ MeV. In calculations for odd $N$ or $Z$ the Fermi level is blocked. All the 136 nuclei included in our study turn out to have either $\gamma =0$ or $\gamma ={60}^{\circ }$, that is, axial symmetry. The deformations are shown in Tables 2 and 3. Also see Note sup ().

The only free parameter remaining is the pairing force coupling constant $G$. A power law for $G$ as a function of $A$ will be fit to $T=0$ doubly even-doubly odd binding energy differences.

For an introduction to the BCS, HB, and RPA theories we refer to textbooks such as the one by Ring and Schuck Ring and Schuck (1980). The calculation of the lowest eigenvalue of the Hamiltonian (3) in the HB + RPA for even $N$ and $Z$ is discussed by Neergård Neergård (2009). Since this formalism is invariant under isorotation, a nucleus with $N\ge Z$ represents the entire multiplet with $T=T_z=(N-Z)/2$. To calculate the energy of the lowest $T=0$ state for odd $N=Z$ we reduce this case to the even one by omitting quartet number $k=(N+1)/2$ from a HB + RPA calculation for $N-1$, $Z-1$ and adding $2{ϵ}_k$.

The HB part of the calculation amounts to the usual BCS theory with quasinucleons annihilators

For a gap $\Delta =0$ the chemical potential $\lambda$ is taken as the limit for $\Delta \to 0$ of the $\lambda$ determined by $N$ or $Z$ for $\Delta >0$.

The RPA part splits into separate equations for an “$n$ space” spanned by ${\hat{\alpha }}_{\overline{k}}{\hat{\alpha }}_k$ and their Hermitian conjugates, a “$p$ space” spanned by ${\hat{\beta }}_{\overline{k}}{\hat{\beta }}_k$ and their Hermitian conjugates, and an “$np$ space” spanned by ${\hat{\beta }}_{\overline{k}}{\hat{\alpha }}_k+{\hat{\alpha }}_{\overline{k}}{\hat{\beta }}_k$ and their Hermitian conjugates. The resulting groundstate energy can be written as $E_{\text{s.p.}}+P$ with $P=P_{\text{BCS}}+P_{\text{RPA}}$ and

Here ${\Delta }_n$ and ${\Delta }_p$ are the neutron and proton BCS gaps and

where $\omega$ denotes an RPA frequency and $E_{k,n\text{or}p}$ are the single-quasinucleon energies. (A term $c$ given by Eq. (35) of Ref. Neergård (2009) vanishes in the present case of $A=2\Omega$.)

For nuclei with $N=Z$, the BCS solutions are equal for neutrons and protons and the $n$, $p$, and $np$ spaces have equal RPA spectra so that $E_{\text{RPA,}n}=E_{\text{RPA,}p}=E_{\text{RPA,}np}$. If ${\Delta }_n>0$ the lowest RPA frequency in the $n$ space vanishes because the Hamiltonian (3) commutes with $N$, and the analogon of this statement holds for protons. For $N>Z$ the lowest frequency in the $np$ space is equal to the difference ${\lambda }_n-{\lambda }_p$ of the neutron and proton chemical potentials because the Hamiltonian commutes with the components of the isospin perpendicular to the $z$ direction.

Bentley and Frauendorf calculated exactly the lowest eigenvalue of the Hamiltonian (3) in spaces with six or seven quartets Bentley and Frauendorf (2013). This allows a comparison of the HB + RPA to an exact calculation. We made this comparison in all the cases displayed in Fig. 7 of Ref. Bentley and Frauendorf (2013). The case shown in Fig. 1

is the one with the largest deviation of the HB + RPA result from the exact energy. The HB + RPA curve is seen to follow closely the exact one except in a small interval about $G=ϵ/6$, which is in this case the critical value $G_{\text{crit.}}$ of $G$ where the BCS solution (equal in this case for neutrons and protons) changes from $\Delta =0$ to $\Delta >0$. While the exact groundstate energy has a smooth variation across $G_{\text{crit.}}$, the HB + RPA curve shows there a prominent cusp.

The origin of this cusp can be traced to the expression (7) for the RPA contributions. Thus notice the plot in Fig. 2

of the two lowest RPA frequencies in the case just considered in any of the $n$, $p$, and $np$ spaces, which have in this case identical RPA spectra because $N=Z$. Both frequencies are seen to go to zero for $G$ going to $G_{\text{crit.}}$ from below, but only one of them stays at zero for $G>G_{\text{crit.}}$ while the other one rises rapidly in this interval. The reason why two frequencies and not only one go to zero for $G$ going to $G_{\text{crit.}}$ from below is that the quasinucleon vacuum becomes in this limit instable against a transition to a vacuum described by a non-zero $\Delta$ with an arbitrary complex phase. This $\Delta$ has two real parameters. In more physical terms, Fig. 1 may be interpreted to display a shortcoming of the RPA, which is a small amplitude approximation, in a region of the parameter $G$ where the equilibrium represented by the quasinucleon vacuum changes rapidly with this parameter.

An example with less symmetry of the single-nucleon spectrum is given in Fig. 3.

It is seen that while the lowest RPA frequencies in the $n$ and $p$ spaces behave as in the preceding case, the lowest RPA frequencies in the $np$ space pass almost smoothly through the critical $G$, which are different in this case for neutrons and protons. This is general for $N>Z$.

These observations suggest that the HB + RPA may be improved by interpolation across the critical regions of $G$. More specifically we have found that one gets a good approximation to the exact groundstate energies by interpolation in the interval from $0.5G_{\text{crit.}}$ to $1.5G_{\text{crit.}}$. This is applied to the terms $E_{\text{RPA,}n}$, $E_{\text{RPA,}p}$, and $E_{\text{RPA,}np}$ for $N=Z$, when neutrons and protons have the same $G_{\text{crit.}}$, and to $E_{\text{RPA,}n}$ and $E_{\text{RPA,}p}$ for $N>Z$, when $G_{\text{crit.}}$ may be different for neutrons and protons. No interpolation is applied if $G_{\text{crit.}}=0$, which occurs when the Fermi level lies within a degenerate shell. The interpolating function is the polynomial of third degree in $G$ which joins smoothly the calculated values at the endpoints of the interpolation interval.

Figure 4

shows the result of using this recipe in the “worst case” of Fig. 1, and Fig. 5

shows its effect for the doubly magic nucleus ${}^{56}$Ni, which resembles the case of Fig. 1 by having its Fermi level (common for neutrons and protons) within a gap in its single-nucleon spectrum. In the case of ${}^{56}$Ni we have no exact calculation for comparison, but the interpolation is seen to remove a certainly unphysical cusp from the curve of the groundstate energy as a function of $G$.

We finish this section with a discussion of the limit $G\to \infty$. This discussion is restricted to the case of even $N$ and $Z$. It is not restricted to $A=2\Omega$. Without loss of generality the centroid of the single-nucleon spectrum is supposed to vanish. First consider the case of degenerate single-nucleon levels, that is, ${ϵ}_k=0$ for all $k$ so that the Hamiltonian (3) has only the second term, the pairing force. Our $4\Omega$-dimensional valence space is then equivalent to a $j$-shell with $j=\Omega -1/2$. From the formulas in Ref. Neergård (2009) one gets in this case

Exactly this expression for the lowest eigenvalue of the pairing force results from the formulas derived by Edmonds and Flowers Edmonds and Flowers (1952) by means of group theory. In other words, for degenerate single-nucleon levels the HB + RPA gives the exact result. Now assume a spreading of the levels ${ϵ}_k$. Due to their vanishing centroid the first term in the expression (3) is then, in the $j$-shell analogy, a sum of spherical tensor components of rank higher than zero. Since its expectation value in the angular momentum zero ground state of the degenerate case then vanishes, its contribution to the energy vanishes in the Born approximation, and the leading term in this contribution in an expansion in powers of $G^{-1}$ is the linear term. In other words, for a general single-nucleon spectrum the HB + RPA result converges asymptotically to the exact one in the limit $G\to \infty$.

The smooth sum ${\tilde{E}}_{\text{s.p.}}$ of single-nucleon energies in Eq. (1) is calculated separately for neutrons and protons by a standard third order Strutinsky smoothing Brack et al. (1972) with smoothing parameter $\Gamma =41A^{-1/3}$ MeV including energies until approximately $\tilde{\lambda }+5\Gamma$, where $\tilde{\lambda }$ is the smooth Fermi level. Corresponding to the spitting of the pairing energy $P$ into a BCS and RPA part discussed in Sec. IV we write

Here the part ${\tilde{P}}_{\text{BCS}}$ is calculated separately for neutrons and protons essentially as suggested by Brack et al. Brack et al. (1972): We relate the paring strength $G$ to a smooth pair gap $\tilde{\Delta }$ by considering a half-filled single-nucleon spectrum of $\Omega$ Kramers doublets with equal distance $1/\tilde{g}\left(\tilde{\lambda }\right)$, where $g\left(\tilde{\lambda }\right)$ is the smooth level density at the smooth Fermi level, and by replacing the sum in the gap equation with an integral. When this integral is evaluated more accurately than in Ref. Brack et al. (1972) one arrives at

The BCS correlation energy $P_{\text{BCS}}$ can be expressed by a sum over the Kramers doublets, which can be approximated in a similar way by an integral. By evaluating also this integral more accurately than in Ref. Brack et al. (1972), one gets

A derivation of an expression for a smooth RPA correlation energy can be based on Eq. (38) of Ref. Neergård (2009). We treat first the case of an $n$ or $p$ space and consider again a half-filled single-nucleon spectrum of $\Omega$ Kramers doublets with equal distance $1/\tilde{g}\left(\tilde{\lambda }\right)$. In Eq. (38) of Ref. Neergård (2009) the sum over Kramers doublets can be replaced with an integral in the propagators $G_0({\hat{P}}_{\tau },{\hat{P}}_{\tau },\omega )$ etc. with ${\hat{P}}_{\tau }={\sum }_k{\hat{P}}_{\tau ,k}$. Some mathematics then leads to

In the derivation of this expression the integral in Eq. (38) of Ref. Neergård (2009) is displaced to the imaginary $\omega$ axis. This is allowed because $G_0({\hat{P}}_{\tau },{\hat{P}}_{\tau },\omega )$ etc. are asymptotically proportional to ${\omega }^{-2}$.

As in Ref. Bentley and Frauendorf (2013), we fit an expression for the pairing force coupling constant $G$ proportional to a power of $A$ to the empirical $T=0$ binding energy differences $2\Delta$ of the even-even and odd-odd $N=Z$ nuclides defined by

for odd $A/2$, cf. the definition of $B(A,T,T_z)$ in the beginning of Sec. II. (Note that this $\Delta$ is different from the BCS pair gap $\Delta$ considered in Secs. IV and V.) We take the groundstate binding energies from the 2012 Atomic Mass Evaluation and the $T=0$ excitation energies for odd $A/2$ from the NNDC Evaluated Nuclear Structure Data Files Tu1 (). The best fit, shown in Fig. 6,

is obtained for

The calculation is seen to reproduce the observed pattern of fluctuations due to shell structure of $2\Delta$ as a function of $A$ about a trend line $2\Delta =24A^{-1/2}$ MeV. For a given $A$, the calculated and observed values thus lie consistently both above or both below this line.

Our further analysis involves Coulomb reduced binding energies. We assume that the electromagnetic contribution to the total energy is given by the last term in Eq. (2) and thus write for the remainder

This reduction is applied to both the calculated and the measured binding energies. The error bars shown in Figs. 6–9 include the uncertainties of the empirical mass differences involved and the uncertainty of $a_c$ in the fit of Eq. (2) to the observed masses. When no error bar is shown, the uncertainty is less than the size of the symbol.

To the extend that the last term in Eq. (16) may be assumed to account for all contributions to the total energy from non-isospin-conserving interactions, $E_S(A,T)$ is independent of $T_z$. On the right hand side of Eq. (16), we mostly choose $T_z=T$. The only exception is that we compare the difference $E_S(A,1)-E_S(A,0)$ for odd $A/2$ measured for $T_z=0$ with the one calculated for $T_z=T$.

Figure 7

shows the measured and calculated values of this difference. The measured $E_S(A,1)-E_S(A,0)$ is the difference in excitation energy of the lowest $T=1$ and $T=0$ states of the $N=Z=A/2$ nucleus and is taken from its Evaluated Nuclear Structure Data File. The observed trend of a shift from $T=0$ ground states of the lighter doubly odd $N=Z$ nuclei in the range $A=$ 26–98 to $T=1$ ground states of the heavier ones as well as the average slope of the excitation energy difference as a function of $A$ are well reproduced. As discussed by Vogel Vogel (2000) and Macchiavelli et al. Macchiavelli et al. (2000), $E_S(A,1)-E_S(A,0)$ can be interpreted as a difference between the $T=1$ symmetry energy and the cost in energy $2\Delta$ of breaking a Cooper pair. Its downslope as a function of $A$ then results from an increase of the latter relative to the former. It may be noticed that $E_S(A,1)-E_S(A,0)$ is the only one of the four combinations of energies displayed in Figs. 6–9 where measured and calculated energies with different $T_z$ are compared and the difference of the energies compared is thus influenced by the Coulomb reduction (16).

We finally consider the quantities $\theta$ and $X$ defined by

where the constant is a function of $A$, and $A/2+T$ is even. The constant $\theta$ can be interpreted as an isorotational moment of inertia Frauendorf and Sheikh (), and $X$ is related to the Wigner energy and may be called the Wigner $X$. In this analysis, $T_z=T$ is chosen throughout on the right hand side of Eq. (16), so only doubly even nuclei, and therefore only groundstate energies, are involved. The empirical groundstate binding energies are taken from the 2012 Atomic Mass Evaluation. The constants $\theta$ and $X$ are extracted from the calculated and measured $E_S(A,T)$ for $T=0,2,4$ in the case of even $A/2$ and $T=1,3,5$ in the case of odd $A/2$.

Figures 8

and 9

show the result of this analysis. The measured $\theta$ is very well reproduced including, generally speaking, the features seen in the experimental values. This is a vast improvement from the calculations in Ref. Bentley and Frauendorf (2013), which considerably underestimated $\theta$. This confirms what was suggested in Ref. Bentley and Frauendorf (2013), namely that the underestimate there originated in the small size of the valence space. For $60<A<80$ the calculation is very accurate, indicating that the deformations used in the calculations are accurate.

As already mentioned in the introduction, several effects contribute according to Neergård Neergård (2009) to the linear term $XT/2\theta$ in the expansion (17). One such contribution, corresponding to $X=1$, comes from the expression $T(T+1)/2\theta$ for an isorotational energy. This is an average contribution, represented in our present model by the linear term in the liquid drop symmetry energy in Eq. (2). Microscopically, a symmetry energy proportional to $T(T+1)$ is shown in Ref. Neergård (2009) to emerge from the HB + RPA in the idealized case of equidistant single-nucleon levels. It is dependent in this case on the spontaneous breaking of the isobaric invariance by the isovector pairing force. In fact, in the case of equidistant single-nucleon levels the symmetry energy is proportional to $T(T+X)$ with $X<1$ in the absence of static pair fields Neergård (2003). The numerical solutions for this case in Fig. 7 of Ref.Bentley and Frauendorf (2013) show that $X$ approaches one from below with increasing pair coupling strength $G$. Equidistant levels are a particular favorable case. The same figure demonstrates that $X=1$ is approached for uneven level distributions as well, because a large $G$ represents the limit of strong deformation in isospace, which results in rigid isorotation.

On top of this isorotational contribution major contributions to the linear term in the symmetry energy arise from shell effects. This is discussed by Neergård Neergård (2009) and further elaborated by Bentley and Frauendorf Bentley and Frauendorf (2013). They relate the deviations of $X$ from one to the deviation of the distance between the last occupied and first free level in the absence of pair correlation (cf. Fig. 8 of Ref. Bentley and Frauendorf (2013)). The realistic pair correlation is too weak to wash out this consequence of the level bunching. Neergård shows, in particular, that $X$ is large when the $T=0$ Fermi level lies within a gap in the single-nucleon spectrum. The reason is that the isospin is then produced by promotion of pairs from proton levels below the gap to neutron levels above the gap, and each such promotion costs approximately the same amount of energy equal to twice the gap energy. That the empirical $X$ is mostly larger than one may thus be seen as the result of the level density at the $T=0$ Fermi level being generally at equilibrium shape lower than corresponding to a uniform spectrum. By analogy with spatial rotation of nuclei it may be seen as evidence for a “softness” of isorotation. The ratio $R_{42}=E_{4^{+}}/E_{2^{+}}$, which is 10/3 for a rigid rotor and 2 for a harmonic vibrator, is commonly used as a measure of how “rotational” a nucleus is. The energy of the first few yrast levels of even-even transitional nuclei can be very well parametrized by the expression $E\left(I\right)=I(I+X)/2J$. This gives $X-1=(10-3R_{42})/(R_{42}-2)$, so $X>1$ is equivalent to $2<R_{42}<10/3$.

As seen from Fig. 9, our model reproduces very well these fluctuations of $X$ due to shell effects. For $A=$ 90, 92, 98, and 100, the calculated valued are markedly below the empirical ones. It should be noted that the corresponding values of $1/\theta$ in Fig. 8 are markedly above the empirical ones. If one would take the product of both numbers to obtain $X/\theta$, which is twice the coefficient of the term linear in $T$ in the expansion (17), a much better agreement would result.

As anticipated from the discussion above, large $X$ occur in the isobaric chains containing the doubly magic nuclei ${}^{40}$Ca, ${}^{56}$Ni, and ${}^{100}$Sn. Both chains with odd $A/2$ neighboring each of these chains with even $A/2$ also have large $X$. This is because in the odd-$A/2$ chains one pair of neutrons is passive in producing isospin. This pair just sits for all $T=1,3,5$ in the last quartet below the gap or the first quartet above the gap while other nucleons are promoted across both the gap and this quartet repeating the mechanism described above. A somewhat similar mechanism gives rise to a large $X$ for $A=30$. There one neutron pair occupies for all $T=1,3,5$ the $2s_{1/2}$ shell while other nucleons are promoted from the $1d_{5/2}$ to the $1d_{3/2}$ shell.

For $A=24$ and $A=48$, deformation is involved in producing a large $X$. In these chains, the $T=0$ nuclei thus have large deformations while the $T=2$ and $T=4$ nuclei are essentially spherical. Since isospin is produced by promotion of nucleons within the $1d_{5/2}$ or $1f_{7/2}$ shell (except that in ${}^{24}$O one neutron pair occupies the $2s_{1/2}$ shell), the two spherical nuclei have roughly equal $E_{\text{s.p.}}$. We therefore have a large increase of $E_{\text{s.p.}}$ from $T=0$ to $T=2$ due to the departure from the deformed shape and essentially no such increase from $T=2$ to $T=4$. This gives rise to a large $X$. It also gives a small $1/\theta$ as seen in Fig. 8.

Tables 2 and 3 show the calculated and measured binding energies of the individual nuclei and components of the calculated ones. Also see Note sup (). The root mean square deviation is 0.95 MeV for the 112 doubly even nuclei with a measured binding energy.

It has been proposed Satuła and Wyss (); Goodman (); B Cederwall et. al. (2011); Qi et al. (2011) that the strong attraction of isoscalar nucleon pairs exhibited by effective shell model interactions, especially in channels with maximally aligned nucleonic angular momenta, could give rise, in $N=Z$ nuclei, to a condensation of such pairs coexisting with or replacing the BCS type of condensation of isovector pairs. We feel that there are points of contention with this proposal: (i) The energy of a condensate depends smoothly on its number of constituents. In the presence of a condensate of isoscalar pairs the mass of the lowest $T=0$ state of an $N=Z$ nucleus should therefore depend smoothly on $A$. In reality, these masses show a staggering with the parity of $A/2$, the doubly odd masses being elevated above the doubly even ones by an amount approximately twice the typical BCS pair gap; see Fig. 6. (ii) Bentley and Frauendorf, in their aforesaid study, examine the effect of adding to the Hamiltonian a schematic interaction of isoscalar pairs of a neutron and a proton in time-reversed orbits with an separable structure similar to that of the isovector pairing force. They find that a weak interaction of this form does not significantly alter their results, while a stronger one would not allow the model to reproduce the data. (iii) In the single-$j$-shell model, seniority zero represents condensation of isovector pairs in the BCS sense. Neergård analyzed the ground states of nuclei with two neutrons and two protons or two neutron holes and two proton holes in the $1f_{7/2}$ or $1g_{9/2}$ shell calculated in the single-$j$-shell approximation with effective interaction from the literature Neergård (2013). He found these states to have by about 80% seniority zero. As pointed out in Neergård’s study, since the seniority zero state has a considerable contingent of isoscalar pairs, the attractive interaction of such pairs stabilizes the seniority zero component of the state vector rather than competing with it.

For these reasons, we do not consider the possibility of condensation of isoscalar pairs in our present work. Some studies, for example Refs. Martínez-Pinedo et al. (1999); Qi et al. (2011), infer a pairing structure of a shell model state from counts of nucleon pairs with given angular momentum. Since Neergård demonstrates in Ref. Neergård (2013) that such a count is not a reliable tool for this purpose, we have not taken such work into account in the discussion in the preceding paragraph.

In Strutinsky-type calculations including schematic isoscalar and isovector pair correlations, Głowacz et al. Głowacz et al. (2004) achieved results for $E(T=1)-E(T=0)$ in doubly odd $N=Z$ nuclei in agreement with the data similar to that of our results displayed in Fig. 7. This is consistent with the findings of Bentley and Frauendorf mentioned in the first paragraph of this section. Although such a mixed scenario cannot be excluded we follow the principle of Occam’s razor, assuming pure isovector pair correlations.

For $A=48$ we made a supplemental $1f_{7/2}$ shell model calculation using the interaction “model I” of Zamick and Robinson Zamick and Robinson () with a normalization to zero of the largest matrix element (two-nucleon angular momentum $J=6$) so as to make the total interaction attractive. Like Satuła et al. Satuła et al. (1997), we switched off successively the interactions in individual $J$ channels. The results for $\theta$ and $X$ are shown in Table 1.

Note that because the $J=6$ matrix element is normalized to zero, it makes no difference whether it is included or not. Like $\theta$ and $X$ the quantity ${ϵ}_W$ displayed in Fig. 2 of Ref. Satuła et al. (1997) is a function of $E_S(A,T)$ for $T=0,2,4$. The relation is ${ϵ}_W=2X/3\theta$. Like Satuła et al., who include in their calculations the shells $2p_{3/2}$, $1f_{5/2}$, and $2p_{1/2}$ and employ an interaction appropriate for this larger valence space, we find that $X/\theta$ and therefore ${ϵ}_W$ decreases when the isoscalar interactions $J=1,3,5,7$ are switched off successively and very much so when the $J=7$ interaction is switched off finally. It is seen, however, that this is due not to a decrease of $X$, which actually increases, but to a decrease of the symmetry energy coefficient $1/2\theta$. The isoscalar shell model interactions thus contribute significantly to the entire symmetry energy and not just its Wigner term. The reduction of the Wigner energy when the isoscalar interactions are switched off is only a side effect of this general reduction of the symmetry energy. The symmetry energy coefficient reaches its minimum when all the isovector interactions $J=0,2,4,6$ and none of the isoscalar interactions $J=1,3,5,7$ are present, and the calculations confirm the well known result derived analytically by Edmonds and Flowers Edmonds and Flowers (1952) that $X$ is exactly one for the pure pairing force, $J=0$. In this case $1/\theta =G$ as seen from Eq. (8).

Our approach makes further simplifying assumptions: (i) The Coulomb interaction can be treated as a first-order perturbation, that is, its contribution to the total energy may be approximated by its expectation value with the wave function determined by the strong interaction only. This contribution can then be incorporated in the form of the Coulomb term in the smooth liquid drop binding energy formula (or subtracted from the experimental binding energies, as done in this paper). (ii) Isospin-breaking terms of the strong interaction Hamiltonian (as the difference between the proton and neutron masses) are neglected. (iii) The difference between proton and neutron mean fields generates only a constant shift of the proton single-particle levels relative to the neutron ones. This constant shift drops out in the shell correction procedure so that one can assume the same single proton and neutron energies from the outset. These assumptions lead to our isospin-invariant Hamiltonian (3) used to calculate shell and pairing corrections.

Sato et al. Sato et al. (2013) studied $N=Z$ nuclei in the framework of the density functional mean field theory. They describe the $T>0$ states by isocranking about an axis in isospace that is tilted with respect to the $z$ axis. The resulting quasiparticles are mixtures of proton and neutron particles and holes. As discussed by Frauendorf and Sheikh Frauendorf and Sheikh (), for an isospin invariant Hamiltonian all mean field solutions that correspond to the same cranking frequency but a different orientation of the cranking axis (the “semicircle” of Ref. Sato et al. (2013)) have the same energy. They can be generated by rotation in isospace from the solution obtained by cranking about the $z$ axis, which has pure proton and pure neutron quasiparticles. The rotation generates a mixing of the proton and neutron quasiparticles (cf. the example of isocranking about the $y$ axis discussed in Ref. Frauendorf and Sheikh ()). Hence if our assumptions hold, it is sufficient to study isosrotation about the $z$ axis, which generates $T_z=T$ solutions and avoids proton-neutron mixing. The $T_z<T$ solutions are given by rotation in isospace. If the mean field theory includes the Coulomb interaction, as the study by Sato et al. Sato et al. (2013) does, the different orientations of the cranking axis are no longer equivalent, and the orientation of minimal energy has to be calculated. However the finding of Sato et al. that such solutions lie with a good accuracy on a shifted semicircle indicates that our assumptions are good approximations.

A model with nucleons in a charge independent potential well interacting by an isovector pairing force has been discussed. For a 24-dimensional valence space, the Hartree-Bogolyubov (HB) plus random phase approximation (RPA) to the lowest eigenvalue of the Hamiltonian was shown to be accurate except near the values $G_{\text{crit.}}$ of the pairing force coupling constant $G$ where the HB solution shifts from a zero to a non-zero pair gap. The HB + RPA was shown to be asymptotically exact in the limit $G\to \infty$. To remedy the inaccuracy of the HB + RPA in the critical regions of $G$ we devised a scheme of interpolation across the these regions. It is described in Sec. V.

The resulting algorithm was used to calculate with a valence space of dimension twice the mass number $A$ pairing corrections in the framework of a Nilsson-Strutinsky calculation. For this purpose we derived in Sec. VI expressions for smooth counterterms to the Bardeen-Cooper-Schrieffer (BCS) and RPA parts of the pair correlation energy. The deformations and corresponding single-nucleon energies for the Nilsson-Strutinsky calculation were taken from a previous Nilsson-Strutinsky plus BCS calculation with the code tac Frauendorf (1993). To enforce charge independence the average of the calculated single-neutron and single-proton energies was employed. Our expression (2) for the macroscopic liquid drop energy was taken from the work of Duflo and Zuker Duflo and Zuker (1995) with the omission of a phenomenological pairing energy and has symmetry energy terms proportional to $T(T+1)$, where $T$ is the isospin. Its five parameters were fit to the empirical masses according to the 2012 Atomic Mass Evaluation Audi et al. (2012) of the 112 doubly even nuclei with a measured binding energy considered in the present study.

In this model we calculated the binding energies of the ground states of the doubly even nuclei with $24\le A\le 100$ and $0\le N-Z\le 10$ and the lowest isospin $T=0$ states of the doubly odd nuclei with $26\le A\le 98$ and $N=Z$, where $N$ and $Z$ are the numbers of neutrons and protons. These calculated binding energies were compared to the empirical ones from the 2012 Atomic Mass Evaluation with $T=0$ excitation energies from the NNDC Evaluated Nuclear Structure Data Files Tu1 (). In terms of both the calculated and the empirical binding energies $B(A,T,T_z)$, where $T_z=(N-Z)/2$, a Coulomb reduced energy $E_S(A,T)$ was defined by Eq. (16) with $T_z=T$ on the right-hand side with one exception to be told later. The following combinations were then extracted and compared: (i) 2$\Delta$ as defined by Eq. (14). (ii) $E_S(A,1)-E_S(A,0)$ for odd $A/2$. In this case the measured $E_S(A,1)$ was defined with $T_z=0$ on the right hand side of Eq. (16). (iii) The constants $\theta$ and $X$ in the expansion (17) for all $A$.

Comparisons of the calculated and measured values of these combinations are shown in Figs. 6–9. The expression (15) adopted for the pairing force coupling constant $G$ was fit to the empirical 2$\Delta$. The present enlargement of the valence space resolved an issue with the constant $\theta$, which was underestimated in the previous exact calculation by Bentley and Frauendorf with a 28-dimensional valence space. The fluctuations of $X$ with $A$ were discussed in Sec. VII. They are well understood from the shell structure. The root mean square deviation of the calculated and measured binding energies of the 112 doubly even nuclei with a measures binding energy is 0.95 MeV.

We anticipate a generalization of the present method to the more realistic case when neutrons and protons move in different potential wells due to the Coulomb force. The chief obstacle to this generalization is the RPA calculation in the $np$ space, which will be more complex because neutron and proton stationary states no longer form time-reversed pairs. Simplifying approximations might be warranted in this step of the procedure, however. The result would be a method for including pairing correlations beyond a mean field approximation (BCS, Hartree-Fock-Bogolyubov, relativistic mean field) that would be simple enough to go on top of any state-of-the-art mean field approach. This would eliminate the need for the phenomenological Wigner term often employed in present mean field calculations such as, for example, those of Refs. Möller et al. (1995); Goriely et al. (2010); Wang and Liu (2013).

This work was supported by the DoE Grant DE-FG02-95ER4093.