Lagrange-Mesh Method for Deformed Nuclei With Relativistic Energy Density Functionals

The application of relativistic energy density functionals to the description of nuclei leads to the problem of solving self-consistently a coupled set of equations of motion to determine the nucleon wave functions and meson fields. In this work, the Lagrange-mesh method in spherical coordinates is proposed for numerical calculations. The essential field equations are derived from the relativistic energy density functional and the basic principles of the Lagrange-mesh method are delineated for this particular application. The numerical accuracy is studied for the case of a deformed relativistic harmonic oscillator potential with axial symmetry. Then the method is applied to determine the point matter distributions and deformation parameters of self-conjugate even-even nuclei from ⁴He to ⁴⁰Ca.

1. Introduction

The theoretical description of nuclei with the help of energy density functionals (EDFs) has advanced steadily during the last decades. Many improvements were introduced in the physical modeling and technical methods were developed continuously. EDFs can be designed using heuristic arguments but their precise form in nuclear physics is often deduced from and guided by systematic approaches employing mean-field approximations. Two major classes can be distinguished: non-relativistic EDFs derived from effective interactions, e.g., Skyrme or Gogny type, in the Hartree–Fock approximation and relativistic or covariant models that start from a field-theoretical formalism with a Lagrangian density. Refer to the review [1] for details and examples of both classes.

The application of EDFs to describe deformed nuclei represents a particular challenge because of the larger complexity of the calculations that is needed to find self-consistent solutions of the coupled equations. This is particular true for relativistic models in which the Dirac equation has to be solved to obtain the nucleon single-particle states and Klein–Gordon or Poisson equation for the mesons and electromagnetic fields [2–9]. In contrast, the Schrödinger equation with given potentials is relevant for non-relativistic EDFs. Different methods have been applied to solve these partial differential equations. They can be discretized on a mesh with suitably distributed grid points, usually in coordinate space using shooting and matching techniques or finite element methods to find solutions. See, e.g., Horowitz and Serot [10], Bonche et al. [11], Reinhard et al. [12], Pöschl et al. [13, 14], Meng [15], Stoitsov et al. [16], and Typel and Wolter [17]. Alternatively, in configuration space, solutions can be expanded in a set of carefully chosen basis functions. Examples are eigenfunctions of the harmonic oscillator, Woods–Saxon potential, or transformed harmonic oscillator functions. For examples, refer to Gambhir and Ring [18], Ring et al. [19], Zhou et al. [20], and Nikšić et al. [21, 22]. Then the diagonalization of more or less large matrices is required. Also B-splines or wavelet representations have been used; see, e.g., Umar et al. [23], Pei et al. [24, 25]. Instead of diagonalizing large matrices, imaginary time-step methods have been exploited with great success to obtain wave functions from non-relativistic EDFs; see, e.g., Davies et al. [26], Reinhard and Cusson [27], and Ryssens et al. [28]. Another method is based on Fourier transformation techniques [29]. A comparative study of the numerical accuracy of different methods in non-relativistic approaches can be found in Bonche et al. [11]. The dependence of the precision on the characteristic parameters of particular numerical methods were also explored [30, 31]. Problems of the imaginary time-step method in solving the Dirac equation on a three-dimensional lattice and the accuracy of the method are discussed in Tanimura et al. [32]. An efficient and precise approach is the Lagrange-mesh method that combines the virtues of the discretization and basis methods previously mentioned. This variational approach has been developed mainly to study non-relativistic problems and found many applications [33–37]. There are only few applications of the Lagrange-mesh method to relativistic systems, mainly for simple cases like hydrogenic atoms [38, 39].

In this work, the Lagrange-mesh method is applied to find solutions of the Dirac equation with non-spherical potentials that appear in the description of deformed but axially symmetric nuclei. The main difference is that the kinetic energy operator in the non-relativistic case is given by a second derivative, whereas only first derivatives occur in the relativistic case. As an example, a relativistic EDF that is derived from a relativistic mean-field model with density-dependent meson-nucleon couplings is considered. The method is easily adapted to other types of relativistic EDFs, e.g., with nonlinear self-interactions of the mesons.

This work starts in section 2 with a short introduction to the relativistic EDF with density-dependent couplings that is used to describe the structure of nuclei. The set of coupled equations that have to be solved self-consistently is derived. In section 3, the numerical techniques are introduced. After summarizing the basic principles of the Lagrange-mesh method, the discussion focuses on how solutions of the Dirac equation and the meson field equations are obtained. Results are presented in section 4. First, deformed relativistic harmonic oscillator potentials are considered. They are used to study the precision of the numerical methods because quasi-analytical results for the energy levels are available. Then the techniques are applied to the description of self-conjugated even-even nuclei from ⁴He to ⁴⁰Ca as they show a large variation of different shapes. Finally, conclusions are given in section 5. In this work, natural units with ℏ = c = 1 are used in all the equations. For the conversion of units, the value ℏc = 197.3269788 MeV fm is utilized [40].

2. Relativistic Energy Density Functional

The general starting point for relativistic models of nuclei and nuclear matter is a covariant Lagrangian density ${L}$. It contains the relevant terms to describe nucleons, mesons, and their interaction, in most cases realized by minimal couplings. From ${L}$, the equations of motion of all degrees of freedom are deduced with the help of the Euler-Lagrange equations and by applying approximations in the subsequent steps. Usually, the mean-field and no-sea approximation are used and meson fields are treated as classical fields. The system of coupled equations is simplified further by presuming certain symmetries. Finally, quantities like the energy density can be obtained from the energy-momentum tensor. This procedure to derive an EDF from the Lagrangian density and its application to nuclear structure calculation has been presented extensively in the literature; see, e.g., Bender et al. [1], Serot and Walecka [2], Ring [3], Serot and Walecka [4], Lalazissis et al. [5], Vretenar et al. [6], Meng et al. [7], Nikšić et al. [8], and Ring et al. [9]. Here, it is sufficient to give only the explicit form of the relativistic EDF itself and the equations of motion.

In the present work, a relativistic EDF is employed assuming stationary systems in a fixed frame of reference. In general, it considers four types of mesons to model the effective in-medium interaction: a scalar σ meson to describe the attraction between nucleons and a vector ω meson for repulsion. Besides these isoscalar mesons, a vector ρ meson and a scalar δ meson are introduced for the isospin dependence of the strong force. Furthermore, the electromagnetic interaction is included. These fields are denoted in the following equations with the symbols σ, ω₀, ρ₀, δ, and A₀, respectively. Because of the assumptions, only the zero component of the Lorentz vector fields are nonzero and there is only one nonzero component of vectors in isospin space. With the spinors Ψ_ik for a proton (i = p) or a neutron (i = n) in a single-particle state k, it is convenient to introduce the kinetic energy densities of the nucleons

$\begin{array}{cc}\begin{array}{l}{\tau }_i={\displaystyle \sum _k}{w}_{ik}{\overline{\text{Ψ}}}_{ik}\overrightarrow{\alpha }·\overrightarrow{p}{\text{Ψ}}_{ik},\end{array}& (1)\end{array}$

the vector densities

$\begin{array}{cc}\begin{array}{l}{n}_i^{(v)}={\displaystyle \sum _k}{w}_{ik}{n}_{ik}^{(v)}={\displaystyle \sum _k}{w}_{ik}{\overline{\text{Ψ}}}_{ik}{\gamma }^0{\text{Ψ}}_{ik},\end{array}& (2)\end{array}$

and the scalar densities

$\begin{array}{cc}\begin{array}{l}{n}_i^{(s)}={\displaystyle \sum _k}{w}_{ik}{n}_{ik}^{(s)}={\displaystyle \sum _k}{w}_{ik}{\overline{\text{Ψ}}}_{ik}{\text{Ψ}}_{ik}\end{array}& (3)\end{array}$

with ${\overline{\text{Ψ}}}_{ik}={\text{Ψ}}_{ik}^{†}{\gamma }^0$ in the usual notation for the relativistic matrices $\overrightarrow{\alpha }$ and γ⁰ = β and the momentum operator $\overrightarrow{p}$. The symbol w_ik is the occupation factor of the state k. In addition, the source densities

$\begin{array}{cc}\begin{array}{l}{n}_m={\displaystyle \sum _{i=n,p}}{g}_{im}{n}_i^{(s)}\end{array}& (4)\end{array}$

for m = σ, δ and

$\begin{array}{cc}\begin{array}{l}{n}_m={\displaystyle \sum _{i=n,p}}{g}_{im}{n}_i^{(v)}\end{array}& (5)\end{array}$

for m = ω, ρ, γ are defined with factors g_nσ = g_pσ = g_nω = g_pω = −g_nρ = g_pρ = −g_nδ = g_pδ = g_pγ = 1 and g_nγ = 0 that specify the couplings between nucleons, mesons, and photons. Then the energy density ε has the explicit form

$\begin{array}{cc}\begin{array}{l}\epsilon \text{}=\text{}{\displaystyle \sum _{i=n,p}\left({\tau }_i\text{}+\text{}{m}_i{n}_i^{(s)}\right)}-{\Gamma }_{\sigma }\sigma n_{\sigma }\text{}-\text{}{\Gamma }_{\delta }\delta n_{\delta }\text{}+\text{}{\Gamma }_{\omega }{\omega }_0{n}_{\omega }\text{}+\text{}{\Gamma }_{\rho }{\rho }_0{n}_{\rho }\\ \text{ }+{\Gamma }_{\gamma }{A}_0{n}_{\gamma }+\frac{1}{2}(\overrightarrow{\nabla }\sigma ·\overrightarrow{\nabla }\sigma +m_{\sigma }^2{\sigma }^2+\overrightarrow{\nabla }\delta ·\overrightarrow{\nabla }\delta +m_{\delta }^2{\delta }^2\\ \text{ }-\text{}\overrightarrow{\nabla }{\omega }_0·\overrightarrow{\nabla }{\omega }_0\text{}-\text{}{m}_{\omega }^2{\omega }_0^2\text{}-\text{}\overrightarrow{\nabla }{\rho }_0·\overrightarrow{\nabla }{\rho }_0-m_{\rho }^2{\rho }_0^2\text{}-\text{}\overrightarrow{\nabla }{A}_0·\overrightarrow{\nabla }{A}_0)\end{array}& (6)\end{array}$

with masses of the nucleons m_i (i = n, p) and those of the mesons m_m (m = σ, ω, ρ, δ). It is a functional of the various nucleon densities, the meson and electromagnetic fields, and their derivatives. The meson-nucleon couplings Γ_σ, Γ_ω, Γ_ρ, and Γ_δ, are assumed to depend on the total vector density $n_{\text{tot}}=\sum _{i=n,p}{n}_i^{(v)}$ of the nucleons in the present work. They determine the effective in-medium interaction. ${\Gamma }_{\gamma }=\sqrt{4\pi \alpha }$ with α = e²/(ℏc) = 1/137.035999139 [40] is the standard electromagnetic coupling constant.

The EDF (6) can be used to derive the field equations of all degrees of freedom using the Euler-Lagrange equations. The time-independent, inhomogeneous, second-order differential equations

$\begin{array}{cc}\begin{array}{lll}-\Delta \sigma +m_{\sigma }^2\sigma & =& {\Gamma }_{\sigma }{n}_{\sigma }\end{array}& (7)\end{array}$

$\begin{array}{cc}\begin{array}{lll}-\Delta \delta +m_{\delta }^2\delta & =& {\Gamma }_{\delta }{n}_{\delta }\end{array}& (8)\end{array}$

$\begin{array}{cc}\begin{array}{lll}-\Delta {\omega }_0+m_{\omega }^2{\omega }_0& =& {\Gamma }_{\omega }{n}_{\omega }\end{array}& (9)\end{array}$

$\begin{array}{cc}\begin{array}{lll}-\Delta {\rho }_0+m_{\rho }^2{\rho }_0& =& {\Gamma }_{\rho }{n}_{\rho }\end{array}& (10)\end{array}$

for the mesons and

