Tensor force and deformation in even-even nuclei

The variational principle is used to build a model which describes open shell nuclei with ground state deformations. Hartree-Fock equations are solved by using single particle wave functions whose radial parts depend on the projection of the angular momentum on the quantization axis. Pairing effects are taken into account by solving Bardeen-Cooper-Schrieffer equations in each step of the minimisation procedure. The Gogny D1S finite-range interaction and an extension of it that includes tensor terms are consistently used in both parts of our calculations. The model is applied to study a set of isotopes with 34 protons and of isotones with 34 neutrons. Total energies, density distributions, their radii and single particle energies are analysed and the results of our calculations are compared with the available experimental data. We focused our attention on the effects of the deformation and of the tensor force on these observables. Our model describes open shell nuclei from a peculiar perspective and opens the possibility of future theoretical developments.


I. INTRODUCTION
The basic approach to describe nuclear deformations in terms of nucleons is the Nilsson model [1]. In this model, the nucleons move independently of each other in a deformed potential and the single particle (s.p.) wave functions depend on the projection of the corresponding angular momenta on the quantization z-axis. For a fixed value of the s.p. angular momentum, the states with smaller absolute values of this projection are more bound in case of prolate deformations, and the inverse happens in oblate nuclei [2,3].
While in the Nilsson model the deformed potential is an external input, in the model we present here this potential is obtained by considering an effective nucleon-nucleon interaction and the variational principle. By searching for the minimum of the energy functional in the Hilbert sub-space formed by Slater determinants, a set of Hartree-Fock (HF) equations is obtained [3] and we make the s.p. wave functions used to build up these Slater determinant explicitly dependent on the z-projection of their total angular momentum. In the minimisation procedure, we take care of the pairing interaction by carrying out Bardeen-Cooper-Schrieffer (BCS) calculations that modify the occupation probabilities of the s.p. states. In our model, which we call HFBCS, the deformation emerges because not all the s.p. states with the same total angular momentum, but with different z-axis projection, are occupied.
In this article, we present our model and apply it to the description of the ground state of even-even nuclei. As a testing ground, we have considered a set of medium-heavy nuclei where the occupation of the s.p. states ends in the f -p shell. Specifically, we have studied a set of Se isotopes, they have 34 protons, and a set of isotones with 34 neutrons. We focused our attention to the emergence of the deformations and, since we have considered effective nucleon-nucleon interactions which include tensor terms, on the role of these tensor terms .
The results obtained within the HFBCS model for the aforementioned even-even nuclei have been compared with those found in deformed Hartree-Fock-Bogoliubov (HFB) calculations and with the available empirical values of their binding energies and their charge radii and distributions. In addition, also the angular momenta of the neighboring odd-even nuclei have been analyzed. The relevance of the use of this new set of s.p. wave functions is studied by comparing our results with those of the spherical HF+BCS model of Ref. [4].
We present in Sec. II the theoretical background of our model, and, in Sec. III, the technical details of our calculations. In Sec. IV we show, and discuss, the results we have obtained, focusing on deformations and the role of the tensor force. The conclusions of our study are summarized in Sec. V.

