hp-Cloud approximation of the Dirac eigenvalue problem: The way of stability

Abstract

We apply hp-cloud method to the radial Dirac eigenvalue problem. The difficulty of occurrence of spurious eigenvalues among the genuine ones in the computation is resolved. The method of treatment is based on assuming hp-cloud Petrov–Galerkin scheme to construct the weak formulation of the problem which adds a consistent diffusivity to the variational formulation. The size of the artificially added diffusion term is controlled by a stability parameter (τ). The derivation of τ assumes the limit behavior of the eigenvalues at infinity. The parameter τ is applicable for generic basis functions. This is combined with the choice of appropriate intrinsic enrichments in the construction of the cloud shape functions.

Introduction

In the last decades, several numerical methods have been derived to compute the eigenvalues of operators. The need of accurate computation of eigenvalues is considered due to their significant applications in many disciplines of science. Mathematically, if a matrix or a linear operator is diagonalized, then by the spectral theorem, it can be analyzed by studying its corresponding eigenvalues, i.e., transforming the matrix or operator to a set of eigenfunctions which can be easily studied. From the physical point of view, the eigenvalues possess a wide range of information about the behavior of the system governed by an operator. This information might be all what is needed to answer many questions regarding the system properties such as stability, positivity, boundedness, asymptotic behavior, etc. In mechanics, eigenvalues play a central role in several aspects such as determining whether the automobile is noisy, whether a bridge will collapse by the water waves, etc. Also, the eigenvalues describe how the quantum state of a physical system changes in time (Schrödinger equation). They also represent the electrons relativistic energies and describe their motion in the atomic levels, this is the well-known Dirac equation, which is the core of the present work.

The calculation of energy levels in Helium-like ions, where the interaction between two electrons takes place, can be determined by studying the electrons correlation which is part of quantum electrodynamic effects (QED-effects). A scheme for calculating QED-effects [31], [35], [40], [42] is based on a basis set of relativistic Hydrogen-like ion wave eigenfunctions (of the Dirac operator). Meanwhile, the numerical computation of the Dirac operator eigenvalues encounters unphysical values (do not match the physical observations) called spurious eigenvalues or spectrum pollution. The spurious eigenvalues result in rapid oscillations in the wave functions, hence, in many cases, affecting the computation reliability of the basis set in the practical atomic calculations.

The spurious eigenvalues are an effect of the numerical methods and are found in the computation of many problems other than the Dirac eigenvalue problem [1], [2], [39], [43]. For general eigenvalue problems, spurious eigenvalues are reported in [48]. The occurrence of the spuriosity is related to mismatching of desired properties of the original solution in the numerical formulation. Also in the computation of electromagnetic problems the spuriosity is perceived [36], [41]. Two leading approaches are derived to solve this difficulty; Shabaev et al. [43] have related the spuriosity to the symmetric treatment of the large and small components of the Dirac wave function. Their approach, for removing the spurious eigenvalues, is based on deriving dual kinetic-balance (DKB) basis functions for the large and small components. Almanasreh et al. [2] have allied the occurrence of spurious eigenvalues to the incorrect treatment of the trial and test functions in the finite element method (FEM). They proposed a stability scheme based on creating diffusivity by modifying the test function so that it includes a gradient-based correction term.

In this work, we apply hp-cloud method [15], [49] to the radial Dirac eigenvalue problem. The technique is known as one of the meshfree methods (MMs) [6], [19], [32], [33], [37] that allows two different enrichments, intrinsic and extrinsic, to be built in the construction of the basis functions. The method was previously applied for several problems, e.g., the Schrödinger equation [10], Mindlin's thick plate model [20], and Timoshenko beam problems [34], etc. Here, we apply hp-cloud method based on the Galerkin formulation. This means that it is required to evaluate the integrals in the weak formulation of the particular equation, thus a background mesh must be employed in the integration techniques. Therefore, the hp-cloud method used here is not really a truly MM. However, all other features of MMs are maintained in our approximation.

In order to treat the spuriosity problem, we stabilize the computation by considering instead an hp-cloud Petrov–Galerkin (hp-CPG) method which is a technique of the general meshfree local Petrov–Galerkin (MLPG) [4], [18], [30]. The stability scheme is based on adding consistent diffusion terms without changing the structure of the equation. The size of the additional diffusivity is controlled by a stability parameter.

Consider the radial Dirac eigenvalue problem $H_{\kappa }\Phi (x)=\lambda \Phi (x)$, where $\Phi (x)={(F(x),G(x))}^t$ is the radial wave function, λ is the electron relativistic energies (eigenvalues), and $H_{\kappa }$ is the radial Dirac operator given by

$$H_{\kappa }=\left(\begin{array}{cc}\hfill m{c}^2+V(x)\hfill & \hfill c(-D_x+\frac{\kappa }{x})\hfill \\ \hfill c(D_x+\frac{\kappa }{x})\hfill & \hfill -m{c}^2+V(x)\hfill \end{array}\right).$$

The constant c is the speed of light, m is the electron mass, V is the Coulomb potential, $D_x=d/dx$, and κ is the spin–orbit coupling parameter defined as $\kappa ={(-1)}^{+\ell +\frac{1}{2}}(+\frac{1}{2})$, where ȷ and ℓ are the total and the orbital angular momentum quantum numbers respectively. The weak formulation of the problem is to find $\lambda \in \mathbb{R}$ and Φ in a specific function space such that for every test function Ψ in some suitable space we have ${\int }_{\Omega }{\Psi }^t{H}_{\kappa }\Phi dx=\lambda {\int }_{\Omega }{\Psi }^t\Phi dx$. The usual hp-cloud Galerkin approximation is to let Ψ to be ${(\psi ,0)}^t$ and ${(0,\psi )}^t$, where ψ is in the same space as of the two components of Φ. To discretize the weak form, the components of the trial function Φ and the test function ψ are chosen from a finite set of hp-cloud basis functions which are constructed using moving least-squares method. Since the radial Dirac operator is dominated by advection (convection) terms, the hp-cloud approximation will be upset by spurious eigenvalues.

To stabilize the hp-cloud approximation, the hp-CPG method is used. In this formulation, the test function Ψ is assumed to belong to a function space different from that of the trial function Φ, in the sense that Ψ is chosen in the form ${(\psi ,\tau {\psi }^{\prime })}^t$ and ${(\tau {\psi }^{\prime },\psi )}^t$ where ψ belongs to the same space as the two components of Φ. The correction term $\tau {\psi }^{\prime }$ is used to add artificial viscosity, controlled by τ, to stabilize the convection terms. The derivation of the stability parameter τ follows the principle used in [2] for the FEM. Two assumptions are considered in deriving τ; (i) the operator limit as the radial variable x tends to infinity, thus obtaining an approximation of the limit point of the eigenvalues (depending on τ) which can be compared to the theoretical limit point eigenvalue [21], (ii) considering the dominant terms with respect to the speed of light (c).

The paper is organized as follows; in Section 2, some preliminaries about the Dirac equation are presented, also we shed some light over the occurrence of the spuriosity. In Section 3, the construction of the hp-cloud functions is provided, also coupling with the FEM to impose essential boundary conditions (EBCs) is explained. The hp-CPG method and the derivation of the stability parameter are treated in Section 4. In the last section, Section 5, we present some numerical results and provide necessary discussion about the stability scheme.