Discontinuous radial potentials, bound states and Dirac equations

Discontinuous, radial ‘square-well’ potentials of scalar and (time-component) vector types require special quantization conditions in Dirac equations transformed to second-order differential equations. Quantization conditions for positive energies are derived from first-order and second-order differential equations using amplitude-phase methods. Analytic condition for no orbital angular momentum of the large spinor component is obtained. Numerical applications relevant to independent-nucleon shell models illustrate level shifts due to shape (square-well/Woods–Saxon) differences.


Introduction
All components in a Dirac spinor are involved in defining energy spectra. Hence, in separated second-order differential equations the spinor components must have a common set of bound states. Separated component equations for scalar-and vector type potentials may look quite different, in particular for discontinuous potentials; see e.g. [1]. These issues may have some relevance for studies of levels in an ideal shell model of nuclei, consisting of noninteracting nucleons in a common potential [2,3]. Ideal shell models of nuclei are based on either second-order (Schrödinger) equations [4] or Dirac equations [5,6], and a square-well approximation of a Woods-Saxon potential must be treated differently in either case.
Constant (square-well) potentials often provide analytic equations and may simplify some basic understanding of quantum mechanics. This is explored in many text books introducing non-relativistic quantum mechanics [7]. In relativistic quantum mechanics most classical text books, if not all, focus on Coulomb potentials [8][9][10][11], pointed out by several authors in recent publications [12]. At present, the 1+1D case of square wells has attracted most interest. Often these studies are very general as to positive and negative energy levels [12]. The present study is modest in this respect, but introduces new methods to study radial discontinuous potentials in 3+1D with second-order Dirac equations.
The amplitude-phase methods introduced here [1] apply to first-order radial Dirac equations of bound states, as well as to second-order radial Dirac equations which require special conditions for discontinuous potentials. An analytic part of the present study is restricted to ground states with vanishing orbital angular momentum. Generalizations to other angular momentum quantum numbers, involving special functions like the Riccati-Bessel functions, should be straight forward.
The famous 'log-derivative' approach (involving the Riccati equation) is trivial to apply in regions where the quantum wave is exponentially behaved. In regions where the wave is oscillating and having zeros the method becomes non-trivial, needing a special algorithm [13]. A comparison of most existing numerical methods is published in [14].
The present amplitude-phase approach is easy to apply and to understand in an oscillating region [15]. It applies to neighboring non-oscillating (classically forbidden) regions as well [16]. A basic idea of the amplitude-phase method is to define a 'slowly varying' amplitude function, which in turn defines the (wave) phase function by quadrature. The complete wave function is formally exact. Indeed, an amplitude of an arbitrary quantal wave is not unique; only a harmonic wave has a well defined (constant) amplitude. The present amplitude-phase method defines a relation between an amplitude function and a phase function in such a way that the complete wave is exactly the same with differing boundary specifications of the amplitude function. Efficient numerical performance is achieved for almost constant, nonoscillating amplitude functions. Historically, the amplitudephase method for calculating bound states originated in the analysis of (second-order) radial Schrödinger solutions [15].
Applied to first-order Dirac equations the amplitudephase approach uses slowly varying complex amplitude functions, and differ from an approach applied to transformed second-order differential equations using slowly varying real amplitude functions.
In this study, the nuclear energy (E) region of interest is restricted by the rest mass energy, i.e. E mc 0 2 < < , where m is the bounded neutron mass and c is the speed of light. This restriction is chosen to make a closer connection between non-relativistic and relativistic models of radial square-well bound states. The centre of potentials is considered fixed. Equations are made dimensionless by relating lengths to the relevant Compton length and energies to the rest mass energy of the neutron. The potentials are chosen to be relevant for singe neutron-nuclei states. Several amplitude-phase formulas are applied and for the simplest ground states an analytic condition is derived.
In section 2 the first-order radial Dirac equations are presented together with an amplitude-phase method for bound states. The approach is applicable to any shape of the vector/ scalar potentials, discontinuous or not. Section 3 introduces the second-order Dirac equations. Here an analytical solution for l=0 is explored and new quantization conditions in case of discontinuous potentials are presented. Numerical results are discussed in section 4. Conclusions are given in section 5.