II. THE MODEL
In our description of an even-even nucleus, composed by A nucleons, the basic ingredient is the set of s.p. wave functions used to build the Slater determinants. We assume that these s.p. wave functions, φ k (x), can be factorized as: where we have indicated with x the generalized coordinate, which includes the position r with respect to the nuclear centre, the spin and the isospin of the considered nucleon. The radial part of the s.p. wave function, is a function of r ≡ |r| and depends not only on the principal quantum number n k , on the orbital orbital angular momentum quantum number l k , on the total angular momentum quantum number j k , and on the isospin third component t k , but also on the projection of j k on the z-axis, m k . The part of the s.p. wave function depending on the angular coordinates, Ω k ≡ (θ k , φ k ), and on the spin third component, s k , is: where Y l k µ k is a spherical harmonic, the symbol | indicates a Clebsch-Gordan coefficient, and χ s k is a Pauli spinor. Finally, in Eq. (1), χ t k indicates the Pauli spinor related to the isospin. We assume time-reversal invariance [2]; this means that, in our approach, R t k n k l k j k , m k (r) = R t k n k l k j k ,−m k (r), and each s.p. level is two-fold degenerated. This also implies that the nucleus has the shape of an ellipsoid whose symmetry axis is the z-axis.
We build a Slater determinant Φ with the s.p. wave functions φ k of Eq. (1), and the application of the variational principle to search for the minimum of the corresponding energy functional, E[Φ], leads to a set of differential equations of the type: As indicated by the first term, that of the kinetic energy, this expression has been obtained by integrating, and summing, on the angular, and spin, coordinates. In the above equation, ǫ k indicates the s.p. energy, and the terms depending on the two-body effective nucleon-nucleon interaction, V (r 1 , r 2 ), are the so-called Hartree potential, the Fock-Dirac term, and a term related to the density dependence of the interaction: The effective force V (r 1 , r 2 ) used in our calculations is a finite-range interaction that includes the four traditional central terms (scalar, isospin, spin and spin-isospin), a zero-range spin-orbit term, a scalar density dependent term, and a tensor and a tensor-isospin dependent terms. More explicit expressions of U, W and K are presented in Appendix A.
In Eqs. (5)- (7), we have indicated with v 2 k the occupation probabilities of the s.p. states. Their values are obtained by solving the set of BCS equations: where λ, the energy gap, is given by: The quantity ∆ k satisfies the relation with |αα 00 indicating a state where the s.p. states φ α are coupled to total angular momentum J = 0 and z-axis projection M = 0. The set of Eqs. (8)- (10) are the BCS equations. The effective nucleon-nucleon interaction enters in the matrix element of Eq. (10) whose detailed expression is presented in Appendix B. As said above, we have adopted a finite-range interaction, and this allows us to use it, without any change, in both the HF and the BCS parts of our calculations [5].
Our HFBCS calculations give a description of the nuclear ground state in terms of the s.p. wave functions. The total energy of the even-even nucleus with A nucleons and Z protons can be expressed as: The density distribution of the system does not have any more spherical symmetry, but it depends also on the angular coordinates. In order to have an estimate of the non-spherical components of the nuclear density, we expand it in multipoles: where the terms of the density expansion are given by: In this last expression, we have used the Wigner 3j symbols [6] instead of the Clebsch-Gordan coefficients. We have considered the proton (p) and neutron (n) root mean square (r.m.s.) radii, which summarize the characteristics of the density distributions and are defined as: Here ρ α 0 (r) indicates the L = 0 multipole of the proton or neutron density, which is calculated by using Eq. (13) but restricting the sum on k to proton or neutron s.p. states only.
Nuclear deformations have been estimated by using the parameter which simplifies the comparison between nuclei with different size and number of nucleons. In the previous equation, R = 1.2 A 1/3 fm and Q 20 = 16π 5 Φ|r 2 Y 20 |Φ = 16π 5 dr r 4 ρ 2 (r) (16) indicates the quadrupole moment of the density distribution, with ρ 2 (r) the L = 2 term of the nuclear density, defined in Eq. (13). We have also calculated the charge radii, which are given by: The charge distribution, ρ charge (r), is obtained by folding the point-like proton density, ρ p (r), with the charge proton form factor. We have used a dipole parameterisation of this form factor [7], having verified that other, more accurate, expressions produce differences smaller than the numerical accuracy of our calculations.

III. DETAILS OF THE CALCULATIONS
The only physics input of our calculations is the effective nucleon-nucleon interaction. We have used a finite-range interaction of Gogny type, specifically the D1S parameterization [8]. In addition, we have considered the interaction D1ST2a [9], obtained by adding to the D1S force a tensor part of the form: In the above expression, τ (i) indicates the Pauli operator acting on the isospin of the i-th nucleon, and S 12 the usual tensor operator (see Eq.(A3)). The values of the parameters of the D1ST2a interaction are those of the D1S force in the common channels. For the parameters of the tensor part of D1ST2a, the values V T = −77.5 MeV, V T τ = 57.5 MeV and µ T = 1.2 fm, have been chosen to reproduce the experimental energy splitting between the 1f s.p. levels of the 48 Ca nucleus, in a HF calculation, and the empirical excitation energy of the first 0 − state of the 16 O nucleus, in a Random Phase Approximation calculation [9].
As indicated in Appendix A, we separate the contribution of the two space coordinates r 1 and r 2 by considering the Fourier transform of the effective nucleon-nucleon interaction. The required integrations in both coordinate and momentum spaces are carried out with the Simpson's technique. A good convergence of up to six significant figures is found by using a set of equally spaced points of 0.1 fm in r space and of 0.5 fm −1 in q space and, respectively, upper integration limits of 15 fm and 10 fm −1 .
The radial HF differential equations are solved by using the plane wave expansion technique described in detail in Refs. [10,11]. The iterative procedure stops when the total energies of two consecutive solutions differ by less than η = 10 −6 MeV. We have used this convergence benchmark in all our calculations.
In principle, after every iteration where the HF equations (4) are solved, the s.p. wave functions just obtained are used in BCS equations in order to modify their occupation probabilities. In practice, we have activated the BCS calculation only when the difference between the total energies of two consecutive HF solutions differ by less than a factor f η. We have found that, in practice, values of f between 500 and 1000 allow us to obtain a stable convergence of the solutions of the whole problem.
The iterative procedure starts by using the s.p. wave functions obtained by solving the Schrödinger equation for a deformed Woods-Saxon potential: where u = (r − R 0 )/a, l and s indicate the s.p. orbital angular momentum and spin operators, respectively, and V C the Coulomb potential. Even though the final result of the iterative procedure is independent of the starting set of s.p. wave functions, an appropriate choice of the values of the parameters U 0 , U so , a, R 0 is crucial to speed up the convergence. In our calculations, we have used the values indicated in Ref. [12], where they have been chosen to reproduce s.p. energies of the odd-even neighbouring nuclei. Pragmatically, we have found that values of |Λ| of the order of tens of MeV are needed to produce deformed solution. Specifically, we have always used Λ = ± 30 MeV. In all the calculations we have carried out, we observed that our procedure finds an energy minimum with the same type of deformation as that of the starting set of s.p. wave functions. For example, when we started with a prolate deformation, by setting Λ > 0, we found an energy minimum that maintained the prolate deformation, and viceversa. For this reason, in Sec. IV, sometimes, we show results obtained for both types of deformations.
Since the energies of the oblate and prolate solutions obtained for a given nucleus are different, we named optimal solution that with the smaller value of E(A, Z), in other words, the solution providing more binding to the system.
The relevance of the use of the new set of deformed s.p. wave functions has been studied by comparing the HFBCS results with those calculated within the HF+BCS approach of Ref. [4], where a spherical approximation is adopted. In this latter model, each s.p. state of angular momentum j is 2j + 1 times degenerated, the occupation of each s.p. state is equally distributed on all the possible z-axis projections and the full system conserves a spherical symmetry. For this reason we have indicated as spherical the results of the HF+BCS calculations. This is the most important difference, from the physics point of view, between the HFBCS and the HF+BCS models.
A second, and more technical, difference between the two approaches is in the treatment of the pairing. In HF+BCS we first carry out a HF calculation and, afterwards, we use the obtained s.p. wave functions to perform a BCS calculation. In HFBCS, the BCS calculations are inserted in the global iterative minimization procedure, connected to the solution of the HF equations.
We have also compared our HFBCS results with those of the deformed HFB calculations performed within a triaxial basis [3,13,14]. In these latter calculations, the solutions are expanded on a orthonormal basis of harmonic oscillators wave functions with oscillator length b 0 . The number of harmonic oscillators considered in each quantisation axis, (n x , n y , n z ), must satisfy the following condition for the energy truncation [14]: Here a x = (qp) 1/3 , a y = q 1/3 p −2/3 and a z = p 1/3 q −2/3 , where p = R y /R x and q = R z /R x are ratios between the semi-axes of the matter distribution. The results we have presented here have been obtained with N 0 = 9, and by using b 0 = 1.01 A 1/6 fm as indicated in Ref. [3]. By adequately choosing the values of p and q, a deformation shape can be selected. In our calculations we have used p = 1, that imposes the symmetry around the z axis, and q = 1.3 or q = 0.7 to select either a prolate or an oblate initial deformation, respectively. As it occurred in the HFBCS calculations, the deformation initially selected is maintained in the final solution provided by the HFB calculations.
We have tested the validity of our HFB results by verifying that the values of the total energies and density r.m.s. radii coincide with those presented in the compilation of Bruyères [15,16].
We close this section by mentioning what is known in the literature as the neutron gas problem [17]. The BCS calculations allow the presence of a long unphysical tail of the nuclear density distributions due to the contributions of slightly bound nucleons. We have shown in Ref. [18] that this is only a formal problem, since in actual calculations the numerical impact of these tails is irrelevant.