$\begin{array}{cc}\begin{array}{l}-\Delta A_0={\Gamma }_{\gamma }{n}_{\gamma }\end{array}& (11)\end{array}$

for the electromagnetic field are found. The single-particle states of the nucleons are obtained from solving the Dirac equation

$\begin{array}{cc}\begin{array}{c}H{\Psi }_{ik}(\overrightarrow{r})=[\overrightarrow{\alpha }·\overrightarrow{p}+\beta (m_i-S_i)+V_i]{\Psi }_{ik}(\overrightarrow{r})=E_{ik}{\Psi }_{ik}(\overrightarrow{r})\end{array}& (12)\end{array}$

with the energies E_ik (including the particle rest mass) and a Hamilton operator H that explicitly contains scalar and vector potentials. They are given by

$\begin{array}{cc}\begin{array}{l}{S}_i=g_{i\sigma }{\Gamma }_{\sigma }\sigma +g_{i\delta }{\Gamma }_{\delta }\delta \end{array}& (13)\end{array}$

and

$\begin{array}{cc}\begin{array}{l}{V}_i=g_{i\omega }{\Gamma }_{\omega }{\omega }_0+g_{i\rho }{\Gamma }_{\rho }{\rho }_0+g_{i\gamma }{\Gamma }_{\gamma }{A}_0+V^{(R)}\end{array}& (14)\end{array}$

with a rearrangement contribution

$\begin{array}{cc}\begin{array}{l}{V}^{(R)}=\frac{d{\Gamma }_{\omega }}{d{n}_{\text{tot}}}{n}_{\omega }{\omega }_0+\frac{d{\Gamma }_{\rho }}{d{n}_{\text{tot}}}{n}_{\rho }{\rho }_0-\frac{d{\Gamma }_{\sigma }}{d{n}_{\text{tot}}}{n}_{\sigma }\sigma -\frac{d{\Gamma }_{\delta }}{d{n}_{\text{tot}}}{n}_{\delta }\delta \end{array}& (15)\end{array}$

that appears in the vector potential V_i because of the dependence of the couplings on the total vector density n_tot. For other choices of the density dependence, see Typel [41]. The set of coupled Equations (7)–(11) and (12) has to be solved self-consistently with the appropriate boundary conditions for a given number of neutrons (N) and protons (Z). Finally, the energy of a nucleus with mass number A = N+Z is obtained by integrating the energy density (6) over all space. It can be written as

$\begin{array}{cc}\begin{array}{l}E(N,Z)={\displaystyle \sum _{i=n,p}}{\displaystyle \sum _k}{w}_{ik}{E}_{ik}+\int d^{3}r{\epsilon }_{\text{field}}(\overrightarrow{r})\end{array}& (16)\end{array}$

with occupation factors w_ik of the single-particle states. The contribution of the fields

$\begin{array}{cc}\begin{array}{l}\begin{array}{ccc}{\epsilon }_{field}(\overrightarrow{r})& =& \frac{1}{2}({\Gamma }_{\sigma }{n}_{\sigma }\sigma +{\Gamma }_{\delta }{n}_{\delta }\delta -{\Gamma }_{\omega }{n}_{\omega }{\omega }_0-{\Gamma }_{\rho }{n}_{\rho }{\rho }_0-{\Gamma }_{\gamma }{n}_{\gamma }{A}_0)\end{array}\\ \begin{array}{c}-V^{(R)}{n}_{tot}\end{array}\end{array}& (17)\end{array}$

contains an explicit rearrangement term and is found from (6) by using partial integration and the field equations. For a comparison with experimental data, corrections have to be added to the energy (16), e.g., for the breaking of translational and rotational symmetries in the numerical calculation. Also the effect of pairing correlations has to be considered. However, these effects are beyond the scope of the present work since the application of the Lagrange-mesh method to the field equations is the focus here.

3. Numerical Techniques

Solutions of the field Equations (7)–(12) can be obtained with different numerical techniques. These are most efficient when they are adapted to the particular system of interest, in particular, if the available symmetries are exploited. For the problem at hand, deformed nuclei with axial symmetry, an expansion of the single-particle wave functions and meson fields using a basis of eigenfunctions of a deformed harmonic oscillator seems to be the natural choice. But these functions drop off much faster with increasing distance than expected for the wave functions of nucleons bound in a nucleus because these basis functions have a functional behavior of polynomial times Gaussian form. This fact necessitates the inclusion of functions up to very high numbers of oscillator quanta leading to a large basis set that has to be used in the diagonalization to find the eigenfunctions. This problem is partly solved by switching to suitably transformed harmonic oscillator functions or by choosing basis functions obtained as solutions of the Dirac equation with Woods–Saxon potentials. An alternative is the discretization of the wave functions and fields on a mesh of well-chosen grid points using cylindrical coordinates. Function values at points between the mesh points are found with suitable interpolation rules. To capture the correct asymptotic behavior, the grid has to be extended to radii far outside the nucleus itself. However, cylindrical coordinates are not well adapted to describe the transition from a strongly deformed system to a nucleus with spherical symmetry.

To avoid some of the problems previously mentioned, the Lagrange-mesh method with a representation of wave functions in spherical coordinates is used in this work. The main ideas of this approach are presented in section 3.1. An extended discussion of the Lagrange-mesh method can be found in the reviews [33, 34, 37]. The specific application to find the single-particle wave functions and fields for axially-symmetric deformed nuclei is given in the following two subsections.

3.1. Lagrange-Mesh Method

Solutions of a differential equation in a variable x can be approximated by a linear combination of mesh functions that are carefully chosen in the Lagrange-mesh method to exhibit some specific properties. In the first step, a set of polynomials p_k(x) with k = 0, 1, 2, …, N on an interval [a, b] is selected that are orthogonal with respect to a certain non-negative weight function w(x). Examples are Legendre, Laguerre, or Hermite polynomials that are defined on intervals [−1, 1], [0, ∞], and [−∞, ∞] with weight functions w = 1, w = exp(−x), and w = exp(−x²), respectively [42]. The polynomials are used to define basis functions

$\begin{array}{cc}\begin{array}{l}{\varphi }_k(x)=\sqrt{\frac{w_k(x)}{h_k}}{p}_k(x)\end{array}& (18)\end{array}$

with positive normalization factors h_k so that the functions (18) form an orthonormal basis, i.e.,

$\begin{array}{cc}\begin{array}{l}{\displaystyle {\int }_a^{b}dx {\varphi }_i(x){\varphi }_j(x)={\delta }_{ij}}\hfill \end{array}& (19)\end{array}$

on the interval [a, b]. From this set of basis functions, the Lagrange-mesh functions

$\begin{array}{cc}\begin{array}{l}{f}_i(x)={\lambda }_i{\displaystyle \sum _{k=0}^{N-1}}{\varphi }_k(x){\varphi }_k(x_i)\end{array}& (20)\end{array}$

with i = 1, 2, …, N are constructed with normalization factors

$\begin{array}{cc}\begin{array}{l}{\lambda }_i={\left({\displaystyle \sum _{k=0}^{N-1}}|{\varphi }_k(x_i){|}^2\right)}^{-1}\end{array}& (21)\end{array}$

where the basis functions ϕ_k are evaluated at the N distinct simple zeroes x_i of the function ϕ_N(x). The definition (20) with the normalization (21) guarantees that the condition

$\begin{array}{cc}\begin{array}{l}{f}_i(x_j)={\delta }_{ij}\end{array}& (22)\end{array}$

is fulfilled. The mesh functions (20) do not constitute an orthonormal set because

$\begin{array}{cc}\begin{array}{l}{\displaystyle {\int }_b^{a}dx f_i(x)f_j(x)={\lambda }_i{\delta }_{ij};}\hfill \end{array}& (23)\end{array}$

instead, the mesh functions constitute an orthogonal basis. With help of the Christoffel–Darboux formula [42], the mesh functions can also be represented as

$\begin{array}{cc}\begin{array}{l}{f}_i(x)=\frac{{\varphi }_N(x)}{(x-x_i){\varphi }_N^{\prime }(x_i)}\end{array}& (24)\end{array}$

for x≠x_i and the normalization factors can be expressed as

$\begin{array}{cc}\begin{array}{l}{\lambda }_i=\frac{1}{x_i{\left[{\varphi }_N^{\prime }(x_i)\right]}^2}\end{array}& (25)\end{array}$

depending only on the basis function ϕ_N and its derivative ${\varphi }_N^{\prime }$. The first derivative of the Lagrange-mesh function (24) is given by

$\begin{array}{cc}\begin{array}{l}{f}_i^{\prime }(x)=\frac{(x-x_i){\varphi }_N^{\prime }(x)-{\varphi }_N(x)}{{(x-x_i)}^2{\varphi }_N^{\prime }(x_i)}\end{array}& (26)\end{array}$

with the particular values

$\begin{array}{cc}\begin{array}{l}{f}_i^{\prime }(x_j)={\lambda }_i{x}_i\frac{{\varphi }_N^{\prime }\left(x_i\right){\varphi }_N^{\prime }(x_j)}{(x_j-x_i)}={(\frac{{\lambda }_i{x}_i}{{\lambda }_j{x}_j})}^{\frac{1}{2}}\frac{1}{(x_j-x_i)}\frac{{\varphi }_N^{\prime }(x_i)}{\left|{\varphi }_N^{\prime }(x_i)\right|}\frac{{\varphi }_N^{\prime }(x_j)}{|{\varphi }_N^{\prime }(x_j)|}\end{array}& (27)\end{array}$

for x_j ≠ x_i and

$\begin{array}{cc}\begin{array}{l}{f}_i^{\prime }(x_i)=\frac{{\varphi }_N^{\prime \prime }(x_i)}{2{\varphi }_N^{\prime }(x_i)}\end{array}& (28)\end{array}$

when the indices of the function and of the zero coincide. These results will be useful in the later application.

The Lagrange-mesh functions f_i are used as a basis to represent any function g(x) on the interval [a, b] by an N-component vector ${(g_1,g_2,\dots ,g_N)}^T$ with coefficients g_i = g(x_i) that are simply the function values at the zeroes x_i of ϕ_N. Because of this condition (22), the expansion

$\begin{array}{cc}\begin{array}{l}g(x)={\displaystyle \sum _{i=1}^N}{g}_i{f}_i(x)\end{array}& (29)\end{array}$

is exact at all zeroes x_i. For x ≠ x_i, values of g(x) are easily determined as usual for a basis expansion without the need for an additional interpolation method that is required for other representations on a grid. If a function g(x) can be written as g(x) = P(x)w(x) with a polynomial P(x) of degree not larger than 2N−1, then the integral

$\begin{array}{cc}\begin{array}{l}{\displaystyle {\int }_a^{b}dx g(x)=\underset{i = 1}{\overset{N}{{{\displaystyle \sum }}^{\text{}}}}{\lambda }_{i}g(x_i)}\hfill \end{array}& (30)\end{array}$

is exactly given by a simple sum. This formula can be considered as a generalized Gaussian quadrature rule.

In the application of this work, the Lagrange-mesh method is used to represent functions of the radius variable r in spherical coordinates. Then Laguerre polynomials $L_k^{(\alpha )}(x)$, which are defined for 0 ≤ x < ∞, are the appropriate choice as the fundamental orthogonal polynomials p_k in (18) with the weight function w^(α) = x^αexp(−x) and the normalization factor $h_k^{(\alpha )}=\Gamma (\alpha +k+1)/k!$. The basis functions (18) in the Lagrange-mesh functions (20) depend on the parameter α, i.e.,

$\begin{array}{cc}\begin{array}{l}{\varphi }_k^{(\alpha )}(x)={\left[h_k^{(\alpha )}\right]}^{-\frac{1}{2}}{L}_k^{(\alpha )}(x)x^{\frac{\alpha }{2}}exp(-\frac{x}{2})\end{array}& (31)\end{array}$

and decrease exponentially for large arguments x. The choice of α will be detailed below. First and second derivatives are given by