First-order radial Dirac equations
The Hamiltonian  defining the time-independent 4×4 Dirac equation for a fermion mass m in terms of its 3-momentum operator p is given by [8]: Here c is the speed of light and s being the vector of the three 2×2 Pauli matrices σ j , j=1, 2, 3 (see [8]). Symbols 2  and 2  are unit and zero 2×2 matrices, respectively. In (2.1) the 4×4 potential  is  being a unit 4×4 matrix, V a time-component radial vector potential and S a radial scalar potential. For central scalar and vector potentials the radial Dirac equations reduce after separating out the angular dependence to two pairs of coupled first-order ordinary differential equations, one pair for each orbital angular momentum quantum number l, given by [5,9] F r r Both potentials vanish as r  +¥. For positive energies E the large spinor component is F. The single particle of mass m satisfying (2.4) has rest energy mc 2 . According to [9], the non-zero integer (Dirac spin) parameter κ in (2.4) is related to the orbital quantum number l in such a way that l 1 The Dirac components F(r) and G(r) can be chosen to be real-valued functions of real r. In equations (2.4) both Dirac components are assumed to vanish at r=0 and at r = +¥.

Dimensionless Dirac equations
The second-order Dirac equations are reduced with the length and energy scales In dimensionless form the first-order equations for F and G with a given value of κ are given by As x  +¥ with vanishing potentials, one has

Amplitude-phase approach
In the present section the amplitude-phase approach is valid for any shape of the well, discontinuous or not. When approaching the original first-order radial Dirac equations the fundamental amplitude-phase solutions F x  ( ) ( ) and G x  ( ) ( ) for analyzing the components F(x) and G(x), respectively, are expressed in exponential form according to [1] (the xdependence is implicit): where the amplitude components u F,G are complex-valued functions of x, and f is a real-valued phase function satisfying By substitution of the amplitude-phase Ansatz (2.10) into equations (2.8) and (2.9), the amplitudes u F,G turn out to satisfy the nonlinear coupled equations In the exterior region of a discontinuity, or more generally as x  +¥, one has The spinors are expected to be oscillatory in a single-well region of an effective potential. The same spinor amplitude functions u F G , from (2.12) are continued into classically forbidden regions, where they are expected to to grow indefinitely. Thereby the phase contributions in the Ansatz (2.10) decrease to finite values as x  ¥.
Equations (2.12) are exact, but not always tractable by numerical methods. In order to find slowly varying amplitude components one starts the numerical integration at some xvalue fulfilling approximate conditions of slowly varying amplitudes. For a radial square well being discontinuous at x=x 0 , one choses the potential value at the oscillating side. In case of a smooth well the integration starts at an x-value 'x m ' satisfying (2.12) with the left hand side being zero. This is assumed also for discontinuous wells, where x m =x 0 . Hence, one solves the approximate eigenvalue problem at some, in general, unspecified x=x m : The requirement that the determinant of the coefficient matrix in (2.15) vanishes implies that When considering x m as an independent variable one can define a function The function K 2 (x) is assumed being positive in some (single) region, so that K(x)>0 exists. K 2 (x) can be understood as being related to an 'effective' semiclassical kinetic energy. By subtracting the asymptotic constant value K 2 +¥ ( ) from K 2 (x) one then obtains something related to an 'effective' interaction potential. A non-relativistic approximation , where ò≈1, suggests that an effective potential near non-relativistic con- . A generalization to accommodate relativistic effects is an effective potential defined by The energy dependence is small if energies and potentials are almost nonrelativistic. The centrifugal term x in the separated second-order differential equations of the components. Note that (2.17) is used only for initiating the integration of (2.12). It is used to specify x m by locating the single maximum of K 2 (x), which easily can be done numerically. In that case the derivative of K 2 (x) with respect to x is zero at x=x m and should ensure locally adiabatic component solutions u F and u G there.
From (2.15) and (2.16) the amplitude components must satisfy: The positive normalization factor N is determined from condition (2.11), yielding The (complex-valued) amplitudes u F and u G are finally determined by integrating the differential equations in two directions from x x m = . This position defines the phase reference point so that x 0 . The Dirac components satisfying (2.8) and (2.9) can be chosen as real functions. To construct a real solution for positive energies, a complex amplitude component, say u F , is written as so that a real representation of the 'large' Dirac components be given by together with its relation to the small component partner is not important for the amplitude-phase quantization condition, but the phase f(x) is to be specified by an additional constant so that F satisfies the boundary conditions. Assuming

