A solvable contact potential based on a nuclear model

We extend previous works on the study of a particle subject to a three-dimensional spherical singular potential including a δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}–δ′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta '$$\end{document} contact interaction. In this case, to have a more realistic model, we add a Coulombic term to a finite well and a radial δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}–δ′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta '$$\end{document} contact interaction just at the edge of the well, which is where the surface of the nucleus would be. We first prove that the we are able to define the contact potential by matching conditions for the radial function, fixing a self-adjoint extension of the non-singular Hamiltonian. With these matching conditions, we are able to find analytic solutions of the wave function and focus the analysis on the bound state structure characterizing and computing the number of bound states. For this approximation for a mean-field Woods–Saxon model, the Coulombic term enables us to complete the previous study for neutrons analyzing the proton energy levels in some doubly magic nuclei. In particular, we find the appropriate δ′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta '$$\end{document} contribution fitting the available data for the neutron- and proton-level schemes of the nuclei 208\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{{208}}$$\end{document}Pb, 40\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{{40}}$$\end{document}Ca, and 16\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{{16}}$$\end{document}O.


Introduction
The present work generalizes the model introduced in a recent publication in which the energy levels of neutrons in doubly magical nuclei were studied [1]. The analysis carried out there is taken a step further by introducing a Coulombic term so that we can study the neutron-and proton-level scheme. Both models are based on the introduction of singular potentials supported on a sphere of radius R, indeed a δ-δ contact interaction which generalizes the Dirac δ. Within the Woods-Saxon model, the aforementioned interaction arises if we make the parameter that describes the surface diffuseness of the nucleus go to zero, which naturally leads to the description of the spin-orbit coupling by the Dirac δ potential. To improve the approximation on surface diffuseness, the δ interaction will be adapted to fit the available data for nucleon energy levels. Within this framework, we can avoid numerical methods for calculating wave functions and we can study the properties of the energy level scheme using the characteristics of the special functions that are the solutions of the corresponding Schrödinger equation. Another approach to obtain analytic solutions is based on the Nikiforov-Uvarov method, applied in [2] to obtain the neutron-level scheme. In this case, the approximation consists in the expansion of the centrifugal and Woods-Saxon terms of the potential in a series of exponentials and keeping the terms up to second order. a e-mail: cesar.romaniega@uva.es (corresponding author) The main advantage of these contact potentials is that they can often be resolved exactly, which gives a good insight of some quantum properties. This is why these interactions have a large number of applications in the modeling of real physical systems, some of which are listed in [1] and the references quoted therein [3][4][5][6][7][8][9][10][11][12][13]. In this regard, it is worth mentioning a recent work in which the vacuum interaction energy of a three-dimensional scalar field is calculated in the presence of two concentric spheres defined by the aforementioned contact interaction [14], so that the δ term leads to positive Casimir interaction energies, which are unattainable only with the Dirac δ.
In the present work, the δ term will be used as part of the phenomenological potential. It is worth noting that there are at least two different definitions of this interaction, so-called local and non-local [15,16]. These definitions are based on the construction of different self-adjoint extensions of the Hamiltonian without the singular term [17]. These extensions are initially defined for a one-dimensional Schrödinger equation but can be adapted to our problem with slightly modifications due to the spherical symmetry of the potential [1]. Here, we will use the local δ interaction, since it is compatible with the unambiguously defined δ potential when both interactions are supported on the same radius [17][18][19][20]. This compatibility condition, absent with the non-local δ , is an unavoidable requirement of the proposed model, as discussed in Sect. 2. In addition, the whole interaction defined in this way can be seen as a generalization of the Robin boundary conditions, which are obtained as a finite limit case [21]. Within this framework, the δ-δ interaction is introduced by imposing appropriate matching conditions for the radial wave functions at the radius R.
Although the formalism developed here is valid for any non-relativistic particle subject to a three-dimensional spherical potential consisting of a Coulombic term, a finite well, and a radial δ-δ contact interaction at the edge of the well, our main motivation is to analyze the already mentioned nuclear mean field potential. In fact, the Dirac δ interaction has been used for the calculation of resonant parameters and energy spectrum in [22], and it has also been employed to study the spectral function of the unbound nucleus 25 O [23]. Following the work [1], we now study the energy levels of protons and compare them with those of neutrons in the previous paper. As we will see, the general properties of the structure of the level scheme are maintained, although there are certain differences due to the Coulombic term.
The article is organized as follows. In Sect. 2, we introduce the potential under consideration, which is written to be well suited for future application to a nuclear model. The radial δ and δ terms appear here explicitly. In Sect. 3, we derive the secular equation for which the solutions give the bound states. In Sect. 4, we study the existence and localization of bound states, and the energy levels, and we will see that the general properties of spherical potentials hold with respect to quantum numbers. These findings are tested with available data of the nucleus 208 Pb in Sect. 5, where we explicitly compare the proton-level scheme with the neutron-level scheme. We also calculate the coupling of the δ interaction that best fits the data available for this nucleus. The same procedure can be applied to other doubly magic nuclei, and therefore, in "Appendix A" we provide useful data for all known doubly magic nuclei and in "Appendix B" the level scheme and the optimized value of the δ coupling for the isotopes 16 O and 40 Ca. The article ends with some final remarks.