$\begin{array}{cc}\begin{array}{l}\begin{array}{ccc}\frac{d}{dx}{\varphi }_k^{(\alpha )}(x)& =& {[h_k^{(\alpha )}]}^{-\frac{1}{2}}[(\frac{\alpha }{2x}-\frac{1}{2})L_k^{(\alpha )}(x)+L_k^{(\alpha )\prime }(x\left)\right]x^{\frac{\alpha }{2}}\end{array}\\ \begin{array}{c}\exp (-\frac{x}{2})\end{array}\end{array}& (32)\end{array}$

and

$\begin{array}{cc}\begin{array}{l}\frac{d^2}{d{x}^2}{\varphi }_k^{(\alpha )}(x)={\left[h_k^{(\alpha )}\right]}^{-\frac{1}{2}}[\left(\frac{\alpha (\alpha -2)}{4x^2}-\frac{\alpha +2k}{2x}+\frac{1}{4}\right)L_k^{(\alpha )}(x)\\ -\frac{1}{x}{L}_k^{(\alpha )\prime }(x)]x^{\frac{\alpha }{2}}\exp \left(-\frac{x}{2}\right),\end{array}& (33)\end{array}$

respectively, where the second-order differential equation for Laguerre polynomials

$\begin{array}{cc}\begin{array}{l}x{L}_k^{(\alpha )\prime \prime }(x)+(\alpha +1-x)L_k^{(\alpha )\prime }(x)+k{L}_k^{(\alpha )}(x)=0\end{array}& (34)\end{array}$

has been used in the derivation of Equation (33).

Any function ψ(x) is represented by its function values ${\psi }_i^{(\alpha )}=\psi (x_i^{(\alpha )})$ or ${\psi }_i^{(\beta )}=\psi (x_i^{(\beta )})$ at the N zeroes $x_i^{(\alpha )}$ or $x_i^{(\beta )}$ of the Laguerre polynomial $L_N^{(\alpha )}(x)$ or $L_N^{(\beta )}(x)$ as

$\begin{array}{cc}\begin{array}{l}\psi (x)={\displaystyle \sum _{i=1}^N{\psi }_i^{(\alpha )}{f}_i^{(\alpha )}(x)\mathrm{\text{or\hspace{1em} }}}\psi (x)={\displaystyle \sum _{i=1}^N}{\psi }_i^{(\beta )}{f}_i^{(\beta )}(x)\end{array}& (35)\end{array}$

with Lagrange-mesh functions $f_i^{(\alpha )}(x)$ or $f_i^{(\beta )}(x)$. The transformation between the two vectors is given by

$\begin{array}{cc}\begin{array}{l}{\psi }_j^{(\beta )}={\displaystyle \sum _{i=1}^N}{U}_{ji}^{(\beta \alpha )}{\psi }_i^{(\alpha )}\end{array}& (36)\end{array}$

with a matrix U^(βα) containing the real matrix elements elements $U_{ji}^{(\beta \alpha )}=f_i^{(\alpha )}(x_j^{(\beta )})$. Integrals of the form

$\begin{array}{cc}\left(\psi \left|\mathcal{O}\right|\chi \right)={\displaystyle {\int }_0^{\infty }d}x{\psi }^{\ast }(x)\mathcal{O}\chi (x)& (37)\end{array}$

with an operator ${O}$ acting on the coordinate x can be written as

$\begin{array}{cc}\begin{array}{l}(\psi \left|{O}\right|\chi )={\displaystyle \sum _{i,j=1}^N}{\left[{\lambda }_i^{(\beta )}{\lambda }_j^{(\alpha )}\right]}^{\frac{1}{2}}{\psi }_i^{(\beta )*}{{O}}_{ij}^{(\beta \alpha )}{\chi }_j^{(\alpha )}\end{array}& (38)\end{array}$

with operator matrix elements

$\begin{array}{cc}{\mathcal{O}}_{ij}^{(\beta \alpha )}={\left[{\lambda }_i^{(\beta )}{\lambda }_j^{(\alpha )}\right]}^{-\frac{1}{2}}{\displaystyle {\int }_0^{\infty }d}x{f}_i^{(\beta )}(x)\mathcal{O}{f}_j^{(\alpha )}(x)& (39)\end{array}$

containing the real Lagrange-mesh functions. They transform as

$\begin{array}{cc}\begin{array}{l}{{O}}_{ij}^{(\beta \alpha )}={\left[{\lambda }_i^{(\beta )}\right]}^{-\frac{1}{2}}{\displaystyle \sum _{k=1}^N}{\left[{\lambda }_k^{(\gamma )}\right]}^{\frac{1}{2}}{U}_{ki}^{(\gamma \beta )}{{O}}_{kj}^{(\gamma \alpha )}\end{array}& (40)\end{array}$

when changing between different sets of basis functions. The prefactors are chosen in (39) such that the unit operator ${O}=1$ is diagonal in the (α, α) representation with ${(1)}_{ij}^{(\alpha \alpha )}={\delta }_{ij}$. Of particular interest are matrix elements of the derivative operator ${O}=\frac{d}{dx}$. They can be expressed as

$\begin{array}{cc}{\left(\frac{d}{dx}\right)}_{ij}^{(\alpha \alpha )}={\left[{\lambda }_i^{(\alpha )}{\lambda }_j^{(\alpha )}\right]}^{-\frac{1}{2}}{\displaystyle {\int }_0^{\infty }d}x{f}_i^{(\alpha )}(x)f_j^{(\alpha )\prime }(x)& (41)\end{array}$

with the integrand having the form of the function g(x) in (30). Thus they are exact in the chosen basis. These matrix elements will appear when the Dirac equation has to be solved; see section 3.2. Using the orthogonality condition (22), the simple expression

$\begin{array}{cc}\begin{array}{l}\begin{array}{ccc}{\left(\frac{d}{dx}\right)}_{ij}^{(\alpha \alpha )}& =& {[{\lambda }_i^{(\alpha )}{\lambda }_j^{(\alpha )}]}^{-\frac{1}{2}}\underset{N}{\overset{k = 1}{{{\displaystyle \sum }}^{\text{}}}}{\lambda }_k^{(\alpha )}{f}_i^{(\alpha )}(x_k^{(\alpha )})f_j^{(\alpha )\prime }(x_k^{(\alpha )})\end{array}\\ \begin{array}{c}={[\frac{{\lambda }_i^{(\alpha )}}{{\lambda }_j^{(\alpha )}}]}^{\frac{1}{2}}{f}_j^{(\alpha )\prime }(x_i^{(\alpha )})\end{array}\end{array}& (42)\end{array}$

is obtained with the specific values

$\begin{array}{cc}\begin{array}{l}{(\frac{d}{dx})}_{ij}^{(\alpha \alpha )}=\frac{{(-1)}^{i-j}}{\left[x_i^{(\alpha )}-x_j^{(\alpha )}\right]}{\left[\frac{x_j^{(\alpha )}}{x_i^{(\alpha )}}\right]}^{\frac{1}{2}}\end{array}& (43)\end{array}$

for i≠j and

$\begin{array}{cc}\begin{array}{l}{\left(\frac{d}{dx}\right)}_{ii}^{(\alpha \alpha )}=-\frac{1}{2x_i^{(\alpha )}}\end{array}& (44)\end{array}$

for equal indices when Equations (25), (27), (28), (32), and (33) are taken into account. The factor (−1)^i−j arises from the alternating sign of the derivatives of the basis functions ${\varphi }_N^{(\alpha )}$ at consecutive zeroes $x_i^{(\alpha )}$.

The practical application of the Lagrange-mesh method for the radial coordinate r is straightforward. A Laguerre polynomial $L_N^{(\alpha )}(x)$ of sufficiently high degree N is chosen and its zeroes $x_i^{(\alpha )}$ for i = 1, …, N are determined. Then the normalization factors ${\lambda }_i^{(\alpha )}$ are given by Equation (25). Every function g(r) is represented by its values $g_i^{(\alpha )}=g(r_i^{(\alpha )})$ at the mesh points $r_i=h{x}_i^{(\alpha )}$ containing a scale factor h. It determines the spatial resolution and has to be considered in the calculation of integrals over the radial coordinate r, e.g., as a prefactor appearing in

$\begin{array}{cc}{\displaystyle {\int }_0^{\infty }d}rg(r)=h{\displaystyle {\int }_0^{\infty }d}xg(hx)=h{\displaystyle \sum _{i=1}^N{\lambda }_i^{(\alpha )}}{g}_i^{(\alpha )}.& (45)\end{array}$

The factor h has to be fixed to an appropriate value depending on the physical problem.

3.2. Solving the Dirac Equation

For a nucleon with mass m moving in vector and scalar fields V and S (suppressing the index i to distinguish between neutrons and protons in the following equations), the single-particle Hamiltonian in the Dirac Equation (12) has the form

$\begin{array}{cc}\begin{array}{l}H=\left(\begin{array}{cc}V+m-S& \overrightarrow{\sigma }·\overrightarrow{p}\\ \overrightarrow{\sigma }·\overrightarrow{p}& V-m+S\end{array}\right)\end{array}& (46)\end{array}$

with a distinct 2 × 2 structure. Correspondingly, it is useful to represent solutions ${\text{Ψ}}_k(\overrightarrow{r})$ of (12) as spinors with an upper and a lower component. In the case of spherical symmetry, when V and S depend only on the radius r, individual states k can be identified uniquely with three quantum numbers k = (n, κ, Ω). The integer κ = ±1, ±2, ±3, …  determines the total angular momentum j = | κ | −1/2 of the single-particle state, the orbital angular momentum l = | κ | −1/2−κ/(2|κ|), and the parity π = (−1)^l of the upper component of Ψ_k. Vice versa, one has κ = (j−l)(2j+1). The quantity Ω = −j, j+1, …, j−1, j is the projection of j on a quantization axis and n = 0, 1, …  denotes the principal quantum number counting states with identical κ and Ω. When the potentials are no longer spherically but axially symmetric, κ is no longer a good quantum number of the eigenstates of the Dirac equation. However, any state can be represented as a sum

$\begin{array}{cc}\begin{array}{l}{\text{Ψ}}_k(\overrightarrow{r})=\frac{1}{r}{\displaystyle \sum _{\kappa }}\left(\begin{array}{c}{F}_{\kappa }(r){{Y}}_{\kappa \Omega }(\widehat{r})\\ i{G}_{\kappa }(r){{Y}}_{-\kappa \Omega }(\widehat{r})\end{array}\right)\end{array}& (47)\end{array}$

with radial wave functions F_κ(r) and G_κ(r) in the upper and lower component of the spinor. Only values of κ with | κ |≥| Ω | +1/2 are allowed. The angular and spin dependence of the state (47) is encoded in the spinor spherical harmonics

$\begin{array}{cc}\begin{array}{l}{{Y}}_{\kappa \Omega }(\widehat{r})={\displaystyle \sum _{m_l{m}_s}}(l{m}_{l}s{m}_s|j\Omega )Y_{l{m}_l}(\widehat{r}){\chi }_{s{m}_s}\end{array}& (48)\end{array}$

with a coupling of the spherical harmonics Y_lml with the spinors χ_sms of spin s = 1/2. Because of the sign change of κ between the upper and lower component, there is a difference of one between the orbital angular momenta. If the potentials V and S in (46) have positive parity, an eigenstate k is characterized by k = (n, π, Ω) and the sequence of components in the sum (47) is limited further to κ = s(|Ω| + 1/2), −s(|Ω| + 3/2), s(|Ω| + 5/2), …with s = (−1)^|Ω|−π/2, and alternating signs for κ. The principal quantum number n is now counting states with identical π and Ω.

The potentials in the Dirac Hamiltonian (46) can be expanded in multipole series

$\begin{array}{cc}V(\overrightarrow{r})={\displaystyle \sum _L{V}_L}(r)P_L\text{ and }S(\overrightarrow{r})={\displaystyle \sum _L{S}_L}(r)P_L& (49)\end{array}$