2.25
The non-negative integer n is referred to as the radial quantum number for the F component. Applications of the condition (2.25) are presented in tables 1 and 2. If, instead, u G is used to derive a condition similar to (2.25) and a real-valued G component defined by a corresponding expression (2.22), then the large component F would be constructed from the relation

Second-order decoupled Dirac component equations
The general first-order Dirac equations (2.4) are transformed to Schrödinger-type second-order differential equations by substituting It is not obvious from first inspection that the two differential equations in (3.2) have the same energy spectrum. A spinsymmetry assumption S=V simplifies R f but not R g . A pseudospin symmetry S=−V simplifies R g but not R f . For step-like potentials the trailing terms of R f and R g in (3.3) respectively (3.4) may develop delta function behaviors at isolated r-values.  in units of MeV, and increasing radial node numbers n of the F component, are shown for two scalar potentials. Energies a c ,  correspond to a constant radial potential well with v 0 =0, s 0 =−0.055 and x 0 =26. Entries ò a are from (2.25), (3.35), and (3.36), while ò c are the Klein-Gordon energies. Energies ò b correspond to a Woods-Saxon potential with the given parameters and with thickness a x 0.09 0 = , and are obtained from (2.25). For κ=−1 condition (3.23) is also used to obtain 1 a  -. The rightmost column shows the order of relativistic levels as the energy increases.  in units of MeV, and increasing radial node numbers n of the F component, are shown for two mixed scalar/vector potentials. Energies ò a,c correspond to a constant radial potential well with v 0 = 0.14, s v 0.055 0 = -and x 0 =26. Entries ò a are from (2.25), (3.35), and (3.36), while ò c are the Klein-Gordon energies. Energies ò b correspond to a Woods-Saxon potential with the given parameters and with thickness a=0.09x 0 , and are obtained from (2.25). For κ=−1 condition (3.23) is also used to obtain ò a −1. The rightmost column shows the order of relativistic levels as the energy increases. In the present study of constant potentials discontinuous at x=x 0 one obtains simplified differential equations in two separate regions. From (2.8) and (2.9) one obtains: where now v and s have constant values. If F G , ¢ ¢-discontinuities are ignored, (3.5) and (3.6) have different sets of energy eigenvalues for a given κ. Note that for spin-less particles (ignoring the small component G), equations in (3.5) would correspond to a Klein-Gordon equation with a radial square-well potential [8].
The potentials enter in quadratic expressions in (3.5) and (3.6). An 'effective' potential W x , F  ( ) for F may be defined as in (2.18). By denoting the coefficient in (3.5) by R F the effective potential is Effective potentials for κ=−1 and κ=−3 are illustrated in the left respectively right subplots of figure 1 for a radial square-well potential and a Woods-Saxon potential, both defined in section 4. The energy dependence of the effective potentials in the figure is specified by the respective groundstate energies.

Discontinuous state derivatives
When derivatives F x d d and G x d d are discontinuous at a radial distance x=x 0 , as in (3.5) and (3.6), one finds from The components F, G are always continuous at x=x 0 . If v=±s, so-called spin-and pseudospin symmetries, one of F, G also has a continuous derivative at x=x 0 . Note that the first-order Dirac differential equations do not see discontinuous potentials as a computational problem. The only conditions required are for F, G at x=0, and x = +¥. In cases of 'local' spin/pseudospin symmetries mentioned above one of the components F, G can be solved by a second-order Dirac equation. The problem occurs when using the second-order Dirac equations for general potentials.
In (3.8) and (3.9) one can substitute for G respectively F from (2.8) and (2.9). As a result (3.8) can be written in two These conditions need two separate solutions of F, (left/ right), being normalized in such a way that F is continuous at The conditions above need to be satisfied by bound states.