Establishing the mathematical model for the physical problem
In this section, we present the model that we intend to analyze using a notation adapted to possible applications in nuclear physics, as we did in [1]. The starting point is the following single-particle three-dimensional Hamiltonian, written with respect to the coordinates of the center of mass of the system, where μ is the reduced mass and U 0 (r ), U so (r ), and U q (r ) are the Woods-Saxon potential [24], its first and second derivatives: The Coulombic potential is the static Coulomb potential of a uniformly charged sphere of radius R and therefore depends on the radial variable differently inside and outside the nucleus: where Z is the number of protons. As already mentioned, only neutron energy levels were studied in Ref. [1], so this Coulombic potential (3) was not included in that Hamiltonian. We now study the connection between the strength which appears in the radial interaction (2) and real physical parameters [25][26][27]. Note that the coefficients of the Woods-Saxon potential (V 0 ) and the spin-orbit term (V so ) must be positive if we want to reproduce the magic numbers available in the literature [28]. However, the sign and value of V q is not determined beforehand and we will use it to better fit the energy levels scheme, as discussed in Sect. 5. For this phenomenological potential, the nuclear radius is written as R = r 0 A 1/3 , being r 0 a constant and A the nuclear mass. Finally, the parameter a is related to the thickness of the nuclear surface, being a = 0.7 fm a typical value [29]. Following the same steps as in [1], first we substitute in (1) the Laplacian in terms of the orbital angular momentum L and the radial coordinates, so that, as usual, the eigenfunctions of the stationary three-dimensional Schrödinger equation are written as a radial part u n j (r )/r times an angular part Y jm (θ, φ), the latter being such that These Y jm (θ, φ) are linear combinations of spherical harmonics Y m (θ, φ), and are simultaneous eigenfunction of the operators L 2 , S 2 , J 2 = (L + S) 2 and J z [29]. The constants If = 0, we have j = 1/2 and this constant vanishes. Consequently, the radial part fulfills Hu n j (r ) = E n j u n j (r ), where If we now take the limit a → 0 + of f (r ) in (2), we get in terms of the Heaviside step function θ(x), the Dirac δ and its first derivative δ . Hence, the radial Hamiltonian (6) in the limit a → 0 + becomes To stress the mathematical rigor of the operations carried out previously, it must be taken into account that the δ-δ perturbation in (8) is defined using the formalism of self-adjoint extensions of Hermitian operators with equal deficiency indices (see the first Appendix of [1]), which gives different options for the meaning of δ term. As mentioned in Sect. 1, the one we adopt here is the local δ , which is defined by taking precise matching conditions at the radius that supports the interaction, see Sect. 3.3. Consequently, noting that R a, we can consider this singular one-dimensional Hamiltonian as a mean-field potential describing the energy levels of protons and neutrons. As we shall see in the following section, the main advantage of this approach is that the corresponding eigenvalue equation can be solved exactly for the wave function in terms of hypergeometric functions.