with Legendre polynomials $P_L=\sqrt{4\pi /(2L+1)}{Y}_{L0}(\widehat{r})$ due to the presumed axial symmetry. Using the identities $$\overrightarrow{\sigma }·\overrightarrow{r}{\mathcal{Y}}_{\kappa \Omega }=-r{\mathcal{Y}}_{-\kappa \Omega }$$ and

$\begin{array}{cc}\overrightarrow{\sigma }·\overrightarrow{p}\frac{f(r)}{r}{\mathcal{Y}}_{\kappa \omega }=\frac{i}{r}\left(\frac{d}{dr}-\frac{\kappa }{r}\right)f(r){\mathcal{Y}}_{-\kappa \Omega },& (50)\end{array}$

the Dirac Equation (12) transforms to a set of coupled first-order differential equations

$\begin{array}{cc}\begin{array}{lll}{\displaystyle \sum _{L,{\kappa }^{\prime }}}(V_L-S_L)C_{\kappa {\kappa }^{\prime }}^{L\Omega }{F}_{{\kappa }^{\prime }}+m{F}_{\kappa }-\left(\frac{d}{dr}+\frac{\kappa }{r}\right)G_{\kappa }& =& E_k{F}_{\kappa }\end{array}& (51)\end{array}$

$\begin{array}{cc}\left(\frac{d}{dr}\text{ }-\text{ }\frac{\kappa }{r}\right)F_{\kappa }-m{G}_{\kappa }\text{}+\text{}{\displaystyle \sum _{L,{\kappa }^{\prime }}\left(V_L\text{}+\text{}{S}_L\right)}{C}_{-\kappa -{\kappa }^{\prime }}^{L\Omega }{G}_{{\kappa }^{\prime }}=\text{ }{E}_k{G}_{\kappa }& (52)\end{array}$

for the radial wave functions after projection on ${{Y}}_{\pm \kappa \Omega }$ states. The factors

$\begin{array}{cc}\begin{array}{l}{C}_{\kappa {\kappa }^{\prime }}^{L\Omega }={\displaystyle \int d}\Omega {\mathcal{Y}}_{\kappa \Omega }^{†}{P}_L{\mathcal{Y}}_{{\kappa }^{\prime }\Omega }\\ ={(-1)}^{s+j^{\prime }+l^{\prime }}\sqrt{(2l^{\prime }+1)(2j^{\prime }+1)}(l^{\prime }0L0|l0)(j^{\prime }\Omega L0|j\Omega )\left\{\begin{array}{ccc}{l}^{\prime }& s& j^{\prime }\\ j& L& l\end{array}\right\}\end{array}& (53)\end{array}$

with Clebsch–Gordan and 6j coefficients satisfy the symmetry $C_{\kappa {\kappa }^{\prime }}^{L\Omega }=C_{-\kappa -{\kappa }^{\prime }}^{L\Omega }$. They have the particular value $C_{\kappa {\kappa }^{\prime }}^{0\Omega }={\delta }_{\kappa {\kappa }^{\prime }}$.

To apply the Lagrange-mesh method, the radial wave functions F_κ and G_κ are represented by vectors of finite dimension N that corresponds to the number of zeroes of the Laguerre polynomial $L_N^{(\alpha )}(x)$. The parameter α is set to 2|κ| so that the basis functions (31) have the correct asymptotic form in the limit r → 0 for the radial wave function in the component with the lower orbital angular momentum. The radial grid points are defined as $r_i^{(\alpha )}=h{x}_i^{(\alpha )}$ for i = 1, …, N with a scale parameter h. Then the vector

$\begin{array}{cc}\begin{array}{l}{X}_{\kappa }=\left(\begin{array}{c}{F}_{\kappa }(r_1^{(\alpha )})\\ ⋮\\ F_{\kappa }(r_N^{(\alpha )})\\ G_{\kappa }(r_1^{(\alpha )})\\ ⋮\\ G_{\kappa }(r_N^{(\alpha )})\end{array}\right)\end{array}& (54)\end{array}$

of dimension 2N for a given channel κ is introduced. The differential Equations (51) and (52) constitute an infinite set of coupled equations that has to be truncated for solving the Dirac equation numerically. Choosing a sequence κ₁, …, κ_Nc of N_c channels, the coupled Equations (51) and (52) become an eigenvalue problem in matrix form

$\begin{array}{cc}\begin{array}{l}\left(\begin{array}{ccc}{H}_{{\kappa }_1{\kappa }_1}& \cdots & H_{{\kappa }_1{\kappa }_{N_c}}\\ ⋮& \ddots & ⋮\\ H_{{\kappa }_{N_c}{\kappa }_1}& \cdots & H_{{\kappa }_{N_c}{\kappa }_{N_c}}\\ \end{array}\right)\left(\begin{array}{c}{X}_{{\kappa }_1}\\ ⋮\\ X_{{\kappa }_{N_c}}\end{array}\right)=E_k\left(\begin{array}{c}{X}_{{\kappa }_1}\\ ⋮\\ X_{{\kappa }_{N_c}}\end{array}\right)\end{array}& (55)\end{array}$

for the state k with quadratic sub-matrices $H_{\kappa {\kappa }^{\prime }}$ of dimension 2N×2N. These have the form

$\begin{array}{cc}\begin{array}{l}{H}_{\kappa {\kappa }^{\prime }}=\left(\begin{array}{cc}{W}_{\kappa {\kappa }^{\prime }}^{(-)}+m{\delta }_{\kappa {\kappa }^{\prime }}& -D_{\kappa }^{(+)}{\delta }_{\kappa {\kappa }^{\prime }}\\ D_{\kappa }^{(-)}{\delta }_{\kappa {\kappa }^{\prime }}& W_{\kappa {\kappa }^{\prime }}^{(+)}-m{\delta }_{\kappa {\kappa }^{\prime }}\end{array}\right)\end{array}& (56)\end{array}$

where $W_{\kappa {\kappa }^{\prime }}^{(\pm )}$ and $D_{\kappa }^{(\pm )}$ are potential and derivative matrices of dimension N×N itself. Since the derivative matrix $D_{\kappa }^{(\pm )}$ only appears in the diagonal sub-matrices of the Hamiltonian matrix in (55), its elements are given by

$\begin{array}{cc}\begin{array}{l}{(D_{\kappa }^{(\pm )})}_{ij}^{(\alpha \alpha )}=\frac{1}{h}{(\frac{d}{dx}\pm \frac{\kappa }{x})}_{ij}^{(\alpha \alpha )}=\frac{1}{h}\left[{(\frac{d}{dx})}_{ij}^{(\alpha \alpha )}\pm \frac{\kappa }{x_i^{(\alpha )}}{\delta }_{ij}\right]\end{array}& (57)\end{array}$

with the representations (43) and (44) of the first derivative. Hence they are analytic in the Lagrange-mesh method. When solving the Schrödinger equation to obtain wave functions using non-relativistic EDFs, an analytic representation of the second derivative would occur instead of the first derivative here. The potential matrices $W_{\kappa {\kappa }^{\prime }}^{(\pm )}$ connect components of the total wave function with equal and different values of κ. Hence, their elements are given in general by

$\begin{array}{cc}\begin{array}{l}{\left(W_{\kappa {\kappa }^{\prime }}^{(\pm )}\right)}_{ij}^{(\alpha {\alpha }^{\prime })}=\sqrt{\frac{{\lambda }_j^{({\alpha }^{\prime })}}{{\lambda }_i^{(\alpha )}}}\left({\displaystyle \sum _L}\left[V_L(r_j^{({\alpha }^{\prime })})\pm S_L(r_j^{({\alpha }^{\prime })})\right]C_{\kappa {\kappa }^{\prime }}^{L\Omega }\right)U_{ji}^{({\alpha }^{\prime }\alpha )}\end{array}& (58)\end{array}$

using the transformation rule (40). In total, a d × d matrix with dimension d = 2N × N_c has to be diagonalized where N_c is the number of channels characterized by different values of κ that are considered. This is achieved numerically by applying a householder reduction of the real, symmetric matrix in (55) to tridiagonal form and subsequently the QL algorithm with implicit shifts to determine eigenvalues and eigenvectors as described in Press et al. [43]. The spectrum of eigenvalues E_k decomposes into two sets with positive and negative values, respectively. Only those that are positive with E_k < m correspond to the bound states relevant in the description of nuclei. For axially symmetric systems, there is an energetic degeneracy of the single-particle states (n, π, Ω) and (n, π, −Ω) and every level of this form can be populated with one particle at most due to the Pauli principle, i.e., the occupation factor w_ik can be 0 or 1.

The wave functions serve as an input to solve the field Equations (7)–(11) of the mesons and the photon. The source densities are obtained from the total scalar and vector densities according to Equations (4) and (5) with the single-particle densities. The single-particle vector and scalar densities (2) and (3) itself can be expanded in multipole series

$\begin{array}{cc}\begin{array}{l}{n}_k^{(v)}(\overrightarrow{r})={\displaystyle \sum _L{n}_{kL}^{(v)}(r)Y_{L0}(\widehat{r})\mathrm{\text{and\hspace{1em} }}}{n}_k^{(s)}(\overrightarrow{r})={\displaystyle \sum _L}{n}_{kL}^{(s)}(r)Y_{L0}(\widehat{r})\end{array}& (59)\end{array}$

with spherical harmonics and radial functions

$\begin{array}{cc}\begin{array}{lll}{n}_{kL}^{(v)}(r)& =& \frac{1}{r^2}\sqrt{\frac{2L+1}{4\pi }}{\displaystyle \sum _{\kappa {\kappa }^{\prime }}}\left[C_{\kappa {\kappa }^{\prime }}^{L\Omega }{F}_{\kappa }^{*}(r)F_{{\kappa }^{\prime }}(r)+C_{-\kappa -{\kappa }^{\prime }}^{L\Omega }{G}_{\kappa }^{*}(r)G_{{\kappa }^{\prime }}(r)\right]\end{array}& (60)\end{array}$

$\begin{array}{cc}\begin{array}{lll}{n}_{kL}^{(s)}(r)& =& \frac{1}{r^2}\sqrt{\frac{2L+1}{4\pi }}{\displaystyle \sum _{\kappa {\kappa }^{\prime }}}\left[C_{\kappa {\kappa }^{\prime }}^{L\Omega }{F}_{\kappa }^{*}(r)F_{{\kappa }^{\prime }}(r)-C_{-\kappa -{\kappa }^{\prime }}^{L\Omega }{G}_{\kappa }^{*}(r)G_{{\kappa }^{\prime }}(r)\right]\end{array}& (61)\end{array}$

that are easily mapped to any given radial grid using the representation with the Lagrange-mesh functions (24). The wave functions are normalized as

$\begin{array}{cc}\begin{array}{l}1=\int d^{3}r{n}_k^{(v)}(\overrightarrow{r})={\displaystyle \sum _{\kappa }}{\displaystyle \sum _{i=1}^N}{\lambda }_i^{(\alpha )}\left[|F_{\kappa }(r_i^{(\alpha )}){|}^2+|G_{\kappa }(r_i^{(\alpha )}){|}^2\right]\end{array}& (62)\end{array}$

using the integral (30).

3.3. Solving the Meson-Field Equations

The field Equations (7)–(11) of mesons and photons have the generic form

$\begin{array}{cc}\begin{array}{l}-\Delta A(\overrightarrow{r})+m^{2}A(\overrightarrow{r})=S(\overrightarrow{r})\end{array}& (63)\end{array}$

where the mass m is zero in case of the electromagnetic potential. Since the system is assumed to be axially symmetric, the series expansions

$\begin{array}{cc}\begin{array}{l}A(\overrightarrow{r})=\frac{1}{r}{\displaystyle \sum _{L=0}^{\infty }{a}_L(r)Y_{L0}(\widehat{r})\mathrm{\text{and\hspace{1em} }}}S(\overrightarrow{r})=\frac{1}{r}{\displaystyle \sum _{L=0}^{\infty }}{s}_L(r)Y_{L0}(\widehat{r})\end{array}& (64)\end{array}$

