Minimization methods for the one-particle Dirac equation

Numerical methods avoiding the problem of variational collapse and the appearance of spurious roots are proposed for the computation of the eigenvalues and eigenstates of one-particle Dirac Hamiltonians. They are based on exact characterizations of the eigenvalues by a direct minimization of the corresponding Rayleigh quotient on a set described by a nonlinear constraint. PACS numbers : 31.30.Jv, 31.15.Pf, 03.65.Pm, 02.60.Cb The variational characterization and the numerical computation of the eigenstates of one-particle Dirac operators are difficult because of the presence of the negative continuum in the spectrum of the free Dirac operator. To deal with the related variational instability, various theoretical approaches based on squared Dirac operators [1,2], min-max formulations [3,4] or even more elaborated methods have been proposed [5], [6] as well as perturbative expansions of non-relativistic models and derivation of effective Hamiltonians. From a numerical viewpoint, the variational collapse and the existence of spurious states [7,8] are serious problems which have been dealt with by taking appropriate projections or imposing additional conditions in a large number of situations. In all those cases, a large quantity of precautions have been necessary to avoid the disagreements related to the total unboundedness of the Dirac operator. Our proposal consists in characterizing the eigenstates of the one-particle Dirac operator by stable and exact variational methods. Then, we will show that these methods can be easily translated into numerically tractable schemes for computing the eigenstates that ignore variational collapse. To illustrate the methods in a simple case, the details of the computations of the ground state are given below for a radial scalar potential V and numerical computations for V (x) = −γ/|x| with β ranging between 0 and 1. We consider the one-particle Dirac Hamiltonian H = α · (−i∇) +mcβ + V , (1) where α and β are the usual Dirac-Pauli matrices. The first natural idea to characterize the eigenvalues of H is based on min-max methods on the corresponding Rayleigh quotients (φ,Hφ)/(φ,φ) . (2) Indeed, the spectrum of the Dirac Hamiltonian H being unbounded from below, it is hopeless to just minimize the Rayleigh quotient without any further precaution, since this would take us to −∞. We refer to [9–13] for a mathematical study of min-max formulations. In [13] it was proved for a large class of potentials that if F+⊕F− is an orthogonal decomposition of a well-chosen space of smooth square integrable functions, the sequence of min-max levels λk = inf G subspace of F+ dim G=k sup φ∈(G⊕F−) φ6=0 (φ,Hφ) ||φ||2 (3) is equal to the sequence of the eigenvalues of H (counted with multiplicity) in the interval (−mc,+mc). An alternative and equivalent approach to min-max characterizations was introduced in [11]. Composing H with its associated positive energy projector Λ would allow to minimize Rayleigh quotients but this idea is formal since the above projector is a priori unknown. Our main idea is to find a (nonlinear) constraint which implicitely defines the projector. To do that, the Dirac equation Hψ = λψ, originally written for 4-spinors, is first reduced to an equivalent equation for 2-spinors. This is usually done to recover Schrödinger-type equations in the non-relativistic limit: the lower 2-spinor, or small component is written in terms of the large component by assuming that E = λ − mc and V are neglegible with respect to mc (see [14] for more details on kinetic balance). Keeping track of E and V can also be done: see [15,14]. More precisely, for any ψ with values in CI , let us write ψ = ( χ ), with φ, χ taking values in CI . Then the equation Hψ = λψ is equivalent to the system