Solving the singular Schrödinger equation
In this section, we determine the eigenfunctions corresponding to the bound states of the singular Hamiltonian (8). Hereinafter, for simplicity and when no confusion arises, we will use as shorthand notation u ≡ u n j : where the Hamiltonian without the singular term is and we have introduced the real parameters α j , associated with the δ term, and β, associated with the δ term: The radial Schrödinger equation (9) is defined on the interval 0 ≤ r < ∞ and can be analyzed for a generic quantum particle subjected to a spherical well, with a δ-δ interaction at the edge and a Coulombic term. Next, we will obtain the wave function in 0 ≤ r < R and R < r , and then apply the matching conditions on r = R.

Wave function inside
For the study of solutions in the region 0 ≤ r < R, we take E > −V 0 and solve the differential equation First, we make the following changes to the parameters that appear in (12), and then, based on the asymptotic behavior of the solutions of (12), we propose the following change in the unknown wave function: With this, Eq. (12) is transformed into a Kummer differential equation: For each value of the orbital angular momentum , the general solution of (15) is where 1 F 1 denotes the Kummer function, A , B are arbitrary constants, and Note that the two solutions are linearly independent since c / ∈ Z for any . In fact, the second solution is not square integrable near zero due to the factor z 1−c . On the other hand, from the definition of 1 F 1 it is obvious that this function is equal to one at the origin, and therefore, it is square integrable on the finite interval considered. Consequently, the admissible solutions are precisely

Neutron case
Let us investigate now the limit Z → 0 as a consistency test of Eq. (18) that can be rewritten as [30] u (r ) = A r +1 e igr 2 /2R 2 ∞ k=0 1 2 where (x) n denotes the Pochhammer symbol. Consequently, we can evaluate On the other hand, the solution obtained in [1] for r < R was where J ν denotes the Bessel function of the first kind. From its definition as a power of series, (21) is equivalent to Therefore, the limit of u (r ) in (29) is essentially the result obtained in [1], as expected.

Wave function outside
To find the bound states in the region r > R is equivalent to solve the Schrödinger equation First of all, we proceed to make the change of variables With this, Eq. (23) becomes a Whittaker differential equation For each particular value of the orbital angular momentum , its general solution is given by where M χ,γ and W χ,γ denote the Whittaker functions [30], with C and D being arbitrary constants. Again, if we are looking for square integrable solutions, we need to know the asymptotic behaviors of these functions for large values of r , which are Consequently, the solution (26) is square integrable if, and only if, C = 0. In this way, the only possible contribution comes from the second term, so that

Neutron case
Again, we investigate the limit Z , λ → 0 as a consistency test of Eq. (28), giving On the other hand, the solution obtained in [1] for r > R was where K ν denotes the modified Bessel function of the second kind, which fulfills [30] 2z so it coincides with the limit just obtained in (29), as expected.