can be introduced. After the separation of radial and angular coordinates, the inhomogeneous second-order differential equation

$\begin{array}{cc}\begin{array}{l}\left[-\frac{d^2}{d{r}^2}+\frac{L(L+1)}{r^2}+m^2\right]a_L(r)=s_L(r)\end{array}& (65)\end{array}$

for the radial functions a_L(r) is obtained for all angular momenta L. The source functions s_L(r) have to be determined by numerical angular integration

$\begin{array}{cc}\begin{array}{l}{s}_L(r)=\int d\Omega S(\overrightarrow{r})Y_{L0}(\widehat{r})\end{array}& (66)\end{array}$

for the meson fields because the source term $S(\overrightarrow{r})$ in (63) is in general a product of a coupling and a source density that both depend on the position in space. The integration over the azimuthal angle φ is trivial and gives a factor 2π. The integration over the polar angle ϑ is performed using a Gauss-Legendre integration rule with N_ϑ = 23 grid points that is sufficient in the applications of this work. Only for the electromagnetic field with constant coupling Γ_γ the angular integration is analytic with $s_L(r)=r{\Gamma }_{\gamma }{n}_{\gamma }^{(L)}(r)$ for $n_{\gamma }(\overrightarrow{r})=\sum _{L=0}^{\infty }{n}_{\gamma }^{(L)}(r)Y_{L0}(\widehat{r})$.

The source functions $S(\overrightarrow{r})$ are rapidly decreasing for radii r that are sufficiently larger than the radius of the nucleus under consideration. They can be assumed to be zero for r larger than a fixed radius R. Since the interaction described by the mesons is short-ranged, the functions a_L(r) can also be considered to vanish for radii beyond R. Hence it is convenient to introduce the expansion

$\begin{array}{cc}\begin{array}{l}{a}_L(r)={\displaystyle \sum _{n=1}^{\infty }}{a}_n^{(L)}{N}_n^{(L)}{z}_L(y_n^{(L)}r/R)\end{array}& (67)\end{array}$

and similar for s_L(r) on the interval [0, R] with coefficients $a_n^{(L)}$ ($s_n^{(L)}$) and Riccati–Bessel functions z_L(y) = yj_L(y). The quantities $y_n^{(L)}$ are the zeroes of the spherical Bessel functions j_L(y) and