dF dx
= for constant potentials with l=0 (κ ¼ À1) The constant potentials are denoted v respectively s. Expressions for quantization conditions for the energy spectrum are then analytic. In particular for the F component with l=0 (κ=−1) and for the G component with l=1 (κ=+1). For general values of κ one has to employ special functions like spherical (Riccati) functions. For l=0, a regular solution F of (3.5) in the exterior region of the discontinuity is proportional to From relation (2.23) it follows that the G component is also defined for x>x 0 , since at x=x 0 : A regular component F in the well region for l=0, may be written and the two constants A and α needed to satisfy boundary conditions. One of the conditions, the vanishing wave at x=0, yields The continuity condition at x=x 0 yields where p is defined in (3.18). The condition (3.23) is used to obtain energies ò in tables 1 and 2.

Amplitude-phase approach
The separated second-order Dirac equations (3.5) and (3.6) require a different amplitude-phase approach than that in section 2. Fundamental (local) solutions of the second-order Dirac equations are and applies to both components F, G. Both amplitude and phase are real on the real x-axis. By insertion, a second-order equation results in a nonlinear (Milne-type) differential equation [15] w R x w w , 3 . 2 5 where R(x) may be continuous, or discontinuous as in (3.5) or (3.6). The amplitude functions w F and w G corresponding to separated second-order differential equations are different from u F and u G of the first-order differential equations. With respect to a discontinuous potential at x=x 0 the amplitudes w F and w G and their derivatives with respect to x are (always) continuous. If there is a single classically allowed interval of the x-axis, where the coefficient function satisfies R(x)>0, an amplitude function is defined at the maximum position x m of R(x). Here it is assumed that x m =x 0 for the discontinuous potentials and R(x 0 ) is the value to the left of the discontinuity. One finds The derivative of amplitude function could as well be defined as w 3.26). Any definition is admissible, but the one in (3.26) is the natural one in the general case of a smooth potential minimum. The integration of the amplitude equation is performed in two steps. The first step is to integrate towards the origin, sufficiently close. The second step consists in integrating across the discontinuity sufficiently far towards x  +¥. Two phases j L (x) and j R (x) are defined by integrating the amplitude function in two directions from x 0 : The phase x R j ( ) converges in the limit x  +¥ since w(x) increases indefinitely as x  +¥. If the effective potential has no discontinuity a bound state may be represented by a wave function defined by a single phase as The normalization constant C may be arbitrary in the present context. The quantization condition results from the assumed vanishing of the wave as r →+0 and as r  +¥. In terms of j L (x) and j R (x) one would have This condition (3.30) is similar to a Bohr-Sommerfeld-type bound-state quantization condition, where n is the number of nodes of the radial wave function in the well. However, n may be different for the Dirac components. The quantization condition (3.30) is used for the F component in tables 1 and 2. It ignores any special condition on the discontinuity of F¢ at x=x 0 . Since the coefficient R(x) contains quadratic terms of energy and potentials condition (3.30) is not a non-relativistic approximation. It corresponds to a bound-state condition for the Klein-Gordon equation, that reduces to the Schrödinger equation under non-relativistic conditions. If the effective potential is discontinuous, the wave function is no longer represented by (3.29). For x=x 0 the Dirac components are still continuous. To describe a correct wave for an arbitrary component, the normalization is chosen such that Ψ(x 0 )=1. A bound-state component wave can then be represented by in the well region, and in the classically forbidden region as x  +¥.
The relevant derivatives of (3.31) and (3.32) at x=x 0 expressed in terms of (3.27) and (3.28) with boundary conditions from (3.26) are With the aid of (3.10) this leads to the quantization condition for F: Note that the amplitudes and phases in (3.35) and (3.36) are unique to respective components, F or G. The quantization conditions are in principle equivalent, but present differing numerical challenges caused by the second term in each condition. If F is the large component, then condition (3.35) has a smaller second term and is numerically preferred. Figures 2 and 3 illustrate the component waves constructed from equations (3.31) and (3.32). Once one of the quantization conditions (3.35) and (3.36) is satisfied the relevant amplitude and phase functions are known and both Dirac components can be plotted. Further details are discussed in the subsequent section.

