\documentstyle[twocolumn,aps]{revtex}
\newcommand{\Frac}{\displaystyle \frac}
\newcommand{\Sup}{\displaystyle \sup}
\def\BbR{{\rm I\!R}}
\def\BbC{{\rm C\!\!\!I}}
\def\BbZ{{\rm Z\!\!\!\sl I}}
\def\BbN{{I\!\!N}}
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Front matter %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Minimization methods for the one-particle Dirac equation}
\author{Jean Dolbeault, Maria J. Esteban, Eric S{\'e}r{\'e} and Michel
Vanbreugel}
\address{CEREMADE (UMR C.N.R.S. 7534)\\
Universit{\'e} Paris IX-Dauphine\\
Place du Mar{\'e}chal de Lattre de Tassigny\\
75775 Paris Cedex 16 - France\\
E-mail: dolbeaul, esteban or sere@ceremade.dauphine.fr,
root@ceremade.dauphine.fr}
\date{\today}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent\abstract{{\bf Abstract. }{\sl Taking into account relativistic
effects in quantum chemistry is crucial for accurate computations
involving heavy atoms. Standard numerical methods can deal with
the problem of {\rm variational collapse} and the appearance of {\rm
spurious roots} only in special cases. The goal of this letter is to
provide a general and robust method to compute particle bound states of
the Dirac equation.}}
\vspace{2mm}
\noindent PACS numbers : 31.30.Jv, 31.15.Pf, 03.65.Pm, 02.60.Cb}
\vspace{4mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The free Dirac operator has a negative continuum in its spectrum. This is
the cause for a variational instability which makes the numerical
computations of one-particle bound states of Dirac equations difficult.
Various approaches based for instance on
squared Dirac operators \cite{[Wal-Ku],[Baylis-Peel]}, min-max formulations
\cite{[Ta],[Da-De]}, special basis \cite{[Dr-Go81],[Ku97]} or even more
elaborated methods have been proposed
\cite{[Hill-Kra],[Dr-Go]}, as well as perturbative expansions of
non-relativistic
models and derivation of effective Hamiltonians. None of these remedies
provides a
complete and satisfactory answer. From a numerical viewpoint, the {\sl
variational
collapse\/} and the existence of {\sl spurious states}
\cite{[Ku84],[Grant]} are
serious problems which have been solved in special cases by taking appropriate
projections or imposing additional conditions, for instance on the
boundary \cite{[Dr-Go81]}.
Here we propose exact and stable variational meth\-ods. First we present
an abstract min-max approach and show how to reduce it to a constrained
minimization problem. The method requires no specific precaution and
allows for instance the presence of several nuclei. It is then
illustrated by three computations of the {\sl ground state\/} for a
particular class of spherical potentials: a {\bf shooting method}
(specific to the case of spherical potentials), which is compared with a
{\bf minimization method} under constraint (see Table I), and a {\bf
direct minimization method} based on a decomposition on a general basis
(see Table II) -- which is equivalent to the minimization method -- but
has a form which is closer to the min-max formulation and is easier to
generalize to non central potentials.
\vspace{2mm}
To compute the bound states of the one-particle Dirac Hamiltonian
\begin{equation}\label{Hamiltonian}
H=\alpha\cdot(-i\nabla)+mc^2\beta+ V\; ,
\end{equation}
where $\alpha$ and $\beta$ are the usual Dirac-Pauli matrices, it is
natural to characterize the energy levels of $H$ by min-max methods applied to
the Rayleigh quotients $(\varphi, H\varphi)/(\varphi,\varphi)$.
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 $-\infty$.
We refer to
\cite{[Es-Se],[Gri-Sie],[Dol-Est-Ser98],[Gri-Lew-Sie],[Dol-Est-Ser99]} for a
mathematical study of min-max formulations.
In \cite{[Dol-Est-Ser99]} it was proved for a large class of potentials that if
$F_+{\displaystyle \oplus}F_-$ is an orthogonal decomposition of a well-chosen
space of smooth square integrable functions, the sequence of
min-max levels \begin{equation}\label{min-max}
\lambda_k =\ \inf_{\scriptstyle G\ {\rm subspace \ of \ } F_+\atop
\scriptstyle{\rm dim}\ G=k}\ \Sup_{\scriptstyle \varphi\in (G\oplus F_-)
\atop \varphi\ne 0 }\ \Frac{(\varphi, H\varphi)}{||\varphi||^2}
\end{equation}
is equal to the sequence of the energy levels of $H$ (counted with
multiplicity) in the
interval $(-mc^2,+mc^2)$.
\vspace{2mm}
This expression has however not much pratical interest and we will now
reformulate it as a minimization problem (see \cite{[Dol-Est-Ser98]}
for mathematical justifications). Composing $H$ with its associated
positive energy projector $\Lambda^+$ would allow to minimize Rayleigh
quotients but this idea is formal since the above projector is {\sl
a priori\/} unknown. Instead we will introduce a (nonlinear) constraint
which implicitely defines the projector. The Dirac equation
$H\psi=\lambda\psi$, originally written for $4$-spinors, is first
reduced to an equivalent equation for $2$-spinors: the lower
$2$-spinor, or
{\sl small component\/} is written in terms of the {\sl large
component\/} (see \cite{[LLBS],[Ku97]}). More precisely, for any $\psi$
with values in $
\BbC^{\,4}$, let us
write
$\psi=({\varphi\atop\chi})$, with $\varphi,\chi$ taking values in
$\BbC^{\,2}$. Then the equation $H\psi=\lambda\psi$ is equivalent to the
system
\begin{equation}\label{System2spinors}
\left\lbrace \begin{array}{l}L\chi \ =\ (\lambda-mc^2- V)\
\varphi\\ L\varphi\ =\ (\lambda+mc^2- V)\ \chi\\ \end{array}\right.
\end{equation}
with $L\! =\! ic(\vec{\sigma}.\vec{\nabla})\! =\!
\sum_{k=1}^3 ic\sigma_k\partial/\partial x_k$, $\sigma_k$ being the
Pauli matrices ($k=1,2,3$). As long as
${\lambda\! +\! mc^2\! -\! V\! \neq\! 0}$,
the system (\ref{System2spinors}) can be written as
\begin{equation}\label{elim}
L\Bigl(\frac{L\varphi}{g_{_E}}\Bigr)+ V\varphi\ =\ E\varphi\
,\quad\chi=
\frac{L\varphi}{g_{_E}}
\end{equation}
where $g_{_E}=E+2mc^2-V$, $E=\lambda-mc^2$.
Then, system (\ref{System2spinors}) written for
$\phi=\varphi/\sqrt{g_{_E}}$ reads
\begin{equation}\label{Eqn2spinor}
H_E\phi = {\scriptstyle(E-V)(E+2mc^2-V)}\phi\end{equation}
where the operator $H_E\phi:=\sqrt{g_{_E}}
\ L(\frac{1}{g_{_E}}L(\sqrt{g_{_E}}\ \phi))$ is symmetric
(and self-adjoint for an appropriate domain). Thus $E$ is a
solution of
$(\phi,\phi)\, E^2+2(\phi,(mc^2-V)\phi)
E-(\phi,(2mc^2-V)V\phi)-(\phi,H_E\phi)=0$,
\begin{equation}\label{Jpm}
E=J^\pm(E,\phi):={\scriptstyle\frac{1}{(\phi,\phi)}\Bigl(\pm\
\sqrt{\Delta(E,\phi)}-
(\phi,(mc^2-V)\phi)\Bigr)}
\end{equation}
where
$\Delta(E,\phi)\! :=\! |(\phi,V\phi)|^{^2}\! +\! (\phi,\phi)
[m^2c^4(\phi,\phi) \! +\! (\phi,H_E\phi)\! -\! (\phi,V^2\phi)]$.
Heuristically (\ref{Jpm}) is one of the analogues
in quantum mechanics of Einstein's relation: $E=\pm
(p^2c^2+m^2c^4)^{1/2}- mc^2+V$.
It turns out that the critical points of
$J^\pm(E,\phi)$ under the constraint $E=J^\pm(E,\phi)$ are exactly (see
\cite{[Dol-Est-Ser98]}) the bound states of the Dirac operator $H$
corresponding to energies in the gap (and those corresponding to $J^+$
converge as $c\to +\infty$ to the nonrelativistic energy levels of the
Schr\"odinger operator). Note that if
$V$ is nonpositive, the range of $J^-$ is contained in
$(-\infty,-mc^2]$ which corresponds to the negative part of the
continuum of the spectrum of the Dirac operator.
We may in that case define the {\sl ground-state\/} as
the lowest energy level in the gap, which turns out to be the minimum of
$J^+(E,\phi)$ under the constraint
$E=J^+(E,\phi)$.
The functional $\,J^+\,$ is well defined only if further conditions
are assumed on the potential. If for instance $V=\frac\zeta{|x|}$ is the
Coulomb potential,
$\zeta=\frac{Ze^2}{\hbar c}$ has to be less than $1$ (which means
$Z<1/\alpha= 137.036...$), otherwise $H$ is not well defined as a
self-adjoint operator: see \cite{[Th],[Dol-Est-Ser98],[Dol-Est-Ser99]} for
more details. The fact that the whole spectrum in the gap is actually
characterized by $J^\pm$ comes from the equivalence of the method with
the min-max characterization of the energy levels in the gap
\cite{[Dol-Est-Ser99]}.
\vspace{2mm}
If $V$ is spherically symmetric (see for instance \cite{[Th],[Bj-Dr]})
the bound states can be expressed in terms of the spherical
harmonics using partial wave Dirac operators acting on the space
$(L^2{\scriptstyle (0,+\infty)})^2$ of the square integrable real
functions on $(0,+\infty)$, which have the form
\begin{equation}\label{radialDirac}
h=\left({\begin{array}{cc} mc^2+V &-c\frac d{dr}+\frac{c\kappa}{r}\\
c\frac
d{dr}+\frac{c\kappa}{r}&-mc^2+V\end{array}}\right)\quad(\kappa=\pm 1,\;
\pm 2,...)\end{equation} From now on we choose a system of units in which
$m=1$ and $c=1$. Problem (\ref{System2spinors})
takes the form
\begin{equation}\label{systemradial}
\left\{\begin{array}{c} u'=(1+\lambda)v-(Vv+\frac\kappa{r}u)\\
v'= (1-\lambda)u+Vu+\frac\kappa{r}v)\end{array}\right.
\end{equation}
The solutions of this system are characterized by two parameters,
$\lambda$ and $\delta=v(1)/u(1)$ for instance, and we shall denote by
$X$ the set of the solutions of (\ref{systemradial}) such that
$u(1)=1$ when $\lambda$ and $\delta$ vary. However, the
condition that
$u$ and $v$ are in $L^2(0,+\infty)$ determines $\lambda$ and
$\delta$. One can show that this condition is equivalent to assuming
that
\begin{equation}\label{shootingCondition}
{\lim_{r\to 0_+}r(|u(r)|^2+|v(r)|^2)=0\; ,\atop\lim_{r\to
+\infty}(|u(r)|^2+|v(r)|^2)=0\; ,}
\end{equation}
thus providing us with a first numerical {\bf shooting
meth\-od} to determine $\lambda$ and $\delta$ (this is of course
valid only for spherically symmetric potentials).
We shall refer to this method by the letter {\bf s}
and use it as a comparison test for the numerical results obtained
by the minimization approach given below. The approximated energy levels
computed with this method will therefore be denoted by $\,\lambda^{\rm
s}\,$ in Table I.
Let us describe now the {\bf minimization me\-thod} based on $J^+$.
Similarly to
(\ref{elim}), $v$ can be eliminated in terms of $u$:
\begin{equation}\label{elimradial}
\frac v{r^\kappa}=(r^{2\kappa}(1+\lambda-V))^{-1}\frac
d{dr}({r^\kappa}u)\; .
\end{equation}
We are looking for the ground-state
so we may choose $\kappa=-1$ and problem
(\ref{systemradial}) is now equivalent to solving
\begin{equation}\label{equationsemiclassic}
h_\lambda\phi=(1+\lambda-V)(1-\lambda+V)\phi
\end{equation}
where $h_\lambda$ is formally a self-adjoint operator:
\begin{equation}\label{huv}
h_\lambda\phi=\sqrt{1+\lambda-V}\frac
d{dr}\biggl[\frac{r^2}{1+\lambda-V} \frac
d{dr} (\sqrt{1+\lambda-V}\phi)\biggr]
\end{equation}
and $\phi(r)=r^{-1} u(r)/\sqrt{1+\lambda-V}$ is now a function
defined on $(0,+\infty)$. Equation (\ref{Eqn2spinor}) is then equivalent
to
\begin{equation}\label{EqnJradial}
(\phi,\phi)\lambda^2-2(\phi,V\phi)\lambda
+(\phi,V^2\phi)-[\Vert\phi\Vert^2+(\phi,h_\lambda\phi)]=0 \end{equation}
where $(.,.)$ is the usual scalar product in
$L^2(\! 0,+\infty\! )$ and $\Vert.\Vert$ the corresponding norm. The
problem is then reduced to finding a critical point of
$J^+(\scriptstyle{\lambda-1},.)$
with $\lambda-1=J^+({\scriptstyle \lambda-1,\phi})$ and
\begin{equation}\label{J+radial}
J^+(\lambda\! -\! 1,\phi)+1=
\frac{\sqrt{\Delta(\lambda-1,\phi)}+(\phi,V\phi)} {\Vert\phi \Vert^2}\;
,\end{equation} \begin{equation}\label{discriminantradial}
\Delta(\lambda\! -\! 1,\phi)\! =\!(\phi,V\phi)^2
\! +\! \Vert\phi \Vert^2[(\phi,h_\lambda\phi)\! +\! \Vert\phi \Vert^2\!
-\! (\phi,V^2\phi)]\; .\end{equation}
To solve this constrained problem numerically, the
natural idea is to introduce a penalization method and to minimize
$J^+({\scriptstyle \lambda-1,\phi})
+A|{\scriptstyle(\lambda-1)}-J^+({\scriptstyle
\lambda-1,\phi})|^2$ for
$A$ large enough. Actually if we assume that $\phi$ is given by (\ref{huv})
with $(u,v)$ in $X$, the condition that $u$ and $v$ are in
$L^2(0,+\infty)$ is equivalent to assuming that the integrals involved in
(\ref{J+radial}) are finite. Of course these integrals are numerically
computed on an interval $(\epsilon, R)$ and the approximate value
$J^+_{\epsilon, R}$ of $J^+$ is finite even if the constraint is not
satisfied, but we observe that $\lim_{(\epsilon, R)\to
(0,+\infty)}J^+_{\epsilon, R}({\scriptstyle
\lambda-1,\phi})=+\infty$ unless
$\lambda-1=J^+({\scriptstyle \lambda-1,\phi})$. A minimization of
$J^+$ (numerically,
of $J^+_{\epsilon, R}$) on the set $X$ without constraint is
therefore equivalent to
a constrained minimization of $J^+$. This method will be referred by
the letter {\bf m}
in Table \ref{table1}, 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 energy levels are computed consists of all the
solutions of
(\ref{systemradial}), with $\delta$ and $\lambda$ to be determined. Once
this is done, as a byproduct of the above minimization method, the wave
functions associated with the computed energy levels are given by
(\ref{systemradial}). This is also true for the shooting method.
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 (\ref{systemradial}). 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,+\infty)$ and 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,\phi)$ will be automatically satisfied. We will therefore call
this method the {\bf direct minimization method}.
As in Equation (\ref{elim}) we may rewrite
(\ref{System2spinors}) as
\begin{equation}\label{elimchi}
\chi=(\lambda+mc^2-V)^{-1}L\varphi\; ,\end{equation}
\begin{equation}\label{eqnphi}
L\biggl(\frac{L\varphi}{\lambda+mc^2-V}\biggr)=(\lambda-mc^2-V)\varphi\; ,
\end{equation}
at least for any $\lambda\in ]-mc^2, +mc^2[$ if $V$ takes negative
values. Multiplying Equation (\ref{eqnphi}) by $\varphi$ and integrating
with respect to $x\in \BbR^3$, we get~:
\begin{eqnarray}\label{eqnlambda}
f_\varphi(\lambda):=&&\int \frac{|L\varphi|^2}{\lambda+mc^2-V}\; dx \\=&&
(\lambda-mc^2)\Vert\varphi\Vert^2_{L^2}-\int V|\varphi|^2\; dx =:
g_\varphi(\lambda)\; .\nonumber\end{eqnarray}
\vspace{5mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\widetext
\begin{table}
\caption{{Comparison of the shooting (s) and the minimization (m)
methods for $\kappa=-1$, $m=c=1$, $V(r)=-\zeta r^{-\beta}$,
$\zeta=0.5$ and $\beta\in (0,1)$. The system (\ref{systemradial}) is
numerically solved with a stepsize adaptative Runge-Kutta method on the
interval $(\epsilon= 10^{-4}, R=15)$. For the shooting method we
minimize the quantity
$\epsilon(|u(\epsilon)|^2+|v(\epsilon)|^2)+\theta(|u(R)|^2+|v(R)|^2) =
\Delta_{{\rm s}}$ for some scale parameter $\theta>0$ (which is chosen to
balance both boundary terms), while for the minimization method, the quantity
$J^+({\scriptstyle
\lambda-1,\phi})$ is directly minimized, the quantity ${|J^+({
\lambda-1,\phi})-(\lambda-1)|^2}$ being computed {\sl a posteriori\/} at
$\,\lambda=\lambda^m$. For
$\beta=1$, the result is known explicitely:
$\lambda_1=[1-\zeta^2]^{1/2}=0.866025...$, $\delta_1=
-[(1-\lambda)/(1+\lambda)]^{1/2}=-0.267949...$}}
\begin{tabular}{cccccccc}
$\beta$ & $\delta^{{\rm s}}$ & $\delta^{{\rm m}}$ &
$\lambda^{{\rm s}}$ & $\lambda^{{\rm m}}$ & $J^+$ &
$|{\scriptstyle(\lambda^m-1)}-J^+({\scriptstyle
\lambda^m-1,\phi})|^2$ & $\Delta^{{\rm s}}$\\
\tableline
1 & -0.267954 & -0.267943 & 0.866034 & 0.866013 &
0.866014 & 1.8$\; $10$^{-12}$ & 0.00029\\
0.9 & -0.235187 & -0.235174 & 0.856725 & 0.856698 &
0.856698 & 2.1$\; $10$^{-14}$ & 0.00053\\
0.8 & -0.207802 & -0.207788 & 0.843181 & 0.843146 &
0.843146 & 5.2$\; $10$^{-14}$ & 0.00063\\
0.7 & -0.183397 & -0.183379 & 0.825877 & 0.825832 &
0.825831 &
4.3$\; $10$^{-13}$ & 0.00076\\ 0.6 & -0.160651 & -0.160627 &
0.804699 & 0.804639 & 0.804639 & 4.1$\; $10$^{-13}$ & 0.00094\\
0.5 & -0.138654 & -0.138619 & 0.779161 & 0.779071 &
0.779070 & 3.4$\; $10$^{-13}$ & 0.0012\\ 0.4 & -0.116645 &
-0.116584 & 0.748381 & 0.748221 & 0.748220 & 3.8$\;
$10$^{-13}$ & 0.0018\\ 0.3 & -0.0938375 & -0.0937016 &
0.710904 & 0.710537 & 0.710536 & 3.5$\; $10$^{-13}$ & 0.0049\\
0.2 & -0.069224 & -0.068798 & 0.664252 & 0.663067 &
0.663067 & 2.4$\; $10$^{-13}$ & 0.0097\\ 0.1 & -0.0412322 &
-0.0392963 & 0.60391 & 0.59833 & 0.59833 & 1.4$\; $10$^{-13}$
& 0.018\\
\end{tabular}
\label{table1}
\end{table}
\narrowtext
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent
Since for a given $\varphi$, $f_\varphi(\lambda)$ is decreasing and
$g_\varphi(\lambda)$ is increasing, if there exits a
$\lambda=\lambda[\varphi]$ such that (\ref{eqnlambda}) is satisfied, then
it is unique (the existence of such a $\lambda$ for all $u$
depends on
the properties of the potential $V$). According to \cite{[Dol-Est-Ser99]},
for those $V$'s,
the ground state is such that
\begin{equation}\label{gsmin}
\lambda_1=\min_{\varphi}\lambda[\varphi]\; .\end{equation}
For a radial potential we may use the radial Dirac equation and
consider (\ref{systemradial}) instead of (\ref{System2spinors}). For
$m=1$ and $c=1$, $\lambda=\lambda_r[u]$ is then the unique solution of
\begin{eqnarray}\label{eqnlambdaradial}
f(\lambda)=&&\int_0^{+\infty}\frac{|(r^\kappa
u)'|^2}{r^{2\kappa}(1+\lambda-V(r))}\; dr \\=&&
(\lambda-1)\int_0^{+\infty}|u(r)|^2\; dr-\int_0^{+\infty} V(r)|u(r)|^2\;
dr\; .\nonumber\end{eqnarray}
To solve it numerically, it is more convenient to rewrite $f(\lambda)$
as
\begin{equation}\label{serielambda}
f(\lambda)=\sum_{k=0}^{+\infty}\biggl[(-1)^k\int_0^{+\infty}
\frac{r^{-2\kappa}|(r^\kappa u)'|^2}{(1-V(r))^{k+1}}\;
dr\biggr]\lambda^k\; .\end{equation}
The solution (with $\kappa=-1$)
is then approximated on a finite basis
$(u_i)_{i=1,2,...n}$:
$u=\sum_{i=1}^n x_iu_i$. If
\begin{equation}f_{ijk}=(-1)^{k-1}\int_0^{+\infty}
\frac{r^2(u_i/r)'(u_j/r)'}{(1-V(r))^k}\; dr\end{equation}
\noindent and
\begin{equation}V_{ij}=\int_0^{+\infty}u_i(r)u_j(r)V(r)\; dr\; ,\end{equation}
\newpage
\noindent
the approximating equation for $\lambda$ corresponding
to (\ref{eqnlambda}) is
\[ %\label{approxlambda}
\sum_{i,j=1}^n \biggl((\sum_{k=1}^mf_{ijk}
\lambda^{k-1})+V_{ij}\biggr)x_ix_j+(1-\lambda)\sum_{i=1}^nx_i^2=0\; ,
\]
where the series in $\lambda$ has been truncated at order $m$. It is actually
more convenient to define
$$ A^{n,m}(\lambda)=\biggl((\sum_{k=1}^mf_{ijk}
\lambda^{k-1})+(1-\lambda)\delta_{ij}+V_{ij}\biggr)_{i,j=1,2,...n}$$
and to approximate $\lambda_1$ by $\lambda_1^{n,m}$ defined as the first
positive root of $\lambda\mapsto\mu(A^{n,m}(\lambda))$ where $\mu(A)$
denotes the first
eigenvalue of the matrix $A$ (see \cite{[Dol-Est-Ser99]} for more details).
Note that $(\lambda_1^{n,m})_{m\geq 1}$ is an alternating sequence (which
essentially converges at a geometric rate): two consecutive eigenvalues
determine an interval containing $\lim_{m\to +\infty}\lambda_1^{n,m}$.
Numerical results for a special basis are given in Table II, with two
types of approximations: a finite basis with $n$ elements is used and the
series in powers of $\lambda$ is truncated at a finite order $m$.
Regarding excited states, they can be obtained by the first and second
methods (shooting and minimization), with $\kappa$ fixed to values
different from $-1$ and for appropriate values for the quantum
numbers corresponding to the decomposition into spherical harmonics.
In the case of the direct minimization method, the
computation of excited states is more straightforward. It is indeed
reduced to the computation of the i$^{th}$ root of $\lambda\mapsto
\mu_i(A^{n,m}(\lambda))$, where $\mu_i(A)$ denotes the
i$^{th}$ eigenvalue of $A$.
In this letter we have proved that computations of one-particle bound
states based on a variational formulation of the Dirac equation are not
subjected to the numerical difficulties due to the negative continuum of
the spectrum, thus showing the possibility of a mathematical foundation
for such numerical methods. These computations are general,
robust and not restricted to central potentials. Optimal numerical
accuracy has not been our
\hbox to\hsize{primary concern, since we were rather interested in}
\vspace{5mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\widetext
\begin{table}
\caption{Direct minimisation method for $\kappa=-1$, $m=c=1$,
$V(r)=-\zeta r^{-\beta}$, $\zeta=0.5$ and $\beta$ close to $1$. The
approximating space is of dimension $n=10$ (we consider the
orthonormal basis generated by the ground
state of the hydrogen atom and $n-1$ Hermite functions, which is probably
not very well adapted) and the series are truncated
at $m=14$ or $m=15$ (the corresponding values $\lambda_1^{10,14}$ and
$\lambda_1^{10,15}$ are respectively a lower and an upper bound of
$\lim_{m\to +\infty}\lambda_1^{10,m}$). This last approximation is
certainly rather crude as shown by the case of the Coulomb potential. As
in Table~I,
$J^+$ is obtained through
a minimization on the set $X$, and ${\,\Delta^{\rm
m}:=({1-\sum_{i=1}^{10}(u^m,
u_i)^2})^{1/2}}\,$ measures the error (in the $L^2$-norm) when
the corresponding solution is
approximated on the finite basis.}
\begin{tabular}{ccccccc}
$\beta$ & 0.90 & 0.93 & 0.95 & 0.97 & 0.99 & 1.00\\
\tableline
$\lambda_1^{10,14}$ & 0.855681 & 0.858516 & 0.860228 &
0.861792
& 0.863200 & 0.863843\\
$\lambda_1^{10,15}$ & 0.858012 & 0.861112 & 0.863004 &
0.864749 & 0.866338 & 0.867071\\
\tableline
$J^+$ & 0.856698 & 0.859984 & 0.861954 &
0.863735 & 0.865310 & 0.866014\\
$\Delta^{{\rm m}}$ & 0.0082 & 0.0058 & 0.0046 & 0.0033 &
0.0020 & 0.0022
\end{tabular}
\label{table2}
\end{table}
\narrowtext
\noindent
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent
showing the feasability of an approach based on the variational structure
of the equation. Hopefully, this new approach will contribute to
the elaboration of more efficient and consistent numerical methods.
\noindent
\begin{references}
\bibitem{[Wal-Ku]} H. Wallmeier, W. Kutzelnigg, Chem. Phys. Lett. {\bf 78}, 341
(1981).
\bibitem{[Baylis-Peel]} W.E. Baylis, S.J. Peel, Phys. Rev. A {\bf 28}, 2552
(1983).
\bibitem{[Ta]} J.D. Talman, Phys. Rev. Lett. {\bf 57}, 1091 (1986).
\bibitem{[Da-De]} S.N. Datta, G. Deviah, Pramana, {\bf 30}, 387
(1988).
\bibitem{[Dr-Go81]} G.W.F. Drake and S.P. Goldman, Phys. Rev. A {\bf 23},
2093 (1981).
\bibitem{[Ku97]} W. Kutzelnigg, Chemical Physics {\bf 225}, 203 (1997).
\bibitem{[Hill-Kra]} R.N. Hill, C. Krauthauser, Phys. Rev. Lett. {\bf
72}, 2151 (1994).
\bibitem{[Dr-Go]} G.W.F. Drake and S.P. Goldman, Adv. Atomic Molecular
Phys. {\bf 25}, 393 (1988).
\bibitem{[Ku84]} W. Kutzelnigg, International J. Quantum Chem. {\bf XXV}, 107
(1984).
\bibitem{[Grant]} I.P. Grant, A.I.P. Conf. Proc. {\bf 189}, 235 (1989).
\bibitem{[Es-Se]} M.J. Esteban, E. S\'er\'e, \it Partial Differential
Equations and Their Applications, \rm CRM Proceedings and Lecture
Notes, {\bf 12}, Eds. P.C. Greiner, V. Ivrii, L.A. Seco and C. Sulem.
AMS (1997).
\bibitem{[Gri-Sie]} M. Griesemer, H. Siedentop, Preprint mp-arc
97-492, to appear in J. London Math. Soc.
\bibitem{[Dol-Est-Ser98]} J. Dolbeault, M.J. Esteban, E.
S{\'e}r{\'e}, Calc. Var. {\bf 10}, 321 (2000).
\bibitem{[Gri-Lew-Sie]} M. Griesemer, R.T. Lewis, H. Siedentop,
Preprint 1999.
\bibitem{[Dol-Est-Ser99]} J. Dolbeault, M.J. Esteban, E.
S{\'e}r{\'e}, J. Funct. Anal. {\bf 174}, 208 (2000).
\bibitem{[LLBS]} R. van Leeuwen, E. van Lenthe, E.J. Baerends, J.G. Snijders,
J. Chem. Phys. {\bf 101}, 1272 (1994).
\bibitem{[Th]} B. Thaller. {\it The Dirac Equation}. Texts
and Monographs in Physics. Springer-Verlag, Berlin, 1st edition (1992).
\bibitem{[Bj-Dr]} J.D. Bjorken, S.D. Drell, Relativistic Quantum
Mechanics, McGraw-Hill (1964).
\end{references}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%