$\begin{array}{cc}{N}_n^{(L)}=\sqrt{\frac{2}{R}}\times \{\begin{array}{lll}{\left|z_{L-1}(y_n^{(L)})\right|}^{-1}\hfill & \text{if}\hfill & L\ge 1\hfill \\ 1\hfill & \text{if}\hfill & L=0\hfill \end{array}& (68)\end{array}$

are normalization factors. After inserting this series into the differential Equation (65), utilizing the differential equation of the Riccati–Bessel functions and the orthogonality relation [42]

$\begin{array}{cc}{\displaystyle {\int }_0^{1}d}t{z}_L(y_n^{(L)}t)z_L(y_m^{(L)}t)=\frac{1}{2}{\left[z_{L-1}(y_n^{(L)})\right]}^2{\delta }_{nm},& (69)\end{array}$

the coefficients in (67) are immediately found as

$\begin{array}{cc}\begin{array}{l}{a}_n^{(L)}=s_n^{(L)}{\left[m^2+{(\frac{y_n^{(L)}}{R})}^2\right]}^{-1}\end{array}& (70)\end{array}$

with

$\begin{array}{cc}{s}_n^{(L)}=N_n^{(L)}{\displaystyle {\int }_0^{R}d}r{s}_L(r)z_L(y_n^{(L)}r/R)& (71)\end{array}$

that have to be calculated numerically.

A similar method can be used for the electromagnetic field, however, with a few modifications. Since the mass m is zero, the terms with $y_n^{(L)}=0$ are excluded from the sum (67). The obtained function a_L(r) vanishes at r = R and a correction has to be implemented by adding a contribution b_L(r) that gives the correct asymptotic behavior of the long-range field, i.e., the electromagnetic potential is written as

$\begin{array}{cc}\begin{array}{l}A(\overrightarrow{r})=\frac{1}{r}{\displaystyle \sum _{L=0}^{\infty }}\left[a_L(r)+b_L(r)\right]Y_{L0}(\widehat{r}).\end{array}& (72)\end{array}$

The additional term is given by

$\begin{array}{cc}{b}_L(r)=\{\begin{array}{ccc}{c}_L{r}^{L+1}& \text{if}& r\le R\\ c_L{r}^{-L}& \text{if}& r>R\end{array}& (73)\end{array}$

with a coefficient

$\begin{array}{cc}{c}_L=\frac{1}{2L+1}{\displaystyle {\int }_0^{R}d}rr{s}_L(r)& (74)\end{array}$

that is related to multipole moments of the charge distribution, e.g., $c_0={\Gamma }_{\gamma }^{2}Z/\sqrt{4\pi }$ with the charge Z of the nucleus.

4. Results

The numerical methods that were introduced in the preceding section are applied to two cases in this work. First, the motion of a fermion in a deformed harmonic oscillator potential is considered to study the numerical accuracy of solving the Dirac equation with the Lagrange-mesh method in spherical coordinates. The relativistic harmonic oscillator was often considered in the connection with pseudo-spin symmetry in different versions [44, 45–47]. Second, the ground states of self-conjugate even-even light nuclei are calculated by a self-consistent solution of the Dirac equation of the nucleons and the field equations for mesons and the photon. Light self-conjugate nuclei are well-studied theoretically, in particular regarding the occurrence of α-particle clusters [48–50].

4.1. Deformed Harmonic Oscillator Potential

Using simple vector and scalar potentials in the Hamiltonian (46), it is possible to find analytic solutions of the Dirac equation. A comparison with numerical results allows to estimate how the accuracy of the calculation depends on the parameters of the method. A particularly simple case is the three-dimensional harmonic oscillator of mass m with different oscillator constants ω_x, ω_y, and ω_z that define the potentials

$\begin{array}{cc}\begin{array}{l}{U}_x=\frac{m}{2}{\omega }_x^2{x}^2{U}_y=\frac{m}{2}{\omega }_y^2{y}^2\mathrm{\text{and\hspace{1em} }}{U}_z=\frac{m}{2}{\omega }_z^2{z}^2\end{array}& (75)\end{array}$

for the spatial coordinates x, y, and z, respectively, in a Cartesian system. Setting

$\begin{array}{cc}\begin{array}{l}V=-S=\frac{1}{2}U=\frac{1}{2}(U_x+U_y+U_z),\end{array}& (76)\end{array}$

the Dirac equation can be written as

$\begin{array}{cc}\begin{array}{l}H\text{Ψ}=\left(\begin{array}{cc}U+m& \overrightarrow{\sigma }·\overrightarrow{p}\\ \overrightarrow{\sigma }·\overrightarrow{p}& -m\end{array}\right)\text{Ψ}=E\text{Ψ}\end{array}& (77)\end{array}$

with the potential U appearing only in the upper left diagonal position. Introducing upper and lower components φ and χ in the wave function as

$\begin{array}{cc}\begin{array}{l}\text{Ψ}=\left(\begin{array}{c}\phi \\ \chi \end{array}\right),\end{array}& (78)\end{array}$

the partial differential equation

$\begin{array}{cc}\begin{array}{l}(\frac{p^2}{E+m}+U)\phi =(E-m)\phi \end{array}& (79)\end{array}$

is found for the upper component by decoupling using

$\begin{array}{cc}\begin{array}{l}\chi =\frac{\overrightarrow{\sigma }·\overrightarrow{p}}{E+m}\phi \end{array}& (80)\end{array}$

and $\overrightarrow{\sigma }·\overrightarrow{p}\overrightarrow{\sigma }·\overrightarrow{p}=p^2$. A separation of variables with φ = φ_x(x)φ_y(y)φ_z(z) and E−m = E_x + E_y + E_z gives the simple ordinary second-order differential equation

$\begin{array}{cc}\begin{array}{l}\left(\frac{1}{2}m{\omega }_x^2+\frac{p_x^2}{E+m}\right){\phi }_x=E_x{\phi }_x\end{array}& (81)\end{array}$

for the oscillation in x direction and similarly for the y and z oscillators. This can be compared with the corresponding equation of a one-dimensional non-relativistic harmonic oscillator

$\begin{array}{cc}\begin{array}{l}(\frac{1}{2}M{\Omega }^2+\frac{p^2}{2M})\Phi =\epsilon \Phi \end{array}& (82)\end{array}$

with mass M and oscillator frequency Ω. It has the eigenenergies ε_n = (n+1/2)ℏΩ with n = 0, 1, 2, … . Identifying $m{\omega }_x^2=M{\Omega }^2$, E+m = 2M, and E_x = ε, the energy eigenvalues E_nxnynz of the three dimensional relativistic oscillator can be found from the implicit equation

$\begin{array}{cc}\begin{array}{l}{E}_{n_x{n}_y{n}_z}-m=\sqrt{\frac{2m}{E_{n_x{n}_y{n}_z}+m}}[\left(n_x+\frac{1}{2}\right)\hslash {\omega }_x+\left(n_y+\frac{1}{2}\right)\hslash {\omega }_y\\ +\left(n_z+\frac{1}{2}\right)\hslash {\omega }_z]\end{array}& (83)\end{array}$

with the three quantum numbers n_x, n_y, and n_z that are non-negative integers. In fact, (83) is a cubic equation for E_nxnynz with only one positive solution that is easily found by iteration.

For the application of the Lagrange-mesh method, it is convenient to transform the potential U = V−S from the Cartesian to the spherical basis. Assuming axial symmetry with ω_⊥ = ω_x = ω_y and ω_|| = ω_z, the potential U can be written as

$\begin{array}{cc}\begin{array}{l}U=\frac{m}{2}\left[{\omega }_{\perp }^2(x^2+y^2)+{\omega }_{\left|\right|}^2{z}^2\right]=\frac{m}{2}\frac{4\pi }{3}({\omega }_{\left|\right|}^2{Y}_{10}^2-2{\omega }_{\perp }^2{Y}_{11}{Y}_{1-1})r^2\end{array}& (84)\end{array}$

where the definition of the spherical harmonics Y_1μ was used. Considering the multiplication theorem for spherical harmonics and the parametrization

$\begin{array}{cc}\begin{array}{l}{\omega }_{\perp }=\omega (1+d)\mathrm{\text{and\hspace{1em} }}{\omega }_{\left|\right|}=\omega (1-2d)\end{array}& (85)\end{array}$

with frequency ω and deformation parameter d, the expression for the potential further reduces to

$\begin{array}{cc}\begin{array}{l}U=\frac{m}{2}{\omega }^2{r}^2\left[(1+2d^2)P_0-2d(2-d)P_2\right]\end{array}& (86)\end{array}$

with Legendre polynomials P_λ. The contribution in U proportional to P₂ is responsible for the coupling of different angular momentum states in the wave function.

The test of the numerical method is performed for an oscillator with ℏω = 20 MeV and mass m = 938 MeV/c². The corresponding oscillator parameter $b=\sqrt{\hslash /(m\omega )}{\left[2m/(E+m)\right]}^{1/4}\approx 1.44$ fm defines the characteristic length scale of the system. In Table 1, the energies $E_{n_x{n}_y{n}_z}^{\text{rel}}-m$ of a relativistic spherical harmonic oscillator with these parameters are compared with those of a non-relativistic oscillator with $E_{n_x{n}_y{n}_z}^{\text{non}-\text{rel}}-m=(n+3/2)\hslash \omega =(20n+30)$ MeV. The energies only depend on the total oscillator quantum number n = n_x+n_y+n_z because of the degeneracy of the spherical system. This degeneracy will be lifted for non-vanishing deformation parameter d. The energies in the relativistic case are lower than those in the non-relativistic case and the relative deviation

$\begin{array}{cc}\begin{array}{l}{\delta }_{n_x{n}_y{n}_z}=\frac{E_{n_x{n}_y{n}_z}^{\text{rel}}-m}{E_{n_x{n}_y{n}_z}^{\text{non}-\text{rel}}-m}-1\end{array}& (87)\end{array}$

is sizable. It is in the order of a few percent and increases with n as can be seen from Table 1.

TABLE 1

Table 1. Energies of the lowest eigenstates of a spherically symmetric harmonic oscillator with -hω = 20 MeV and m = 938 MeV/c² in the non-relativistic and relativistic case and relative difference (87).

In general, numerical solutions of the Dirac equation depend on the deformation parameter d and three parameters that define the Lagrange-mesh method: the scale parameter h, the number of zeroes N of the Laguerre polynomial in (24) and (25), and the number of channels N_c in the expansion of the wave function as in (55). As baseline for the following calculations and comparisons, the values h = 0.1 fm, N = 25, and N_c = 11 are chosen. Because the coupling potential has positive parity, there are two distinct sequences of κ values in the expansion (47) of the single-particle state that do not couple as discussed in section 3.2. They are {1, −2, 3, −4, …} and {−1, 2, −3, 4, …} for states of positive and negative parity, respectively. From these sets N_c consecutive values κ₁, …, κ_Nc are chosen as the considered channels with the lowest value |κ₁| = |Ω|+1/2. In Tables 2, 3, the exact energies $E_{n_x{n}_y{n}_z}^{\text{exact}}-m$ of the deformed relativistic harmonic oscillator with deformation parameters d = 0.1 and d = −0.1, respectively, are given for the lowest states with n_x+n_y+n_z ≤ 3. The corresponding parities π, principal quantum numbers n, and angular momentum projections Ω are also shown. A degeneracy of states with identical values for n_x+n_y and n_z is clearly seen because of the rotational symmetry around the z axis. The last column in Tables 2, 3 gives the absolute deviations

$\begin{array}{cc}\begin{array}{l}{\Delta }_{n_x{n}_y{n}_z}=|E_{n_x{n}_y{n}_z}^{\text{LMM}}-E_{n_x{n}_y{n}_z}^{\text{exact}}|\end{array}& (88)\end{array}$

of the energies $E_{n_x{n}_y{n}_z}^{\text{LMM}}$, which were determined with the Lagrange-mesh method using the baseline parameters, from the exact values $E_{n_x{n}_y{n}_z}^{\text{exact}}$, which are shown in the next-to-last column. The highest accuracy is found for the states of lowest energy. It decreases with increasing principal quantum number n for fixed |Ω| and π. Nevertheless, the level of accuracy is impressive as the absolute scale of the energies in the eigenvalue problem is given by the mass m, i.e., of the order of 1 GeV.

TABLE 2

Table 2. Exact energies of the eigenstates of a prolate relativistic harmonic oscillator with d = 0.1, -hω = 20 MeV, and m = 938 MeV/c² and absolute deviation (88) of a calculation using the Lagrange-mesh method with scaling factor h = 0.1 fm, number of grid points N = 25, and number of channels N_c = 11.

TABLE 3

Table 3. Exact energies of the eigenstates of an oblate relativistic harmonic oscillator with d = −0.1, -hω = 20 MeV, and m = 938 MeV/c² and absolute deviation (88) of a calculation using the Lagrange-mesh method with scaling factor h = 0.1 fm, number of grid points N = 25, and number of channels N_c = 11.

The dependence of the accuracy of the calculation on the parameters h, N, and N_c will be studied in the following sections. The parameters h and N are specific to the Lagrange-mesh method since they define the position of the discretization points for the radius. In contrast, the parameter N_c appears only because a partial-wave expansion of the three-dimensional wave fuctions was introduced in the present work. Hence, the convergence with respect to h and N is rather independent from that of N_c. Alternatively, different versions of the Lagrange-mesh method based on other basis functions using cylindrical or Cartesian coordinates could be used.

4.1.1. Variation of the Scale Parameter

The scale parameter h connects the radial grid points $r_i=h{x}_i^{(\alpha )}$ in a given channel κ with α = 2|κ| to the zeroes $x_i^{(\alpha )}$ of the Laguerre polynomial $L_N^{(\alpha )}(x)$. It has to be selected so that the details and the gross structure of the wave functions can be well represented. The effect of changing h can be clearly demonstrated by looking at the radial wave function of the ground and first excited state with orbital angular momentum l = 0 of a spherical relativistic harmonic oscillator. They are depicted in the left and right panel of Figure 1, respectively. In this case, only a single channel with κ = 1 has to be considered and the matrix Hamiltonian in (55) has dimension 2N × 2N with N = 25 radial grid points. The exact result for the upper component is given by a simple Gaussian φ(r) = [2/(π^1/4b^3/2)]exp[−r²/(2b²)] for the ground state and φ(r) = (3/2−r²/b²)[2^3/2/(3^1/2π^1/4b^3/2)]exp[−r²/(2b²)] for the first excited state with oscillator parameter b. Of course, the choice of Lagrange-mesh functions (20) with an exponential decay for large radii is not very well adapted to represent Gaussian-type functions and some deviation can be expected. In Figure 1, the modulus |Cφ(r)| of the exact and numerically determined wave functions are depicted. The factor C is determined such that the radial wave function of the upper component is normalized to one. Hence C = 1 for the exact result and C = [∫dr r²|φ(r)|²]^−1/2 with φ(r) = F_κ(r)/r in case of the numerical calculation. For small radii r, the exact wave function (full line) practically coincides with the numerically determined wave function if the scale parameter h varies between 0.05 fm and 0.80 fm. Deviations show up for radii larger than approximately 4 fm if h = 0.05 (red squares). This fact is not surprising because the highest grid point is located at $h{x}_N^{(\alpha )}\approx 4.45$ fm in this case and the Gaussian decrease of the exact wave function cannot be described well with the asymptotic form of the Lagrange-mesh functions. For h = 0.10 fm, the exact wave function is very well reproduced up to about 9.5 fm in the calculation (green circles) with a value of about 10⁻¹⁰ fm^−3/2 for the ground state and 10⁻⁹ fm^−3/2 for the first excited state, much smaller than that at the center. For larger values of h, the tail of the wave function develops an oscillatory behavior indicating that the location of the grid points is not chosen properly and the number of grid points within the major part of the wave function is too small. For example the highest grid point is at a radius of approximately 71.2 fm for h = 0.8 fm, much above the relevant range, and there are 15 points beyond 10 fm, leaving only 10 points for the region of lower radii since N = 25 in the baseline case. The tails of the wave functions of the ground and first excited state show a very similar pattern, except that the values of the excited state are about one to two order of magnitude larger. In addition, the node in this wave function is clearly visible.

FIGURE 1

Figure 1. Modulus of the scaled radial wave function of the ground (A) and first excited (B) state with orbital angular momentum l = 0 in a spherically symmetric relativistic harmonic oscillator with ℏω = 20 MeV and m = 938 MeV/c² in the exact case (full black line) and the calculation with the Lagrange-mesh method for different scale parameters h (colored symbols) and constant number of grid points (N = 25). The normalization factor C is explained in the text.

The dependence on the scale parameter can also be studied by looking at the deviations (88) between the calculated energies and the exact energies. They are depicted in Figure 2 for selected values of h between 0.05 fm and 1.00 fm and deformation parameters d = 0.1, d = 0.0, and d = −0.1 in (85). Only states with n_x + n_y + n_z ≤ 2 are considered. If h is in the range of 0.08 fm to 0.3 fm, there is almost no change in the accuracy of the calculation. Only for h values outside this interval, a deterioration is observed with stronger effects at low h. In the deformed calculations also, a systematic trend with the number of excited oscillator quanta is found. The energies of higher-excited states are less well reproduced. This is related to the fact that several channels are needed for a good representation of the wave functions; see section 4.1.3. For the spherical case, the effect is much less obvious because here a single κ channel is sufficient for the calculation.

FIGURE 2

Figure 2. Deviation (88) of the calculated lowest energies of a relativistic harmonic oscillator with ℏω = 20 MeV and m = 938 MeV/c² from the exact results for constant number of grid points (N = 25) and constant number of channels (N_c = 11). Results are shown as a function of the scaling factor h for the spherical case (full symbols), prolate deformation with d = 0.1 (open diamonds), and oblate deformation with d = −0.1 (open squares).

4.1.2. Variation of the Number of Grid Points

The choice of the order N of the Laguerre polynomial in (24) and (25) determines the number of radial grid points of the mesh. The effect of changing N can again be demonstrated by looking at the radial wave functions of the ground and excited states with l = 0 of a spherical relativistic harmonic oscillator in comparison to the exact result as in section 4.1.1. In Figure 3, the modulus |Cφ(r)| of the exact and numerically determined wave functions are depicted for constant scale parameter h = 0.1 fm varying the number of grid points. The ground and first excited states for orbital angular momentum l = 0 are considered in panels A and B, respectively. A clear improvement in the approximation of the exact wave function is found with increasing N with an increase of the radial range where the wave function is well represented. N = 10 grid points are obviously not sufficient for a good representation of the wave functions.

FIGURE 3

Figure 3. Modulus of the scaled radial wave function of the ground (A) and first excited (B) state with orbital angular momentum l = 0 in a spherically symmetric relativistic harmonic oscillator with ℏω = 20 MeV and m = 938 MeV/c² in the exact case (full black line) and the calculation with the Lagrange-mesh method for different number of grid points N (colored symbols) and constant scale parameter (h = 0.10 fm). The normalization factor C is explained in the text.

Keeping h = 0.1 fm and N_c = 11 fixed, the deviations (88) of the energies are depicted in Figure 4 as a function of the number of grid points N in each channel. The deviation is largest for the lowest values of N, and no accurate results for the single-particle energies are obtained from the application of the Lagrange-mesh method. The accuracy increases with N, and above approximately N = 20, a stabilization of the deviations occurs. A further increase of the number of grid points does not result in smaller differences between exact and numerical eigenvalues.

FIGURE 4

Figure 4. Deviation (88) of the calculated lowest energies of a relativistic harmonic oscillator with ℏω = 20 MeV and m = 938 MeV/c² from the exact results for constant scaling factor (h = 0.1 fm) and constant number of channels (N_c = 11). Results are shown as a function of the number of grid points N for the spherical case (full symbols), prolate deformation with d = 0.1 (open diamonds), and oblate deformation with d = −0.1 (open squares).

4.1.3. Variation of the Number of Channels

The third parameter that is relevant in the application of the Lagrange-mesh method is the number of channels N_c in the expansion of the single-particle wave function (47). In Figure 5, the deviations (88) of the energies of the lowest states are shown for values of N_c ranging from 4 to 14. As in Figures 2, 4, the relativistic harmonic oscillator with parameters ℏω = 20 MeV and m = 938 MeV/c² is chosen as test case with the baseline values for h and N. The spherical harmonic oscillator displays no dependence of Δ_nxnynz on N_c because there is no coupling between different κ channels and the accuracy of the calculation only depends on h and N as discussed in sections 4.1.1 and 4.1.2. In contrast to this case, the accuracy of the single-particle energies exhibits a clear trend of improvement with increasing number of channels in the calculations for harmonic oscillator potentials with prolate and oblate deformations. The ground state energy seems to converge for N_c reaching the baseline value but an even larger number of channels are required to obtain a similar accuracy in the calculation of excited states.

FIGURE 5

Figure 5. Deviation (88) of the calculated lowest energies of a relativistic harmonic oscillator with ℏω = 20 MeV and m = 938 MeV/c² from the exact results for constant scaling factor (h = 0.1 fm) and constant number of grid points (N = 25). Results are shown as a function of the number of channels N_c for the spherical case (full symbols), prolate deformation with d = 0.1 (open diamonds), and oblate deformation with d = −0.1 (open squares).

4.2. Deformed Nuclei

In the description of nuclei that employs a relativistic EDF, two coupled problems have to be solved self-consistently. The single-particle wave functions of the nucleons are found as solutions of the Dirac Equation (12) with vector and scalar potentials that are determined by the meson and electromagnetic fields. These are calculated with the help of the field Equations (7)–(11) from the source densities (4) and (5) that in turn depend on the single-particle wave functions. An iterative approach is usually selected with a predefined accuracy that has to be reached for some quantities.

In the present application, the numerical methods of section 3 for the description in spherical coordinates are employed. The iteration of the coupled equations is stopped when the single-particle energies of the occupied levels do not change more than 10⁻⁶ MeV in consecutive steps. The number of grid points N in the Lagrange-mesh method and the number of channels N_c are set to the baseline values, i.e., N = 25 and N_c = 11. In contrast, the scale parameter h changes from nucleus to nucleus. With the radius estimate R = 1.2·A^1/3 fm for a nucleus with mass number A and the largest zero $x_N^{(\alpha )}$ for the Laguerre polynomial $L_N^{(\alpha )}(x)$ with α = 0, the scale parameter is defined as $h=0.95(R+10\mathrm{\text{fm}})/x_{25}^{(0)}$ with $x_{25}^{(0)}\approx 85.26$. This choice provides a suitable discretization of the wave functions and densities in the range that is relevant for the calculation including the asymptotic behavior at large radii. Since a defined parity is imposed on the single-particle wave functions, only even angular momenta appear in the expansions of potentials (49), densities (59), fields, and sources (63). Other possibilities, e.g., octupole distributions, that would require to include odd angular momenta, are not considered here.

The EDF selected in this work is based on a relativistic mean-field model with density-dependent meson-nucleon couplings [17]. It has been applied successfully with different parametrizations in the description of nuclear structure and nuclear matter. Here, the parameter set DD2 [51] is used; this includes σ, ω, and ρ but not δ mesons to describe the effective interaction. The masses of neutrons and protons are fixed to m_n = 939.565379 MeV/c² and m_p = 938.272046 MeV/c², respectively. The nuclear matter parameters calculated with this parametrization are very reasonable and consistent with present constraints [52].

As examples for the application of the Lagrange-mesh method, the ground states of self-conjugate even-even nuclei from ⁴He to ⁴⁰Ca are calculated with mass numbers that are multiples of four. They show large structural differences and the nucleon distribution changes between spherical, oblate, and prolate shapes. The rather small number of nucleons allows a fast calculation and the choice of the baseline parameters N = 25 and N_c = 11 in the Lagrange-mesh method is sufficient to obtain faithful results. It is expected that the exponential decrease of the radial Lagrange-mesh functions gives a better approximation of the single-particle wave functions in a nucleus than in a harmonic oscillator. The obtained accuracy is fully adequate for a representation of single-particle energies and density distributions in the manner of the figures in this work. The matter radii and deformation parameters are expected to change at most by one unit in the last digit given in Table 4 if the precision is improved by increasing the number channels N_c as trial calculations have shown.

TABLE 4

Table 4. Point matter radius r_matter and deformation parameters β₂, …, β₁₀ for nuclei ⁴He, …, ⁴⁰Ca.

In Figures 6, 7, the single-particle energies of neutrons and protons are shown, respectively. Each circle represents a pair of states (n, π, Ω) and (n, π, −Ω). Degenerate levels are clearly identified in the spherical closed-shell nuclei ⁴He, ¹⁶O, and ⁴⁰Ca where the orbital angular momentum is an additional good quantum number that characterizes a level. These levels belong to p, d, and f states. There is no such degeneracy found in the remaining deformed nuclei. Owing to the Coulomb repulsion, the proton levels are less bound than the neutron levels with increasing charge number of the nucleus. Shell gaps are easily identified, in particular, for the less strongly deformed nuclei.

FIGURE 6

Figure 6. Neutron levels in the self-conjugated even-even nuclei ⁴He - ⁴⁰Ca. Occupied levels (with two neutrons) are denoted by a full red circle and unoccupied levels by an open blue circle.

FIGURE 7

Figure 7. Proton levels in the self-conjugated even-even nuclei ⁴He - ⁴⁰Ca. Occupied levels (with two protons) are denoted by a full red circle and unoccupied levels by an open blue circle.

The main observable considered here is the point nuclear matter distribution in the intrinsic frame of the nucleus. It is obtained by adding the vector densities as given in (59) of all occupied single-particle states. This intrinsic density depends on the three Cartesian coordinates xyz in the body-fixed frame of the nucleus. The nuclear matter distributions in the laboratory frame, obtained from a projection on good total angular momentum J, would be spherical since only even-even nuclei in their ground state with J = 0 are considered here. Instead of showing radial density distributions that would depend on the the polar angle, intrinsic matter distributions are depicted for all 10 nuclei in Figures 8–11 using lines of constant density in the xz plane to show the deformation clearly. Owing to the axial symmetry, the density distributions are invariant for a rotation around the z axis and hence the lines of constant density in the yz plane would look identical to those in the xz plane or any plane that contains the z axis. As expected for doubly-magic nuclei, ⁴He, ¹⁶O, and ⁴⁰Ca are spherically symmetric. In the case of ⁴He and ⁴⁰Ca, the maximum density is at the center of the nucleus, whereas ¹⁶O exhibits a central depression; see Figure 9. The nuclei ⁸Be, ²⁰Ne, ²⁴Mg, and ³²S show a prolate point matter distribution, however, with distinct differences. A pronounced two-alpha cluster structure is visible in ⁸Be. A pronounced bulge in the xy plane is found in ²⁰Ne but not in ²⁴Mg. The prolate deformation is less obvious for ³²S and the density maxima lie even in the xy plane and not on the z axis as for the other prolate nuclei. An oblate density distribution of the ground state is predicted for ¹²C, ²⁸Si, and ³⁶Ar. A bulge along the z axis is obtained for ²⁸Si and density maxima on this symmetry axis are seen for ³⁶Ar in contrast to ¹²C and ²⁸Si. Looking at the sequence of nuclei with increasing A, an increasing extension of the density distribution can be also traced that will be reflected in the matter root-mean-square radius.

FIGURE 8

Figure 8. Lines of constant intrinsic matter density ρ of the nuclei ⁴He (A) and ⁸Be (B). Points of maximum density are indicated with a × .

FIGURE 9

Figure 9. Lines of constant intrinsic matter density ρ of the nuclei ¹²C (A), ¹⁶O (B), and ²⁰Ne (C). See Figure 8 for the legend. Points of maximum density are indicated with a × or a full black line in case of ¹⁶O.

FIGURE 10

Figure 10. Lines of constant intrinsic matter density ρ of the nuclei ²⁴Mg (A), ²⁸Si (B), and ³²S (C). See Figure 8 for the legend. Points of maximum density are indicated with a × .

FIGURE 11

Figure 11. Lines of constant intrinsic matter density ρ of the nuclei ³⁶Ar (A) and ⁴⁰Ca (B). See Figure 8 for the legend. Points of maximum density are indicated with a × .

In Table 4, the point matter root-mean-square radii r_matter of the nuclei are explicitly given in the second column. They are calculated as

$\begin{array}{cc}\begin{array}{l}{r}_{\text{matter}}={\left[\frac{1}{A}\langle r^2\rangle \right]}^{\frac{1}{2}}={\left[\frac{1}{A}\int d^{3}r{r}^2{\displaystyle \sum _{i=n,p}}{n}_i^{(v)}(\overrightarrow{r})\right]}^{\frac{1}{2}}\end{array}& (89)\end{array}$

with the vector densities (2) of neutrons and protons. The radii r_matter increase with the mass number A except for the step from A = 8 to A = 12 because the nucleus ⁸Be has a particular large radius due to the pronounced cluster structure with two α particles. In reality, ⁸Be is even unbound for the decay into two ⁴He. It is only a resonance in this channel at an energy of 91.84 keV [53]. Besides the point matter radii, Table 4 also gives the deformation parameters β_L for L = 0, 2, 4, 6, 8, and 10. They are defined as

$\begin{array}{cc}\begin{array}{l}{\beta }_L=\frac{4\pi }{3A{R}^L}\langle r^L{Y}_{L0}(\widehat{r})\rangle \end{array}& (90)\end{array}$

with factors A and R in the denominator to obtain a dimensionless quantity that is approximately independent of the size of the nucleus and thus comparable between different nuclei. The radius estimate R = 1.2·A^1/3 fm is defined as above in the determination of the scale parameter h. The value of β_L for all L vanishes for the spherical nuclei ⁴He, ¹⁶O, and ⁴⁰Ca. Prolate nuclei have a positive β₂ and oblate nuclei have a negative β₂ but the signs of β_L with L > 2 are not strictly correlated with the occurrence of a prolate or oblate deformation. The moduli of the deformation parameters usually become smaller with increasing L. The decrease is fast in the case of ¹²C, ³²S, and ³⁶Ar where higher-order deviations from a quadrupole distribution are not very significant. ²⁸Si is an intermediate case and the most extreme variation of the distribution is found for ⁸Be and ²⁰Ne. Substantial values of β_L for large L indicate that channels with large differences of the κ values couple strongly in the diagonalization of the single-nucleon Dirac Hamiltonian. Hence a large number of channels is required for the convergence of the results.

Binding energies of the nuclei are not given here since the energy (16) has to be corrected by a proper projection on center-of-mass momentum and specific total angular momentum or by applying approximate treatments. Methods to perform these corrections will be discussed elsewhere. The angular momentum projection will profit from the representation of the single-particle wave functions in spherical coordinates as proposed in this work.

5. Conclusion

The Lagrange-mesh method is a very efficient variational approach to solve differential equations numerically with high precision. It has been used extensively in non-relativistic problems but in this work it has been formulated for the application to the Dirac equation to find bound states of fermions moving in deformed but axially symmetric potentials. This problem appears in the description of nuclei employing relativistic energy density functionals. The use of spherical coordinates leads to an eigenvalue problem with coupled equations of the radial wave functions. They are expanded in a set of Lagrange-mesh functions derived from Laguerre polynomials. Finally, a Hamiltonian matrix of moderate size has to be diagonalized.

The accuracy of the Lagrange-mesh method depends on three parameters: the scale factor, the number of grid points, and the number of channels. The dependence of the numerical precision on these parameters was studied by comparing analytical solutions for a deformed relativistic harmonic oscillator potential with numerical results obtained with the proposed formulation of the Lagrange-mesh method. With appropriate choices of the parameters, very accurate solutions are found without extensive computational efforts.

The description of ground states of self-conjugated even-even nuclei from ⁴He to ⁴⁰Ca was chosen as an application of the Lagrange-mesh method to a problem in nuclear structure physics. The equations of motion of nucleons, mesons, and the electromagnetic field were derived from a relativistic energy density functional based on a relativistic mean-field approach with density-dependent nucleon-meson couplings. The field equations were solved self-consistently by iteration, yielding the single-particle states of the nucleons and the intrinsic density distributions. A strong variation of the shapes and deformation parameters is found with spherical, oblate, and prolate distributions.

The present work focused on formulating the basic application of the Lagrange-mesh method to relativistic nuclear structure calculations. The presented approach can be extended easily in the future. The restriction to intrinsic states of given parity can be lifted allowing, e.g., also the description of octupole deformations. Proper projections, e.g., on given angular momentum and parity, are needed for a comparison of densities, radii, and energies with experimental data. Constraints on different multipole moments can be introduced and beyond mean-field calculations, e.g., using the generator-coordinate approach, are foreseeable.

Author Contributions

ST wrote and tested the computer program, performed all calculations, prepared the tables and figures, and wrote the manuscript.

Funding

This work was partly supported by the Deutsche Forschungsgemeinschaft (DFG) through grant no. SFB 1245.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Bender, M, Heenen, PH, and Reinhard, PG Self-consistent mean-field models for nuclear structure. Rev Mod Phys. (2003) 75:121–80. doi: 10.1103/RevModPhys.75.121

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Serot, BD, and Walecka, JD The relativistic nuclear many body problem. Adv Nucl Phys. (1986) 16:1–327.

PubMed Abstract | Google Scholar

3. Ring, P Relativistic mean field in finite nuclei. Prog Part Nucl Phys. (1996) 37:193–263. doi: 10.1016/0146-6410(96)00054-3

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Serot, BD, and Walecka, JD Recent progress in quantum hadrodynamics. Int J Mod Phys. (1997) E6:515–631. doi: 10.1142/S0218301397000299

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lalazissis, GA, Ring, P, and Vretenar, D Extended Density Functionals in Nuclear Structure Physics. Lecture Notes in Physics 641. Berlin; Heidelberg: Springer (2004).

PubMed Abstract | Google Scholar

6. Vretenar, D, Afanasjev, AV, Lalazissis, GA, and Ring, P Relativistic Hartree Bogoliubov theory: static and dynamic aspects of exotic nuclear structure. Phys Rept. (2005) 409:101–259. doi: 10.1016/j.physrep.2004.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Meng, J, Toki, H, Zhou, SG, Zhang, SQ, Long, WH, and Geng, LS Relativistic Continuum Hartree Bogoliubov theory for ground state properties of exotic nuclei. Prog Part Nucl Phys. (2006) 57:470–563. doi: 10.1016/j.ppnp.2005.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Nikšić, T, Vretenar, D, and Ring, P Relativistic nuclear energy density functionals: mean-field and beyond. Prog Part Nucl Phys. (2011) 66:519–48. doi: 10.1016/j.ppnp.2011.01.055

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ring, P, Abusara, H, Afanasjev, AV, Lalazissis, GA, Nikšić, T, and Vretenar, D Modern applications of covariant density functional theory. Int J Mod Phys. (2011) E20:253. doi: 10.1142/S0218301311017570

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Horowitz, CJ, and Serot, BD Selfconsistent Hartree description of finite nuclei in a relativistic quantum field theory. Nucl Phys. (1981) A368:503–28. doi: 10.1016/0375-9474(81)90770-3

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bonche, P, Flocard, H, Heenen, PH, Krieger, SJ, and Weiss, MS Self-consistent study of triaxial deformations: application to the isotopes of Kr, Sr, Zr and Mo. Nucl Phys. (1985) A443:39–63. doi: 10.1016/0375-9474(85)90320-3

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Reinhard, PG, Rufa, M, Maruhn, J, Greiner, W, and Friedrich, J Nuclear ground state properties in a relativistic meson field theory. Z Phys. (1986) A323:13–25. doi: 10.1007/BF01294551

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Pöschl, W, Vretenar, D, Rummel, A, and Ring, P Application of finite element methods in relativistic mean-field theory: spherical nucleus. Comp Phys Commun. (1997) 101:75–107. doi: 10.1016/S0010-4655(97)84583-3

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pöschl, W, Vretenar, D, and Ring, P Relativistic Hartree-Bogoliubov theory in coordinate space: finite element solution for a nuclear system with spherical symmetry. Comp Phys Commun. (1997) 103:217–50. doi: 10.1016/S0010-4655(97)00042-8

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Meng, J Relativistic continuum Hartree-Bogoliubov theory with both zero range and finite range Gogny force and their application. Nucl Phys. (1998) A635:3–42. doi: 10.1016/S0375-9474(98)00178-X

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Stoitsov, M, Ring, P, Vretenar, D, and Lalazissis, GA Solution of relativistic Hartree-Bogolyubov equations in configurational representation: spherical neutron halo nuclei. Phys Rev. (1998) C58:2086–91.

PubMed Abstract | Google Scholar

17. Typel, S, and Wolter, HH Relativistic mean field calculations with density dependent meson nucleon coupling. Nucl Phys. (1999) A656:331–64. doi: 10.1016/S0375-9474(99)00310-3

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Gambhir, YK, Ring, P, and Thimet, A Relativistic mean field theory for finite nuclei. Ann Phys. (1990) 198:132–79. doi: 10.1016/0003-4916(90)90330-Q

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ring, P, Gambhir, YK, and Lalazissis, GA Computer program for the relativistic mean field description of the ground state properties of even-even axially deformed nuclei. Comp Phys Commun. (1997) 105:77–97. doi: 10.1016/S0010-4655(97)00022-2

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhou, SG, Meng, J, and Ring, P Spherical relativistic Hartree theory in a Woods-Saxon basis. Phys Rev. (2003) C68:034323. doi: 10.1103/PhysRevC.68.034323

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Nikšić, T, Vretenar, D, Tian, Y, Ma, Zy, and Ring, P 3D Relativistic Hartree-Bogoliubov model with a separable pairing interaction. Phys Rev. (2010) C81:054318. doi: 10.1103/PhysRevC.81.054318

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Nikšić, T, Paar, N, Vretenar, D, and Ring, P DIRHB - A relativistic self-consistent mean-field framework for atomic nuclei. Comput Phys Commun. (2014) 185:1808–21. doi: 10.1016/j.cpc.2014.02.027

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Umar, AS, Strayer, MR, Wu, JS, Dean, DJ, and Guclu, MC Nuclear Hartree-Fock calculations with splines. Phys Rev. (1991) C44:2512–21. doi: 10.1103/PhysRevC.44.2512

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pei, JC, Stoitsov, MV, Fann, GI, Nazarewicz, W, Schunck, N, and Xu, FR Deformed coordinate-space Hartree-Fock-Bogoliubov approach to weakly bound nuclei and large deformations. Phys Rev. (2008) C78:064306. doi: 10.1103/PhysRevC.78.064306

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Pei, JC, Kortelainen, M, Zhang, YN, and Xu, FR Emergent soft monopole modes in weakly-bound deformed nuclei. Phys Rev. (2014) C90:051304. doi: 10.1103/PhysRevC.90.051304

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Davies, KTR, Flocard, H, Krieger, S, and Weiss, MS Application of the imaginary time step method to the solution of the static Hartree-Fock problem. Nucl Phys. (1980) A342:111–23. doi: 10.1016/0375-9474(80)90509-6

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Reinhard, PG, and Cusson, RY A comparative study of Hartree-Fock iteration techniques. Nucl Phys. (1982) A378:418–42. doi: 10.1016/0375-9474(82)90458-4

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ryssens, W, Hellemans, V, Bender, M, and Heenen, PH Solution of the SkyrmeHF+BCS equation on a 3D mesh, II: a new version of the Ev8 code. Comput Phys Commun. (2015) 187:175–94. doi: 10.1016/j.cpc.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Maruhn, JA, Reinhard, PG, Stevenson, PD, and Umar, AS The TDHF code Sky3D. Comput Phys Commun. (2014) 185:2195–216. doi: 10.1016/j.cpc.2014.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Schunck, N, McDonnell, JD, Sarich, J, Wild, SM, and Higdon, D Error analysis in nuclear density functional theory. J Phys. (2015) G42:034024. doi: 10.1088/0954-3899/42/3/034024

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Rodríguez, TR, Arzhanov, A, and Martínez-Pinedo, G Toward global beyond-mean-field calculations of nuclear masses and low-energy spectra. Phys Rev. (2015) C91:044315. doi: 10.1103/PhysRevC.91.044315

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Tanimura, Y, Hagino, K, and Liang, HZ 3D mesh calculations for covariant density functional theory. PTEP (2015) 2015:073D01. doi: 10.1093/ptep/ptv083

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Baye, D, and Heenen, PH Generalised meshes for quantum mechanical problems. J Phys A Math Gen. (1986) 19:2041. doi: 10.1088/0305-4470/19/11/013

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Baye, D, Hesse, M, and Vincke, M The unexplained accuracy of the Lagrange-mesh method. Phys Rev E (2002) 65:026701. doi: 10.1103/PhysRevE.65.026701

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Imagawa, H, and Hashimoto, Y Accurate random phase approximation calculation of low lying states on a three-dimensional Cartesian mesh. Phys Rev. (2003) C67:037302. doi: 10.1103/PhysRevC.67.037302

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hashimoto, Y Time-dependent Hartree-Fock-Bogoliubov calculations using a Lagrange mesh with the Gogny interaction. Phys Rev. (2013) C88:034307. doi: 10.1103/PhysRevC.88.034307

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Baye, D The Lagrange-mesh method. Phys Rep. (2015) 565:1–107. doi: 10.1016/j.physrep.2014.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Baye, D, Filippin, L, and Godefroid, M Accurate solution of the Dirac equation on Lagrange meshes. Phys Rev E (2014) 89:043305. doi: 10.1103/PhysRevE.89.043305

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Filippin, L, Godefroid, M, and Baye, D Relativistic polarizabilities with the Lagrange-mesh method. Phys Rev A (2014) 90:052520. doi: 10.1103/PhysRevA.90.052520

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Mohr, PJ, Newell, DB, and Taylor, BN CODATA recommended values of the fundamental physical constants: 2014. Rev Mod Phys. (2016) 88:035009. doi: 10.1103/RevModPhys.88.035009

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Typel, S Relativistic mean-field models with different parametrizations of density dependent couplings. Particles (2018) 1:1. doi: 10.3390/particles1010002

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Abramowitz, M, and Stegun, IA Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. New York, NY: Dover Publications, Inc. (1965).

PubMed Abstract | Google Scholar

43. Press, WH, Flannery, BP, Teukolsky, SA, and Vetterling, WT Numerical Recipes. Cambridge: Cambridge University Press (1987).

PubMed Abstract | Google Scholar

44. Lisboa, R, Malheiro, M, de, Castro AS, Alberto, P, and Fiolhais, M Pseudospin symmetry and the relativistic harmonic oscillator. Phys Rev. (2004) C69:024319. doi: 10.1103/PhysRevC.69.024319

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ginocchio, JN U(3) and pseudo-U(3) symmetry of the relativistic harmonic oscillator. Phys Rev Lett. (2005) 95:252501. doi: 10.1103/PhysRevLett.95.252501

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Guo, JY, Fang, XZ, and Xu, FX Pseudospin symmetry in the relativistic harmonic oscillator. Nucl Phys. (2005) A757:411–21. doi: 10.1016/j.nuclphysa.2005.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

47. de, Castro AS, Alberto, P, Lisboa, R, and Malheiro, M Relating pseudospin and spin symmetries through charge conjugation and chiral transformations: the case of the relativistic harmonic oscillator. Phys Rev. (2006) C73:054309. doi: 10.1103/PhysRevC.73.054309

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Ikeda, K, Takigawa, N, and Horiuchi, H The systematic structure-change into the molecule-like structures in the self-conjugate 4n nuclei. Prog Theor Phys Suppl. (1968) E68:464–75. doi: 10.1143/PTPS.E68.464

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Tohsaki, A, Horiuchi, H, Schuck, P, and Röpke, G Alpha cluster condensation in C-12 and O-16. Phys Rev Lett. (2001) 87:192501. doi: 10.1103/PhysRevLett.87.192501

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Girod, M, and Schuck, P α-particle clustering from expanding self-conjugate nuclei within the Hartree-Fock-Bogoliubov approach. Phys Rev Lett. (2013) 111:132503. doi: 10.1103/PhysRevLett.111.132503

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Typel, S, Röpke, G, Klähn, T, Blaschke, D, and Wolter, HH Composition and thermodynamics of nuclear matter with light clusters. Phys Rev. (2010) C81:015803. doi: 10.1103/PhysRevC.81.015803

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Dutra, M, Lourenço, O, Avancini, SS, Carlson, BV, Delfino, A, Menezes, DP, et al. Relativistic mean-field hadronic models under nuclear matter constraints. Phys Rev. (2014) C90:055203. doi: 10.1103/PhysRevC.90.055203

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wang, M, Audi, G, Kondev, FG, Huang, WJ, Naimi, S, and Xu, X The AME2016 atomic mass evaluation (II). Tables, graphs and references. Chin Phys C (2017) 41:030003. doi: 10.1088/1674-1137/41/3/030003

PubMed Abstract | CrossRef Full Text | Google Scholar