IV. RESULTS
In this section we present some selected results of our calculations with the goal of pointing out the combined effects of deformation and tensor force. We have considered 16 even-even Z = 34 isotopes, from 64 Se to 94 Se, and 10 even-even N = 34 isotones: 52 Ar, 54 Ca, 56 Ti, 58 Cr, 60 Fe, 62 Ni, 64 Zn, 66 Ge, 68 Se, and 70 Kr.

A. Comparison with HFB and HFBCS calculations
In the first step of our study, we have tested the reliability of our calculations. For this purpose, we have compared our HFBCS results with those of well established nuclear models. In our case, we refer to the deformed HFB results of the Bruyères compilation [15,16] and to those we have obtained by using the HFB approach of Ref. [14] described in Sect. III.
The total energies per nucleon, E/A, with E given by Eq. (11), of the nuclei considered are shown in Fig. 1. Our HFBCS results are indicated by the full circles and the benchmark HFB results by the empty squares. Since, as pointed out in Sect. III, in both calculations the deformation initially selected is maintained in the iterative procedure between the corresponding results are shown. The comparison between the energies obtained with the two nuclear models is very satisfactory. The relative differences ∆E HFBCS HFB are smaller than 1% for both type of deformations. We observe that most of the HFB energies in the prolate solutions (upper panel) are smaller than those of the HFBCS calculations. This is clearly indicated by the negative values of relative differences shown in the inset. The oblate solutions (lower panel) exhibit the opposite behavior.
We show in Fig. 2 the differences between the absolute values of the total energies per nucleon of the prolate and oblate solutions: Positive (negative) values of δ min indicate that the optimal solution has prolate (oblate) deformation. For HFBCS (full circles), these differences are, at most, 80 keV, in absolute value, and for HFB (empty squares) they are even smaller. These numbers are close to the numerical uncertainty of the calculations, which is of a few tens of keV. This fact has been named shape coexistence [19]. It means that, in practice, the optimal solution, and consequently the nuclear shape, is not so well defined in many cases.
According to these results, mean-field like approaches might not be enough to have a good description of these nuclei and dynamical, beyond mean-field correlations could be important. Within this context, angular momentum and particle number projections and fluctuations of the collective deformation parameters might play a relevant role [20][21][22]. The results of Fig. 2 also indicate that, often, the deformation of the optimal solution is not the same in HFB and HFBCS calculations. This happens for all the nuclei up to A = 62 in the N = 34 chain, and for quite a few Se isotopes. We have found that, in some cases, this also occurs when we compare our HFB results and those of Bruyères [15].
It is worth pointing out that the total energies of the optimal solutions obtained in HFBCS are slightly lower (0.5% at most) than those of the HFB for all cases considered, with the exception of the five nuclei with A = 62 − 66.
If the agreement between the total energies obtained with HFBCS and HFB is satisfactory, the situation is quite different when the deformation parameters are considered. We show in Fig. 3 the β 2 values obtained in HFBCS (full circles) and HFB (empty squares) calculations. Positive and negative values of β 2 refer to prolate and oblate deformations, respectively. The values of the β 2 of the HFBCS solutions are always remarkably smaller, in absolute value, than those found with HFB, although their overall behavior, as A varies, is similar in both calculations. The more relevant exception is that of the 62 Ni which is spherical in HFBCS and deformed in HFB, especially in the case of the oblate solution. On the other hand, the 54 Ca and 84 Se nuclei are clearly spherical in both approaches.
Finally, the effects of the deformation have been evaluated by comparing the total energies obtained in the HFBCS calculations with those calculated within the HF+BCS approach of Ref. [4], where a spherical approximation is adopted, as we have briefly discussed in Sect. III. The results of the two types of calculations are very similar, the relative differences being ±1.5% at most. In general, the deformed HFBCS produces optimal solutions that are more bound than the spherical ones. We have found only three exceptions: 54 Ca, 62 Ni and 64 Zn.
The impact of the tensor force on the HFBCS total energies is shown in Fig 4, where we compare the results obtained with the D1ST2a (empty circles) and the D1S (full circles) forces. Also in this case, we present separately the E/A values corresponding to the prolate (upper panel) and oblate (lower panel) deformations. We show in the  In the insets we give the relative differences between D1ST2a and D1S results, in percentage, calculated as indicated by Eq. (23).
insets the relative differences between the results obtained with the two interactions: In general, the tensor interaction generates less binding. The only exceptions we have found are those of the 54 Ca, for both deformations, and of the 84 Se, for the oblate solution. On the other hand, and in agreement with our previous studies [9], we notice that the role of the tensor force on this observable is rather small. We identify values of the relative differences which are 2% at most. These differences are, in any case, larger than those between the HFBCS and HFB energies obtained with the D1S interaction (see Fig. 1).
In Fig. 5 we show the differences δ min , defined in Eq. (22), calculated with the D1ST2a (empty circles) and the D1S (full circles) forces. As already pointed out, a change in the sign of δ min implies a change of the deformation of the optimal solution and one can see that this occurs in about half of the nuclei investigated.
The differences |δ min | between prolate and oblate energies are, in average, smaller when the tensor interaction is considered. We have obtained an average value of about 17 keV for D1ST2a against about 34 keV for the D1S. In the case of the HFB calculations with the D1S force (empty squares in Fig. 2), the average value found is ∼ 7 keV.
The effects of tensor force on the deformation are shown in Fig. 6, where the values of β 2 calculated with the D1S (full circles) and D1ST2a (empty circles) interactions are shown. Also in this case, we show separately the results obtained for prolate and oblate deformations.
The values of β 2 for the prolate deformation are almost the same for the two interactions. In the case of the oblate solutions, we have found noticeable differences in many of the nuclei investigated, the most noteworthy are those of the 62 Ni nucleus, which loses its spherical shape when the D1ST2a force is used, and of the 74 Se and 90 Se nuclei, which, on the contrary, become spherical.  Table I: Values of the deformation parameter β2 for the optimal solutions of the nuclei studied in the present work.
Up to now, we have presented our results separately for both prolate and oblate deformations. For the comparison with the experimental data we consider only the results of the optimal solutions. We show in Table I the corresponding β 2 values calculated with the D1S and D1ST2a interactions. We remark that, quite often, the type of the deformations of the optimal solutions for the two interactions are different. Out of the 26 nuclei considered, only eight, two of them are spherical, maintain the same type of deformation when the tensor terms are included in the interaction.
The total energies per nucleon corresponding to the optimal solutions obtained in HFBCS with the D1S (full circles) and D1ST2a (empty circles) interactions, are compared in Fig. 7 with the experimental data of Ref.
[23] (full black squares). We show in the inset the relative differences between our results and the experimental values: The differences with the experimental energies are quite small: they are, at most, 1% for the D1S interaction and reach values of about 2% for the D1ST2a. The HFBCS optimal solutions obtained with both interactions are less bound than the experimental values. From the point of view of the variational principle, this is quite reassuring.
The effect of the tensor force consists of a reduction of the nuclear binding. This has been already indicated for prolate and oblate solutions separately, and Fig. 7 shows that it is also true for the optimal solutions. A critical discussion of these results is in order, since the better agreement of the D1S results with the experimental data can be misinterpreted. The values of the parameters of the D1S interactions have been optimised to reproduce experimental binding energies and charge r.m.s. radii [8,24]. As said above, the tensor terms in D1ST2a have been added without changing the other parameters of the force. For this reason, it is plausible that the the inclusion of a new term in the interaction is worsening the overall agreement with the experiment.
In the Se isotope chain, the nuclei with the largest binding energy per nucleon are the 78 Se and 80 Se. Our results indicate that 78 Se is more bound than 80 Se in agreement with the experiment. In fact, 78 Se is the nucleus showing the smallest difference between our HFBCS (D1S) results and the experiment. In the N = 34 isotones the experimental |E(A, Z)/A| value is slightly smaller in 60 Fe than in 62 Ni while the opposite occurs in the HFBCS calculation. We have investigated the two-neutron separation energies for the N = 34 isotones and the two-proton separation energies for the Se isotopes. We have not found relevant effects due to the presence of the tensor force. On the other hand, the agreement with the experiment is good.
We present in Fig. 8 the values of β 2 of the optimal solutions obtained in the HFBCS calculations by using the D1S (full circles) and D1ST2a (empty circles) forces, and we compare them with some empirical estimations. In the upper panels the comparison is done with the values of the semi-empirical model of Möller et al. [25] (blue full squares). Even though no specific patterns are observed, our deformations are smaller, in absolute value, than those of Ref. [25]. Quite often, the deformation of the optimal solution obtained with the D1ST2a interaction has different sign than that of the D1S interaction, but also in this case we cannot identify any trend.
In the lower panels of Fig. 8 our HFBCS results are compared with the empirical data of Ref.
[23] (full black squares). In this case, we have considered |β 2 |, since the sign of the deformation in these empirical data is usually undetermined. Again, it is evident that our approach generates smaller |β 2 | values than the experimental ones. However, also this comparison should be considered with caution. The β 2 values of Ref. [23] have been obtained by assuming that the first 2 + excited state is due to a rotation of the deformed nucleus described with a sempi-empirical liquid drop model. The assumptions of this procedure are quite strong, and they lead to assign ground state deformations ground state deformations even to nuclei that are well know to be spherical: for example, β 2 = 0.353 is quoted for 16 O due to the presence of a 2 + state at 6.917 MeV.   In this section, we present the results of our HFBCS model concerning proton, neutron and charge density distributions, and their r.m.s. radii.
We have tested the reliability of our study by comparing our HFBCS results with those of HFB calculations carried out with the D1S interaction. The relative differences between the corresponding r.m.s. radii of the proton and neutron density distributions are smaller than 3% and 1.5% for the oblate and prolate solutions, respectively. On the other hand, the relative differences obtained with the D1S and D1ST2a interactions in the HFBCS approach are even smaller: at most 1.1% for the prolate solutions and 1.7% for the oblate ones. Contrary to what we have found for the total energies, in this case, the effect of the tensor is smaller than the differences between HFB and HFBCS results.
In order to study the effect of the deformation, we have calculated the r.m.s. radii with the HF+BCS approach. We have found that the largest relative difference with the HFBCS r.m.s radii of the optimal solutions is about 1%.
The situation is well summarized in Fig. 9, where we show the r.m.s. charge radii of the HFBCS optimal solutions of the nuclei studied, calculated, according to Eq. (17), with both the D1S (full circles) and D1ST2a (empty circles) interactions. These results are compared to the experimental values taken from the compilation of Ref. [26] (black full squares) and to those obtained in HFB with the D1S force (empty squares). These latter charge radii have been calculated by using the proton distributions of Ref. [16].
For the D1S interaction, the HFB radii are larger than those of our HFBCS calculations by 3% at most. Also the charge radii obtained with the D1ST2a force in HFBCS are slightly larger than those found with the D1S interaction, but the differences are smaller than 1.5%.
The comparison with the experimental data is limited to only seven nuclei. These few data are well described by all the three types of calculations, even though, globally, the best agreement is obtained for the HFBCS calculations with the D1S interaction.
We have analyzed in detail the density distributions and, as example of this study, we show in Fig. 10   for the two nuclei of our set of isotopes and isotones whose charge distributions are available in the compilation of Ref. [27]: 62 Ni and 64 Zn.
In the upper panels of the figure, we compare the HFBCS charge distributions with the empirical ones. The agreement between them is excellent at the surface, and this explains the good description of the experimental charge radii, which are mostly sensitive to this part of the distributions. Remarkable differences are evident in the nuclear interior, the region where correlations of various types, long-and short-ranged, are most effective [28]. The oscillations of the distributions obtained with the D1S force are smoothed by the presence of the tensor force, which produces charge densities closer to the empirical ones.
In order to frame the above discussion in a proper perspective, we remember that the empirical charge densities are tailored to fit elastic electron scattering cross sections. These experimental data have been measured within a restricted range of momentum transfer values, which, in our case, for both nuclei considered, goes up to q max = 2.2 fm −1 [27]. In the lower panels of Fig. 10 we show the elastic electron scattering cross sections calculated in Distorted Wave Born Approximation [29] by using the charge densities shown in the upper panels. We have assumed an incident electron energy of 300 MeV. In this kinematic conditions the value of q max is reached at θ max = 93.79 deg, which is indicated in the figure by the vertical dashed lines. As a consequence, the comparison between theoretical and empirical cross sections is meaningful only for θ < θ max , where a good agreement with the experiment is observed.  The effects of the deformation on the density distributions are related to the presence of terms with L > 0 in the expansion of Eq. (12). In Fig. 11 we show the L = 0 and L = 2 components of the proton and neutron HFBCS densities of 62 Ni and 72 Se. In each panel of the figure, we compare the results corresponding to the optimal solution obtained with the D1S force (solid curves) with the analogous components of both the prolate (dashed curves) and oblate (dotted curves) solutions obtained with the D1ST2a interaction. We indicate with a star the D1ST2a optimal solutions.
At the nuclear surface, all the ρ 0 distributions have, essentially, the same values. The differences show up in the interior, where the results obtained with the D1S interaction present larger oscillations than the other ones. In the case of the 72 Se nucleus the behaviour of the ρ 0 distributions remains the same, because the tensor force only produces a damping of the oscillations (see the dashed curves in Figs. 11e and 11g)). The case of the 62 Ni is more complex. The prolate D1ST2a solution shows proton and neutron densities similar to the D1S ones, even though the oscillations are damped. On the contrary, the densities of the oblate minimum, which is, by the way, the optimal solution, present a completely different trend (dotted curves in Figs. 11a and 11c).
The optimal solution found for the 62 Ni in the HFBCS calculation with the D1S force is spherical: the corresponding proton and neutron ρ 2 components must be multiplied by 100 and 20 to be seen at the scale of the figure (see Figs.  11b and 11d). When the tensor is added, the optimal solution becomes oblate, with β 2 ∼ −0.1. It is worth pointing out that the solution obtained with the D1ST2a interaction by starting with a set of prolate Wood-Saxon s.p. wave functions (dashed curves in the left panels) has spherical symmetry: it is necessary to use multiplicative factors of 5000 and 1000 in order to show the proton and neutron L = 2 terms at the scale of the figure.
For the 72 Se nucleus, the optimal solution is oblate, with β 2 ∼ −0.075, when the D1S interaction is considered, and changes its shape after including the tensor terms in the force, becoming a prolate nucleus, with β 2 ∼ 0.032. The ρ 2 obtained with the D1ST2a interaction for the oblate solution are of the same order of magnitude than those found with the D1S force, in particular the proton one that remains almost unaltered.