Numerical aspects and computational results
Several quantization conditions are derived in the preceding sections. They provide cross controls of the level energy calculations. The analytical quantization condition (3.23) for l=0 in the second-order differential equations for F is part of this control. The quantization conditions for the second-order differential equations are new, at least in context of amplitude-phase methods.
Computations are done in single precision on a laptop using a standard university environment. All quantization conditions agree to the number of decimals presented, but involve different computational difficulties. Both amplitudephase methods (first-/second-order Dirac equations) need some care when amplitude equations are integrated into classically forbidden regions. This relates to the existence of a classical forbidden region near the origin for either F or G. If integration is not truncated conveniently the corresponding amplitude function w(x) then tends to infinity as an irregular solution, while the phase function j(x) in (3.24) tends to zero towards the origin. Cancellation problems when using/plotting explicit waves occur in multiplications like for the regular representation w x x sin j ( ) ( ) in the limit x →+0. The integration of the amplitude equations from the discontinuity towards the origin is thus truncated at some x not too small, but where the phase function has converged.  (2.8) and (2.9). Note that F is positive and normalized to 1 at the discontinuity r x 0.210 02 = fm with x=26. (Right subplot) The small component G is obtained from (3.6) satisfying the condition (3.12). The large component F is obtained from G and G¢ of the equations (2.8) and (2.9). Here G is positive and normalized to 1 at the discontinuity r x 0.210 02 = fm with x=26.
Next, physical implications of the calculations are considered. Since radial square-well potential models are discontinuous approximations of Woods-Saxon potentials describing single neutron-nuclei interactions it is of interest to illustrate how a potential shape difference affects the energy spectrum. Shape differences are illustrated in figure 1. The figure shows how relativistic low-energy levels are shifted to lower values by the square-well shape. They also show that states near dissociation come close together and that levels due to the square-well shape lie closer to threshold than those of the Woods-Saxon shape. Figures 2 and 3 illustrate how spinor components appear when using the new quantization conditions (3.35) and (3.36) for second-order Dirac equations with radial square-well potentials. The condition (3.35) normalizes F to be unity at the discontinuity, while (3.36) normalizes G to be unity there. The two quantization conditions both provide separate numerical information sufficient to reconstruct the component waves, except possibly near the origin. Numerical cancellation errors may occur as discussed further below. It is obvious from figures 2 and 3 that the potential discontinuity affects the small component the most. The large component appears graphically as a smooth function. Figure 2 illustrates the same ground-state components, but obtained from (3.35) in the left subplot and (3.36) in the right subplot. The left subplot is obtained by firstly finding the lowest level energy ò with κ=−1. In this process all corresponding amplitude and phase functions are obtained and are used to reconstruct the F component (large red curve) and its derivative (not shown). The G component is finally reconstructed from F and F¢. The right subplot in figure 2 starts with solving (3.36) for the energy and follows the same logic of constructing the components. One notices that the small component is positive because of the normalization at the discontinuity.
While both components in figure 2 have the same number of nodes, their node numbers differ in figure 3. The large component of the ground state for κ=1 has no node, but the small component has one node. For κ=1 there is no centrifugal barrier in the second-order differential equation for G. Obviously, there is no G component for κ=1 with zero node number. Tables 1 and 2 show computed binding energies, ò−1, obtained by different quantization conditions for discontinuous and smooth potential models. Results for the F component when ignoring the discontinuity condition, ò c −1, are also shown for comparison. The results ò a of tables 1 and 2 obtained by several alternative quantization conditions for discontinuous potentials are in numerical agreement. For l=0 (κ=−1) the analytic condition (3.23) supports the results.
Since in the non-relativistic limit the effective Schrödinger potential is the sum of the vector and scalar potentials, tables 1 and 2 show two possible combinations of vector and scalar potentials with the same sum of them. Single neutronnuclei spectra shown in text books and recent articles indicate that scalar-type potential wells are important for the order of states in the spin-orbit level splittings.
In table 1 the vector potential is zero. A discontinuous attractive scalar potential and a similar Woods-Saxon potential is used to compare binding energies relevant for single neutron-nuclei shell states. The scalar Woods-Saxon potential is  table 1 shows that the binding energies of the discontinuous potential are more negative than those for the Woods-Saxon model if energies are far below threshold. In some neighborhood of threshold it seems that this ʼshape-ordering' of states is reversed; see entry κ=−1, n=2. This state corresponds to the term s 3 1 2 , where the notation here is n l 1 j + ( ) with l the quantum number of angular momentum and j the total angular momentum. Spinorbit ordering of states ( s  p  p  1 , 1 1 2 ) with increasing energies for a given potential is the same. If an equally attractive vector potential is used, the energy order changes to s p p 1 , 1 3 2  for both discontinuous and smooth  (2.8) and (2.9). Note that F is positive and normalized to 1 at the discontinuity r x 0.210 02 = fm with x=26. (Right subplot) The small component G is obtained from (3.6) satisfying the condition (3.12). The large component F is obtained from G and G¢ of the equations (2.8) and (2.9). Here G is positive and normalized to 1 at the discontinuity r=0.210 02x fm with x=26. models. The splittings of spin-orbit pairs for both potential shapes are similar close to threshold. Table 2 shows results for mixed scalar/vector potentials defined below, one is discontinuous and one is continuous. As reviewed by [5], referring to field-theory results, a combination of vector/scalar potentials may be relevant for single nucleon-nuclei spectra. It is suggested that the attraction in this combined vector/scalar case is due to the scalar potential, the vector potential being repulsive. In this spirit Guo et al [6] present graphical results for a given Woods-Saxon potential modeling shell states of a neutron-Pb nucleus. The present Woods-Saxon potential has a similar ratio of thickness/ range, but the range is somewhat shorter. The Woods-Saxon potential model is specified by Energy states in table 2 are generally shifted upwards compared to those of the scalar potential in table 1, but the order of states is the same. Some more states are closer to threshold and one observes the ʼshape-order' shift (as in table 1) ). This is likely due to the wider Woods-Saxon well near threshold, when compared with the sharp discontinuous well. The splittings of spin-orbit pairs for both potential shapes are similar close to threshold.
In both tables the Klein-Gordon type binding energies obtained from (3.30) are also shown. Their computations are based on the second-order F-equation for positive energies and ignores discontinuities of F¢. Energies obtained by (3.30) fail to explain spin-orbit splittings, but seem to provide close upper bounds to the discontinuous model not too close to threshold.

Conclusions
The present study suggests alternative numerical methods relevant for Dirac equations with sharp-step/discontinuous potentials. As a smooth potential well approaches a sharp-step potential, energy levels are generally shifted to lower values. However, levels near dissociation are shifted in the opposite direction.
A general amplitude-phase method applicable to firstorder Dirac differential equations is already known [1], but is also presented here. For sharp potential shapes an approximate, discontinuous model may be of interest. An amplitudephase method for second-order Dirac differential equations to handle such potentials is derived.
Computations suggest that a small spinor component at positive energies is more sensitive than the large component to the discontinuity of the potential model. For almost nonrelativistic conditions the second-order Dirac equations for the large component can thus be approximated by ignoring the discontinuity condition. Resulting energy levels are close upper bounds to the relevant spin-orbit level pairs for the discontinuous potential.
To summarize, the present study introduces a discontinuity in a potential model to approximate a smooth step in the 'original/actual' potential model of second-order Dirac equations. The conclusion is that: (a) a non-relativistic Schrödinger problem with a sharp (Woods-Saxon type) potential well can be accurately approximated by a squarewell potential and standard boundary conditions. (b) However, this is not always true for a corresponding Schrödinger equation derived from a Dirac problem without modifying quantization conditions.