The variational characterization and the numerical computation of the eigenstates of one-particle Dirac operators are difficult because of the presence of the negative continuum in the spectrum of the free Dirac operator.To deal with the related variational instability, various theoretical approaches based on squared Dirac operators [1,2], min-max formulations [3,4] or even more elaborated methods have been proposed [5], [6] as well as perturbative expansions of non-relativistic models and derivation of effective Hamiltonians.From a numerical viewpoint, the variational collapse and the existence of spurious states [7,8] are serious problems which have been dealt with by taking appropriate projections or imposing additional conditions in a large number of situations.In all those cases, a large quantity of precautions have been necessary to avoid the disagreements related to the total unboundedness of the Dirac operator.Our proposal consists in characterizing the eigenstates of the one-particle Dirac operator by stable and exact variational methods.Then, we will show that these methods can be easily translated into numerically tractable schemes for computing the eigenstates that ignore variational collapse.To illustrate the methods in a simple case, the details of the computations of the ground state are given below for a radial scalar potential V and numerical computations for V (x) = −γ/|x| β with β ranging between 0 and 1.
We consider the one-particle Dirac Hamiltonian where α and β are the usual Dirac-Pauli matrices.The first natural idea to characterize the eigenvalues of H is based on min-max methods on the corresponding Rayleigh quotients (ϕ, Hϕ)/(ϕ, ϕ) . (2) Indeed, the spectrum of the Dirac Hamiltonian H being unbounded from below, it is hopeless to just minimize the Rayleigh quotient without any further precaution, since this would take us to −∞.We refer to [9][10][11][12][13] for a mathematical study of min-max formulations.In [13] it was proved for a large class of potentials that if F+⊕F− is an orthogonal decomposition of a well-chosen space of smooth square integrable functions, the sequence of min-max levels is equal to the sequence of the eigenvalues of H (counted with multiplicity) in the interval (−mc 2 , +mc 2 ).
An alternative and equivalent approach to min-max characterizations was introduced in [11].Composing H with its associated positive energy projector Λ + would allow to minimize Rayleigh quotients but this idea is formal since the above projector is a priori unknown.Our main idea is to find a (nonlinear) constraint which implicitely defines the projector.To do that, the Dirac equation Hψ = λψ, originally written for 4-spinors, is first reduced to an equivalent equation for 2-spinors.This is usually done to recover Schrödinger-type equations in the non-relativistic limit: the lower 2-spinor, or small component is written in terms of the large component by assuming that E = λ − mc 2 and V are neglegible with respect to mc 2 (see [14] for more details on kinetic balance).Keeping track of E and V can also be done: see [15,14].More precisely, for any ψ with values in C I 4 , let us write ψ = ( ϕ χ ), with ϕ, χ taking values in C I 2 .Then the equation Hψ = λψ is equivalent to the system with L = ic( σ. ∇) = 3 k=1 icσ k ∂/∂x k , σ k being the Pauli matrices (k = 1, 2, 3).As long as λ+mc 2 −V = 0, the system (4) can be written as where In a second step, we consider φ = ϕ/ √ g E which solves where ∆(E, φ) . Heuristically (7) is one of the analogues in quantum mechanics of Einstein's relation: It turns out that the critical points of J ± (E, φ) under the constraint E = J ± (E, φ) are exactly the eigenfunctions of the Dirac operator H corresponding to eigenvalues in the gap, and the eigenvalues corresponding to J + converge as c → +∞ to the nonrelativistic eigenvalues of the Schrödinger operator, at least at a formal level.Note that if V is nonpositive, the range of J − is contained in (−∞, −mc 2 ] which corresponds to the negative part of the spectrum of the Dirac operator. In any case, we may define the ground-state as the lowest positive eigenvalue, which turns out to be the minimum of J + (E, φ) under the constraint E = J + (E, φ) (and converges as c → +∞ to the ground-state of the Schrödinger operator) whatever the values of V are, and is compatible with the nonrelativistic limit (but may not correspond to the lowest eigenvalue in the gap).From now on we assume that m = 1 and c = 1.We may notice that for the problem to be well defined, further conditions have to be assumed on the potential.If for instance V is the Coulomb potential, γ = Ze 2 hc has to be less than 1 (which means Z < 1/α = 137.036...) unless H is not well defined as a self-adjoint operator any more: see [16,11,13] for more details.The fact that the whole spectrum in the gap is actually characterized by J ± comes from the equivalence of the method with the min-max characterization of the eigenvalues in the gap [13].
Since V is radial (see for instance [16,17]) the eigenfunctions can be expressed in terms of the spherical harmonics: any square integrable spinor ψ defined on IR 3 with values in C I 4 can therefore be written as where j = 1 2 , 3 2 , ..., mj = −j, −j + 1, ... + j, κj = ±(j being the usual spherical harmonics.The radial Dirac operator acting on the space L 2 (0,+∞) of the square integrable real functions on (0, +∞) is With c = m = 1, the eigenvalue problem takes the form The solutions of this system are characterized by two parameters, λ and δ = v(1)/u(1) for instance, and we shall denote by X the set of the solutions of ( 13) such that u(1) = 1 when λ and δ vary in IR.However, the condition that u and v are in L 2 (0, +∞) determines uniquely λ and δ.One can show that this condition is equivalent to assuming that thus providing us with a first (and simple) numerical ("shooting") method to determine λ and δ.We shall refer to this method by the letter "s" and use it to compare the numerical results with the numerical minimization method given below.The approximated eigenvalues computed with this method will be denoted by λ s in Table I below.
Let us describe now the minimization method based on J + .Similarly to (5), v can be eliminated in terms of u: We are looking for the ground-state so we may choose κ = −1 and the eigenvalue problem ( 13) is now equivalent to solve where h λ is a formally self-adjoint operator: and φ(r) = r −1 u(r)/ √ 1 + λ − V is now a function defined on (0, +∞).Equation ( 6) is then equivalent to where (., .) is the usual scalar product in L 2 (0, +∞) and .the corresponding norm.The problem is then reduced to finding a critical point of J + (λ−1,.)with λ − 1 = J + (λ−1,φ) and To solve this constrained problem numerically, the natural idea is to introduce a penalization method and to minimize J + (λ−1,φ) + A|(λ−1) − J + (λ−1,φ)| 2 in the limit A → +∞.Actually if we assume that φ is given by ( 17) with (u, v) in X, the condition that u and v are in L 2 (0, +∞) is equivalent to assuming that the integrals involved in the expression (19) are finite.Of course these integrals are numerically computed on an interval (ǫ, R) and the approximate value J + ǫ,R of J + is finite even if the constraint is not satisfied, but we observe that lim (ǫ,R)→(0,+∞) J + ǫ,R (λ−1,φ) = +∞ unless λ − 1 = J + (λ−1,φ).A minimization of J + (numerically of J + ǫ,R ) on the set X without constraint is therefore equivalent to a constrained minimization of J + .This method will be referred by the letter "m" in Table I, which contains the results of our computations using both the shooting and the minimization methods described above.The set of functions over which the approximated eigenvalues are computed consists of all the solutions of ( 13), with δ and λ to be determined.
The main advantage of the minimizing setting described above is that it can be extended to nonsymmetric situations (non central potentials), but of course for a minimizing set which is larger than X.Indeed, above we have minimized J + only among the set of functions which are already solutions of a radially symmetric system (13).In the rest of this letter, for convenience, we still assume that the potential is radial, but consider a general basis of L 2 (0, +∞) (of course well chosen).Next, we introduce a third formulation, which is intermediate between the abstract min-max theory and the minimization of J + described above.Its main advantage is that the constraint E = J + (E, φ) will then be automatically satisfied.We will therefore call this method the "direct minimization method".
Multiplying Equation ( 22) by ϕ and integrating with respect to x ∈ IR 3 , we get :

TABLE I.
Comparison of the shooting (s) and the minimization (m) methods for κ = −1, m = c = 1, V (r) = −γr −β , γ = 0.5 and β ∈ (0, 1).The system ( 13) is numerically solved with a stepsize adaptative Runge-Kutta method on the interval (ǫ = 10 −4 , R = 15).For the shooting method we minimize the quantity ǫ(|u(ǫ ∆s for some scale parameter θ > 0 (which is chosen to balance both boundary terms), while for the minimization method, the quantity J + (λ−1,φ) is directly minimized, the quantity −0.267949....For practical reasons, the results given here correspond to parameters taken in a neighborhood of (δ1, λ1).The results correspond therefore to the branch (δ β , λ β ) starting from (δ1, λ1) at β = 1 and parametrized by β.Since for a given ϕ, fϕ(λ) is decreasing and gϕ(λ) is increasing, if there exits a λ = λ[ϕ] such that (23) is satisfied, then it is unique (the existence of such a λ for all u depends on the properties of the potential V ).According to [13], for those V 's, the ground state is such that For a radial potential we may use the radial Dirac equation and consider (13) instead of (4).For m = 1 and c = 1, λ = λr[u] is then the unique solution of To solve it numerically, it is more convenient to rewrite f (λ) as Numerically, the solution (with κ = −1) is approximated on a finite basis (ui)i=1,2,...n: u = n i=1 xiui.If and the approximating equation for λ corresponding to (23) is then where the series in λ has been truncated at order m.It is actually more convenient to define and to approximate λ1 by λ n,m 1 defined as the first positive root of λ → µ(A n,m (λ)) where µ(A) denotes the first eigenvalue of the matrix A (see [13] for more details).Note that (λ n,m 1 ) m≥1 is an alternating sequence (which essentially converges at a geometric rate): two consecutive eigenvalues determine an interval containing limm→+∞ λ n,m 1 .The results in Table II have been obtained by taking an orthonormal basis generated by the ground state of the hydrogen atom and n − 1 Hermite functions, with n = 10 and then minimizing λr in the space generated by that basis, without taking any further precaution.Our purpose in this numerical computation is not to provide very accurate results but just to prove the feasibility of such a numerical approach.Clearly, depending on the specific properties of the potential, the choice of a well suited basis should greatly improve the accuracy of the computation.
In this letter we have described three different ways of computing the ground state of the one-particle Dirac operator: a first one is based on a "shooting" argument (specific to the case of central external potentials) and serves us as test for the variational computations.The second method is based on the minimization (in a set defined by means of a nonlinear constraint) of a functional which is bounded from below.Hence, no particular precaution has to be taken to avoid variational collapse, the method is automatically free of that "disease".The third method consists in the minimization (in the space of all upper spinors) of a function which is again bounded from below, and which is defined by an explicit identity.These approaches, which can be applied to general external potentials, give rise to numerical schemes ("minimization" in Table I, "direct minimization" in Table II

1 and λ 10 1 )
), which are based on simpleTABLE II.Direct minimisation method for κ = −1, m = c = 1, V (r) = −γr −β , γ = 0.5 and β close to 1.The approximating space is of dimension n = 10 and the series are truncated at m = 14 or m = 15 (the corresponding values λ 10,14 and an upper bound of limm→+∞ λ 10,m .As in Table I, J + is obtained through a minimization on the set X, and ∆ m := 1 − 10 i=1 (u m , ui) 2 1/2measures the error (in the L 2 -norm) when the corresponding solution is approximated on the basis (with n = 10 elements) used for the direct minimisation method. β the exact eigenvalues and eigenstates of Dirac operators.So, if implemented with large enough basis sets, they can yield the eigenvalues and the eigenfunctions to any desired degree of approximation.