D. Single particle energies
As said above, in our model, the deformation is obtained by breaking the degeneracy of the s.p. states with the same n, l, and j, quantum numbers since states with different z-axis projection m have different energy. Because we assumed time-reversal invariance, those s.p. states with the same |m| remain degenerated. In order to study the  Table II. combined effects of the deformation and of the tensor force on the s.p. energies we have considered two quantities: the spread, Σ α nlj , and the centroid splitting, ∆Γ α nl . The spread is defined as the difference: According to this definition, the spread is zero for spherical nuclei, positive for prolate solutions and negative for the oblate ones.
We have found a strong relation between the spread of the s.p. energies and the deformation. This is evident in Fig. 12 where we present the s.p. energy spreads of the 1d 3/2 , 1d 5/2 , 1f 5/2 and 1f 7/2 states against the values of β 2 . These results are those of the optimal solution for each nucleus considered and have been obtained with the D1S (full circles) and D1ST2a (empty circles) interactions. The values for the nuclei of the N = 34 isotone chain (red symbols) refer to neutron s.p. states, while those of Z = 34 isotope chain (green symbols) stand for proton s.p. states. Finally, the straight lines fit the D1S (dashed) and D1ST2a (full) results separately.
It is evident that Σ nlj and β 2 are linearly correlated. This is confirmed by the large values of correlation coefficients obtained in the linear fits to the data and shown in Table II. Even more, this correlation is the same independently of the isotonic or isotopic chain analyzed. The slopes of the dashed and full lines are very similar, indicating a small effect of the tensor force. On the other hand, the lines fitting the data of the s.p. states with j = l + 1/2 (Figs. 12a and 12c) are steeper than those corresponding to j = l − 1/2 (Figs. 12b and 12d) Table II: Parameters of the linear fits of Σ nlj and ∆Γ nlj , as a function of β2, shown in Figs. (12) and (14). In both cases, the fitting function is y = a + bβ2. The uncertainties of the parameters and the linear correlation coefficients, r, are also given.
The second quantity that we have used in our study of the s.p. energies is the centroid splitting: where we have indicated with the centroid of the s.p. energies of the multiplet with quantum numbers n, l and j.
In Fig. 13 we show the values of ∆Γ α obtained in our HFBCS calculations for the 1d (upper panels) and 1f (lower panels) multiplets. The left panels indicate the results for the neutron s.p. states of the N = 34 isotones and the right panels those of the proton s.p. states of the Z = 34 isotopes. The results obtained with the D1S and D1ST2a forces are indicated by the full and empty circles, respectively.
The tensor force reduces the value of ∆Γ α . This behaviour is similar to the well-known effect that has been pointed out, discussed and explained for spherical systems by Otsuka and collaborators [30,31]. In that case, the tensor produces a reduction of the splitting between spin-orbit partners, which is precisely a quantity equivalent to the centroid splitting defined in Eq. (26) for the deformed nuclei. We have checked that this effect also occurs in the results of the spherical HF+BCS calculations that we have performed for all the nuclei here considered.
Since we have found a good correlation between s.p. energy spread and deformation, see Fig. 12, we repeated an analogous study also for the centroid splitting. We present in Fig. 14 the values of ∆Γ α against the deformation parameter β 2 . The results are those obtained for the optimal solutions of each nucleus considered, in the cases of the 1p (Fig. 12a), 1d (Fig. 12b), and 1f (Fig. 12c) s.p. states. We show the results for the neutron states in the case of the N = 34 isotones (red symbols) and for the proton states in the case of the Z = 34 isotopes (green symbols). Full and empty circles indicate the values obtained by using D1S and D1ST2a interactions, respectively. The data do not show any evident correlation with β 2 . The linear fits of the D1S (dashed lines) and D1ST2a (full lines) data remain almost constant against the changes of the deformation parameter. The absence of correlation between ∆Γ nlj and β 2 is quantitatively defined by the low values of the correlation coefficients given in Table II.
The tensor terms of the force modify the sequence of s.p. states in the deformed nuclei. As an example of this effect, we show in Fig. 15 the proton and neutron s.p. spectra of the optimal solutions of the 60 Fe and 90 Se nuclei obtained with the D1S and D1ST2a interactions.
The 60 Fe nucleus is oblate for D1S, with β 2 = −0.032, and prolate for D1ST2a, with β 2 = 0.050. In the case of the 90 Se nucleus, the optimal solution obtained with the D1S force is oblate, with β = −0.07, while for the D1ST2a interaction we obtain a value of β 2 = 8.7 · 10 −5 , indicating an essentially spherical shape.
The effect of the tensor is quite evident in 60 Fe. There is a slight increase of the spreading for the states with the same (n, l, j) values, and the inversion of the order of the levels with different |m|, due to the change of shape. The result of these two combined effects is that for the 1d states a level with j = 3/2, the 1d 3/2,1/2 state, has an energy lower than a state with j = 5/2, the 1d 5/2,5/2 level.
A said above, the optimal solution obtained with the D1ST2a interaction for 90 Se has a spherical shape. This is evident because all the s.p. levels with the same n, l, and j quantum numbers and different |m|, converge to an unique energy value. The deformed results obtained with the D1S force show, in some case, an inversion of the levels with different j by placing states with j = l − 1/2 below states with j = l + 1/2. This happens for the 1f 5/2,5/2 level whose energy is smaller than that of the three 1f 7/2,|m| states, and for the the 1f 7/2,7/2 level whose energy lies between those of the 1d 5/2,1/2 and 1d 5/2,3/2 levels.
Single particle states with their specific characteristics are the basic entities of our model. In order to have the possibility of comparing our predictions with some empirical observation we exploit the Koopman's theorem [3] which establishes that, in mean-field models, the global properties of odd-even nuclei are fully determined by those of the s.p. states of the unpaired nucleon.
By considering this approach, we have evaluated the angular momenta and parities of the ground states of those nuclei having one proton less than those considered in the N = 34 chain of isotones. For the case of the Z = 34      isotopes, we considered those odd-even nuclei with one neutron less than the even-even isotope partner. The angular momenta and parities of these odd-even nuclei are presented in Table III and compared to the empirical data taken from the compilation of Ref. [23]. The parentheses in some of these experimental assignments indicate that they are not fully identified. Since our model considers partial occupations of the s.p. levels, sometimes the definition of the Fermi level is quite ambiguous, and this generates uncertainty in the definition of the last occupied level, which, on the contrary, would be well identified in pure HF calculations.
By using the data of Table III we can better analyze the results of Fig. 15. The spin-parity of the ground state of the 59 Mn obtained in our calculation is 7/2 − while the experimental value indicated in [23] is 5/2 − , even though there are uncertainties on the spin assignement. In our calculations, both with and without tensor force, the energies of the proton 1f 5/2 states are always larger than those of the 1f 7/2 states, therefore our model does not account for that spin. From the neutron point of view, we remark that the empirical value of the 59 Fe nucleus given in [23] is 3/2 − , properly predicted by the calculation with the D1ST2a interaction.
The experimental spin-parity assignment of the 89 Se ground state is (5/2) − , correctly described by our calculation performed with the tensor force. Since the spin-parity of the 89 As is unkown, we cannot make a comparison with our predictions for the neutrons.
The results shown in Table III indicate that experimental spin-parity assignments are reproduced in only about half of the cases. Our calculations are unable to provide adequate results in all the Se isotopes up to 81 Se with the only exception of the 69 Se nucleus, for the D1S interaction. In nine of the nuclei investigated the results obtained with the D1S and D1ST2a interactions are different.