Matching conditions at the border
Once we have obtained the inner and outer solutions, we need to link both at r = R in a proper way. It is known [18] that there are specific conditions for u (r ) at the point r = R that fix a self-adjoint determination of the operator H r of Eq. (10), defining the key Hamiltonian of Eq. (9). These requirements can be written in terms of a SL(2, R) matrix as [21,[31][32][33] where f (R ± ) denotes the limit approaching R on the right or left, respectively. The function u (r ) was obtained in previous subsections: Eqs. (18) and (28). There is a rigorous discussion on the self-adjointness of H r with Z = 0 in Appendix A of [1]. In the latter, the appropriate domain in which the one-dimensional operator (10) with Z = 0 admits a four parameter family of self-adjoint extensions for a fixed value of is found. With this in mind, this Hamiltonian plus the δ-δ interaction was defined using the formalism of self-adjoint extensions of symmetric operators with equal deficiency indices. The same procedure is valid for Z > 0. To show that adding V C (r ) to the aforementioned Hamiltonian in [1] does not change the self-adjointness, we use the Kato-Rellich theorem. More precisely, since V C (r ) is bounded and Hermitian, we can state that the complete Hamiltonian is self-adjoint [34]. It should be noted that for β = ±2, the above matching conditions are no longer applicable, and in this case, the respective self-adjoint extensions of the radial Hamiltonian are given by the following boundary conditions at r = R [31] Using (18) and (28), the matrix relation (32) yields the following secular equation: where we have introduced the dimensionless quantities v 0 , w j , ε, and the constant g, already defined in Eq. (13): Then, previously handled variables can be rewritten based on these as We will denote the left side of (34) by ϕ (ε) and the right side by φ (ε) (which also depends on j and β), that is, Note that Im[φ ] = g since the Whittaker functions are real, and therefore, Im[ϕ ] = g for the left side of (34). In the situation, we are analyzing the first difference with the neutron case [1] is precisely the imaginary part that appears in the secular equation (34), although, as we will see, this fact does not alter the solutions obtained from the real part. In order to prove it, we first derive the following result on Kummer's functions.

Proposition 1
For any x, y, z ∈ R, the following relation holds Proof Using Kummer's first transformation formula [30], we have where w * denotes complex conjugation. Consequently, ln M − ln M * = −i z and Therefore, (2022) 137:33 Page 9 of 24 33 The proof is concluded noting that the logarithmic derivative of this function satisfies [30] 2i d dz In our case, taking the left side of the secular equation can be written as Finally, from Proposition 1: Consequently, the energy of the bound states will be determined by This equation does not admit solutions in closed form, but they will be analyzed in the next section.

General properties of the bound state structure
In the previous section, we have established the matching conditions that radial wave functions must fulfill so that the δ and the local δ interactions are well defined. From that, in this section we consider the complete Hamiltonian (9) and we will draw some relevant conclusions as to the existence and properties of bound states from the secular equation (34). In this secular equation, divergent terms appear for values of the dimensionless energy ε ∈ (0, 1) when either W λ, 1 2 + (η R) = 0 or when 1 F 1 (a, + 3 2 ; −ig) = 0. Therefore, it is necessary to analyze the zeros of these functions. Next, we will show that φ (ε), essentially a combination of Whittaker functions, has no singularities for ε ∈ (0, 1) because the Whittaker function in its denominator has no zeros, and we will also give a good approximation for the zeros of the Kummer function that appears in the denominator of ϕ (ε).

Zeros of the hypergeometric functions
First we want to show that φ (ε) in (38) is a continuous function when ε ∈ (0, 1), for which we must show that W λ, 1 2 + (η R) = 0 has no solution in the mentioned interval. Since the Whittaker function can be written in terms of the Tricomi function as the zeros of W χ,γ (z) are those of U 1 2 + γ − χ, 1 + 2γ, z . To prove the absence of positive zeros we use the first theorem of [35], which states that there are no positive zeros z ν of U (a 1 , a 2 ; z ν ) = 0, if the two real parameters satisfy 2a 1 − a 2 > −1. In the present case, we have W λ, + 1  ∈ (0, 1).
On the other hand, the function ϕ (ε) defined in (37) will not be continuous in the interval ε ∈ (0, 1) due to the zeros of 1 F 1 (a, + 3 2 ; −ig), the dependence on ε entering only through a defined in (36). As a good approximation to the a-zeros of this function, 1 F 1 (a n , c; z) = 0 is given by [30] even for low values of n. In the present case, c = + 3 2 , z = −ig, and a tedious calculation takes us from (36) and the truncated form of (48) to the following definition Therefore, the function ϕ (ε) will present a finite number of discontinuities in ε ∈ (0, 1), which are approximately given by ε n . The typical behavior of the real parts of the ϕ (ε) and φ (ε) functions is shown in Fig. 1, where the ε values at the intersections of the blue curves with the orange curve correspond to the solutions of the secular equation (34). The monotonicity properties of the curves shown in Fig. 1 are due to the fact that The first inequality holds for all the intervals delimited by the zeros of the function 1 F 1 (a n , + 3 2 ; −ig) = 0. We have shown the inequalities of (50) for the particular case g = 0 in Appendix B.1 of [1], where we have also proved that ϕ (ε) is one to one and onto as a function being I n the interval delimited by two consecutive zeros of 1 F 1 (a, + 3 2 ; −ig). With the approximation given in (49), we have I n ≡ (ε n+1 , ε n ), the interval (ε 1 , 1) not being considered, see Fig. 1. In all the examples shown in this work, particularly those necessary for the nuclear model, these properties have been verified (although not rigorously demonstrated). From now on, we will assume that the two inequalities in Eq. (50) are valid. For simplicity, in this case we have used the same value of w j for all in order to illustrate the dependence of ϕ (ε) on . Although this could be valid for a generic quantum particle subject to the singular potential described, note that the dependence of w j on and j should be considered for nucleons within the proposed phenomenological model, as discussed in Sect. 5

Angular momentum dependence
The objective of this section is to study some properties of the number of bound states with angular momentum , denoted by n . Strictly speaking, the total number of bound states is given by (2 + 1)n when we take degeneracy into account. As shown in Fig. 2, n +1 ≤ n . In addition, there appears to be a certain max such that n max > 0 and n = 0 ∀ ≥ max + 1.
This behavior is typical of spherical potentials [36] and was proven for g = 0 in [1]. In our case, we will use the a-zeros given in (48). Specifically, we analyze the behavior of ε n with respect to . We can state that for any ≤ max the number of bound states is given by where M is the number of finite discontinuities of ϕ (ε) for ε ∈ (0, 1) and This result, analogous to Theorem 1 in [1], can be proved using the approximation ε n in (49). First of all, note that the cases ϕ (0) = φ (0) and ϕ (1) = φ (1) are excluded since we are considering energies ε ∈ (0, 1). Due to the properties of ϕ (ε) and φ (ε), the number of bound states is primarily determined by the number of discontinuities. For a fixed , these are given by the last value of n such that ε n > 0. Therefore, solving ε n = 0 we obtain where we have assumed −24g 2 + 9( + 2 + v 2 0 ) ≥ 0 and · denotes the floor function so ε n ≥ 0. Note that ϕ (0) > φ(0) implies the existence of an additional bound state whose energy is the closest to ε = 0. On the other hand, if φ (1) < ϕ (1), there is one less bound state.
The above result is based on approximation (48). In Table 1, we compare approximation (49) with the value obtained by numerically solving 1 F 1 (a(ε), + 3 2 ; −ig) = 0. In Table 5 of "Appendix A," we have included similar information for most of the known doubly magic nuclei, verifying that ε n gives a correct approximation. Furthermore, we have verified for these nuclei that the inequality − 24g 2 + 9( assumed in the previous proof, is also true. In Eq. (53), we have assumed the existence of an upper bound for the angular momentum, max . Bargmann inequality [37] n < ensures the existence of this upper bound for spherically symmetric potentials satisfying As explained in [1], due to the δ interaction, it is not clear how to interpret (57) and Bargmann result cannot be used. In any case, using ε n we can give an approximate value for the maximum value of the angular momentum. Specifically, we investigate the first zero for a given value of the angular momentum, that is, we solve ε 1 = 0 for . We got max = max 6 − 6π 2 + 2 9π 2 − 36 v 2 0 + 9 2π 2 + 1 − 24 π 2 − 4 g 2 3 π 2 − 4 , 0 .
This result is analogous to Theorem 2 of [1]. It is worth noting that an increase in v 0 , V 0 , and R tends to increase max . However, if we boost g, the maximum value of the angular momentum diminishes due to the repulsive character of the Coulombic potential. The existence of max was proved for g = 0 in [1]. For the previous example, illustrated in Fig. 2 and Table 1, we have max = 13.34 = 14. As shown in Fig. 3, in this case the approximation is exact. Note that we have assumed that the inequality holds, which is verified in Table 6 of "Appendix B" for all doubly magic nuclei. Now, let us focus on the order of the bound states. For spherically symmetric potentials, it can be proved using spectral properties of the Hamiltonian [36] that the energies of bound states characterized by n and satisfy We extend this result to our case, in which the additional quantum number j is present and, from now on, will be displayed using nuclear spectroscopic notation n j ≡ n j . If there are bound states with relative energies ε n j , ε (n+1) j , and ε n( +1) j for n, ∈ N 0 , the following inequalities hold: a result that was proved for g = 0 in [1]. Note that the first inequality is trivial. The second one, which only applies to j = + 1/2, follows from the behavior of ϕ (ε) and φ (ε) with respect to . Firstly, if we have 1 = and 2 = + 1, we need j 1 = 1 + 1 2 = + 1 2 and j 2 = 2 − 1 2 = + 1 2 . Consequently, since V so > 0 [29]. Now, we study how ϕ (ε) and φ (ε) vary with respect to . The generic behavior is illustrated in Fig. 4 for a particular case. We analyze the points of intersection between ϕ (ε) and y = y 0 > 0 for two values of : Note that for all the proton cases considered in this paper we have φ (ε) > 0 in ε ∈ (0, 1). This also holds for the nucleus 132 Sn used in [1]. Due to the monotonicity and the behavior shown in Fig. 4a, we conclude that < implies that ε > ε . Now, taking into account of Eq. (63) and that φ +1 (ε) > φ (ε) for any and ε ∈ (0, 1), see Fig. 4b, it follows that the inequality −ε n j < −ε n( +1) j is necessarily fulfilled. Note that in Fig. 4b we assumed a constant w j . For 2 , w 2 j 2 >0 and for 1 , w 1 j 1 < 0, then the quantity φ +1 (ε) − φ (ε) is even greater than that shown in Fig. 4b. Finally, the third inequality in (62) is easily proved considering that the dependence on j only enters through w j /(2 − β) 2 , being w j − > 0, w j + < 0, where we are considering > 0 and j ± ≡ ± 1/2.

Low-lying energy states
As it can be seen from the examples shown in Figs. 1 and 2, ε n can be reasonably well approximated by the value of the bound state, specially for the first values of n and : We must bear in mind that this approach only considers the a-zeroes of 1 F 1 (a n , + 3 2 ; −ig), so it does not depend on w j , that is, on j. However, this dependence can be introduced by considering the complete secular equation ϕ (ε) = φ (ε) so Eq. (65) can be used to give an easy and analytical proof of the inequalities (62): • The inequality −ε n j < −ε (n+1) j can be shown from (65), as it is clear that ε n+1 < ε n since the term decreases with n. Consequently, since we are considering ε n ε n j , the inequality is proved.
• The inequality −ε n j < −ε n( +1) j is also trivial from Eq. (65) since the terms depending on are and this quantity decreases with . • Let us consider finally the inequality −ε n +1/2 < −ε n −1/2 . As mentioned before, the dependence on j only enters through w j /(2 − β) 2 , being w j > 0 for j = − 1 2 and w j < 0 for j = + 1 2 . Since we are considering > 0, it is clear that we can write Note that the approximation of Eq. (65) could be no longer valid if w j 0 such that φ (1) ≈ 0.

Proton and neutron energy levels of 208 Pb
This section is dedicated to applying the above framework to a realistic physical situation. To do this, we choose the free parameters values of our potential that have been shown in order to mimic certain nuclei. In particular, for this phenomenological potential the usual parametrization is [29] where the + sign is for a proton and the − sign for a neutron. This model is particularly useful for doubly magic nuclei. Consequently, we will exemplify the procedure with the nucleus 208 Pb. In "Appendix B," the same procedure is applied to 40 Ca and 16 O. We have illustrated with Pb since it has the largest number of bound states. Our objective is to show that the low-lying bound states-level scheme obtained from our model is quite reasonable for the isotopes considered. A high-precision numerical fit is beyond the purpose of the model. To begin with, in Table 2 we show the physical constants and dimensionless quantities needed for these three nuclei. The parameter β, the intensity of the δ term, is not fixed and in fact we will use it to adjust the results to fit the data available for the above isotopes.
Next, for 208 Pb we determine the neutron and proton energy levels for = 0, 1 setting β = 0. In Fig. 5, we show the points of intersection in this case.
Using nuclear spectroscopic notation, we compare the data obtained here with the neutron and proton energy levels found in [38], where the usual parameters of the Woods-Saxon (WS) model are optimized to obtain the best fit with the available experimental data. The results are shown in Table 3.  Table 3   Table 3 Proton-and neutron-level scheme of 208 Pb arising from our model with β = 0 compared to that of the optimized Woods-Saxon model developed in [38] State We also observe that the energy values decrease with n. In general, for all the isotopes considered in this work we have |E(n −1/2 )| < |E(n +1/2 )|, as expected from Eq. (62). It should be mentioned that this is the opposite behavior to that found for electrons. In fact, the nucleon that is aligned with the orbital angular momentum is attracted more strongly. The explanation of all the experimental magic number was finally possible when the spin-orbit interaction was included in the mean-field shell model. It is also worth noting that the energy values of neutron are deeper than those of protons as shown in Fig. 5. This can be partly explained by Coulomb's repulsive term in the potential. However, since N > Z for this nucleus, note that V 0 is larger for protons, as can be see in (69). The latter is the result of the average nucleon-nucleon potential, which is stronger for the proton-neutron interaction than the neutron-neutron or proton-proton case. Consequently, V so is also different for protons and neutrons. For other nuclei without a neutron excess, 40 Ca or 16 O, protons will not feel a stronger potential well than neutrons. Now we are going study the effect of the δ term for 208 Pb. From a phenomenological point of view, it can be used to improve the level scheme given by our model with just the Dirac δ. We denote by E β (n j ) the energy of the state characterized by {n, , j, β}, being E WS (n j ) the corresponding energy arising from the optimized Woods-Saxon model in [38]. To measure the amount of dispersion for this set of values, we define where N is the number of states. In Fig. 6, we have represented this quantity as a function of β and it is obvious that there are values of β that minimize this variation. We have focused on the values that minimize σ β and the results are shown in Table 4. Note that the values β = 0 and β = 1 correspond to V q = 0 and V q = −20.83 MeV fm 2 , Table 4 Proton-and neutron-level scheme of 208 Pb arising from our model with β = 0 and the optimal β compared to that of the optimized Woods-Saxon model developed in [38] State respectively. It is clear that the introduction of a non-trivial β allows us to better reproduce the level scheme found with the Woods-Saxon model.

Concluding remarks
In this work, we have extended a previous study on a nuclear model based on an interaction due to a three-dimensional Dirac δ and local δ , allowing the presence of the Coulombic term due to the electromagnetic interaction between protons. From a mathematical point of view, this contact potential has been precisely defined using suitable matching conditions satisfied by the radial wave functions. We have obtained general properties relative to the number and behavior of bound states of the aforementioned nuclear potential. As expected from what we know about spherically symmetric potentials, the eigenstates are less bound if we increase n or the angular momentum, the other quantum numbers being fixed. We have also shown that there is a maximum value of , so that consequently there is a finite number of bound states, since the number of them, for a given value of , is finite. We have given an estimate of these two quantities, max and n . We have analyzed the effect of the spin-orbit interaction, obtaining that the state with j = + 1/2 is more bounded than the one with j = − 1/2, as expected from considerations on the shell model. For the lowest energy states, those near the bottom of the well, all these properties have been analytically tested. In the last section, we used this model to approximately describe the single-particle neutron and proton energy levels of the doubly magic nucleus 208 Pb. We have shown that the Hamiltonian (8), based on the Woods-Saxon potential with the surface diffuseness going to zero, provides a fair approximation for the low-lying bound states. The δ term gives the nuclear spin-orbit contribution. The aim of the additional interaction, given by the local δ , is providing a correction such that the results of the proposed model better fit to the available data. In these sense, we have determined which is the best value of the δ coupling, β, for the nucleon and proton energy levels in the three nuclei 208 Pb, 40 Ca and 16 O. In particular, we have focused on 208 Pb, comparing our level scheme for nucleons with a Woods-Saxon model in which the free parameters as a and r 0 are not the usual, but those optimized in order to fit with the available experimental data. As we have seen, the introduction of the δ interaction does not change the fundamental properties of the level scheme. This interaction is initially justified for the construction of a self-adjoint extension of the Hamiltonian, but, as we have seen, it can be used to improve numerical precision. This justifies the choice of so-called local δ interaction. Note that the other possible choice for the δ interaction is not compatible with the Dirac δ if both are supported at r = R. We have proved that the latter is a requirement when the starting point is Woods-Saxon potential.
This simplified model could be used to gain insight into the proton and neutron energy levels for doubly magic nuclei. In this sense, the case studied in [1] could be seen as a particular case, Z = 0, of the present model. In both cases, the main advantage over the Woods-Saxon potential is that the eigenfunction equation can be solved exactly, obtaining analytic properties of the spectrum using well-known features of special functions: for Z = 0 the solution is given in terms of hypergeometric functions with complex arguments, while for Z = 0 the solution reduces to the usual Bessel functions. In the latter, it is easier to take advantage of known results in order to prove the desired properties of the bound state structure, as was done in [1]. However, hypergeometric functions with complex arguments are not as well studied in the scientific literature, so a rigorous analysis of the Z = 0 case has not been possible and a more qualitative study has been carried out instead. Both depend on the dimensionless quantities v 0 and g, which are also calculated for these nuclei. Note that both radicands r 1 ≡ −24g 2 + 9( + 2 + v 2 0 ), r 2 ≡ 9π 2 − 36 v 2 0 + 9 2π 2 + 1 − 24 π 2 − 4 g 2 , (72) must be positive, which is in effect shown in Table 6 for all isotopes considered here. and neutrons is the same for these two isotopes, so we use the same potential depth, V 0 = 51 MeV, for both level schemes. From this table, we can also obtain the dimensionless quantities necessary for the calculation of the bound states: 40 Ca: v 0 ≈ 6.797, g ≈ 1.730 and w j = − 8μR h 2 V so ξ j ≈ −18.723 ξ j , 16 O: v 0 ≈ 5.008, g ≈ 0.939 and w j = − 8μR h 2 V so ξ j ≈ −13.795 ξ j .
In this case, we will compare with the data obtained from the usual Woods-Saxon model [29] and not with an optimized one as discussed in Sect. 5. We can evaluate σ β , which was defined in Eq. (71) and is represented in Figs. 7 and 8. In both cases, σ β is similar for protons and neutrons. Unlike 208 Pb, the best value of β for neutrons is obtained before the peak at β = 2 and after this peak for protons. Remember that for the values β = ±2, the matching conditions (32) are ill defined and we have to use the conditions (33).
It is observed that, as happened in Sect. 5, negative values of β give rise to greater dispersion. Once we have obtained the best values of β, we can compare the resulting level scheme with the one in the absence of the δ term (β = 0) and with the one that arises from the regular Woods-Saxon model, which is shown in Tables 7 ( 40 Ca) and 8 ( 16 O).
As with the previous results, note that for 208 Pb it is also possible to test the above inequality for n = 1, see Table 3. This is why, we have used this nucleus to test our model based on the δ-δ contact interaction in Sect. 5.