V. SUMMARY AND CONCLUSIONS
In this article, we have presented a model describing open shell nuclei. This model is based on the variational principle and uses Slater determinants built with s.p. wave functions whose radial part depends on m, the projection of the total angular momentum j on the quantisation axis z. This feature automatically introduces a deformation in the many-body state. Each step of the iterative procedure minimising the energy functional with these trial wave functions consists of two different calculations. In the first one, we solve the HF equations (4) generating the s.p. wave functions, and in the second calculation, we solve a set of BCS equations which modify the occupation probabilities of the s.p. states. Since the solution of HF and BCS equations is considered in each step of the minimisation procedure, we named HFBCS our model, to distinguish it from the approach of Ref. [4] where the solution of the BCS equations is carried out after the full solution of the HF equations has been found. We called HF+BCS this latter approach which, by the way, uses spherical s.p. wave functions. We consistently use the same finite-range interaction to carry out both HF and BCS calculations. We have considered the Gogny type D1S interaction [5] and an extension of it, the D1ST2a force, which also contains tensor terms [9].
The iterative procedure starts from trial wave functions which already have a deformation. In all our calculations we have found that the type of deformation is conserved until convergence is reached. This means that prolate or oblate trial wave functions lead to final results with the same type of deformation. We found the same feature in the HFB calculations of Ref. [14].
In general, the two, oblate and prolate, solutions, obtained for each nucleus considered, have very similar total nuclear energies. We have called optimal solution that with the smallest energy value.
The aim of our study was the investigation of how the deformation of the nuclear ground state emerged in our model, the effects on observable quantities, the effects of the tensor terms of the effective interaction and the, eventual, relation between deformation and tensor force. We presented results of our HFBCS model regarding energies, density distributions and single particle properties of medium-heavy nuclei belonging to the N = 34 isotone chain and to the Se, Z = 34, isotope chain.
Deformation effects on energies and radii are rather small. The total energies obtained in our HFBCS calculations are lower, at most of 1.5%, with respect to those of the spherical HF+BCS results. The differences between the r.m.s. radii obtained with these two different approaches are even smaller.
The effects of the tensor force are more evident. Calculations carried out with the D1ST2a interaction produce nuclei which are slightly less bound than those described by the D1S force. The relative differences between the results of the two calculations are smaller than 2% for all the cases considered, but the effect is clear and consistent in all nuclei investigated. This effect is worsening the agreement with the experimental energies which are smaller than those obtained without tensor by 1% at most. All these facts are in compliance with the variational principle which provides upper limits of the correct energy eigenvalues, and with the fact that the global fit to the experimental energies and radii carried out to select the parameters of the interaction has been done for the interaction without tensor.
Also the effects of the tensor force on the density radii are quite small; relative differences with the results obtained without tensor are smaller than 1%. These effects have the same sign in almost all the nuclei we have studied. Calculations carried out with tensor terms in the interaction produce r.m.s. radii larger than those obtained without them.
A detailed investigation of proton, neutron and charge density distributions has shown that the most evident differences between the results of the various calculations show up in the interior of the nucleus. Densities without tensor present a rather oscillating behaviour in the nuclear interior. The tensor is smoothing these oscillations and, for the cases where results are available, we found a better agreement with the empirical charge densities.
The s.p. energies are the quantities most affected by deformation and tensor force which, both, generate a reordering of the s.p. level scheme. The deformation destroys the 2j + 1 degeneracy of the spherical s.p. states characterised by the n, l and, obviously, j quantum numbers. We have assumed rotational symmetry around the z axis, and time-reversal symmetry, therefore, we obtain different s.p. energies for each value of |m|. Each spherical, and 2j + 1 degenerated, s.p. state is split in j + 1/2 different states. We have defined the spread as the difference between the s.p. energies of the two extreme states with the same j (those with |m| = j and |m| = 1/2) and we found a strong linear correlation between its values and those of the deformation parameter β 2 . This correlation is present for both oblate and prolate solutions obtained with or without tensor force.
The tensor force changes the type of deformation of the optimal solution, therefore the last occupied proton or neutron s.p. state. A comparison between the measured angular momenta of odd-even nuclei and those of the last occupied states does not show any specific trend and does not provide a real preference between calculations carried out with or without tensor. The only clear facts are that the two type of calculations produce different results in some of the cases analyzed, and only in half of the nuclei considered one the two calculations is able to predict the experimental values.
The parameter which better summarizes the information on the deformation is β 2 , defined in Eq. (15). Our calculations generate β 2 values remarkably smaller, in absolute value, than those obtained in HFB calculations. Also the comparison with the values obtained by an empirical model [25] and those indicated as experimental data [23] shows that our results are smaller, in absolute value. The size of the deformation is essentially the same for calculations carried out with and without tensor force, even though, in general, the optimal solutions with tensor force are less deformed than those without it.
The results of our study clearly indicate that the present accuracy of the experimental data on binding energies, charge radii and distributions imposes a new global fit of a force containing tensor terms in such a way that all the force parameters will be modified.
Our HFBCS approach proposes a peculiar manner to describe open-shell nuclei which automatically generates deformations in nuclear ground states. Under many aspects this approach is simpler than that of the HFB model, and shows s.p. properties that are still well recognizable. The extension of this approach to describe odd-even nuclei is under way. The set of s.p. wave functions with their occupation probabilities is the starting point to build up a Deformed Quasi-Particle Random Phase Approximation.
In the previous equations, the isospin term is given by: T p W = δ t k ,ti , p = 1, 3, 5 , 2 − δ t k ,ti , p = 2, 4, 6 . (A11) The contributions of the Coulomb force are, for the Hartree term, Here e is the elementary charge and we have indicated, respectively, with r < and r > the smaller and the larger values between r 1 and r 2 .
The two zero-range components of our interaction are the density dependent term and the spin-orbit one. The direct, Z ≡ U, and exchange, Z ≡ W, contributions of the latter are given by In the case of the density dependent term we have used for the function P (ρ) the expression: P (ρ) = ρ 0 (r 1 ) + ρ 0 (r 2 ) 2