%%%%%%%%% TIENE 13 PAGINAS COMPILADAS --- PAPER 2 --- %%%%%%%%%
%%%%%%%%% Arvesu

\documentstyle[11pt,twoside,epsf]{iwop96}

\textwidth= 16.3cm
\textheight=24cm
\oddsidemargin .5cm
\evensidemargin -.8cm
\topmargin -1cm

\newcommand\cyr{%
\renewcommand\rmdefault{wncyr}%
\renewcommand\sfdefault{wncyss}%
\renewcommand\encodingdefault{OT2}%
\normalfont
\selectfont}
\DeclareTextFontCommand{\textcyr}{\cyr}
\def\cyr{\fontencoding{OT2}\fontfamily{wncyr}\selectfont}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\ry}{\longrightarrow}
\newcommand{\bq}{\begin{equation}}
\newcommand{\eq}{\end{equation}}
\newcommand{\bt}{\begin{tabular}}
\newcommand{\et}{\end{tabular}}
\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}
\newcommand{\br}{\begin{eqnarray}}
\newcommand{\er}{\end{eqnarray}}
\newcommand{\brn}{\begin{eqnarray*}}
\newcommand{\ern}{\end{eqnarray*}}
\newcommand{\benu}{\begin{enumerate}}
\newcommand{\eenu}{\end{enumerate}}
\newcommand{\bite}{\begin{itemize}}
\newcommand{\eite}{\end{itemize}}
\newcommand{\vc}[1]{{\overline{#1}}}
\newcommand{\rc}[1]{{\sqrt{#1}}}
\newcommand{\bpr}{\begin{pr}}
\newcommand{\epr}{\end{pr}}
\newcommand{\dst}{\displaystyle}

\newcommand{\btd}{\bigtriangledown}
\newcommand{\btu}{\bigtriangleup}
\newcommand{\mesp}{\hspace{.3cm}}
\newcommand{\marg}{\vspace*{0.3cm}}
\newcommand{\esp}{\hspace{1cm}}

\newcommand{\fuu}{{\cal U}}
\newcommand{\fuc}{{\cal C}}
\newcommand{\ful}{{\cal L}}
\newcommand{\fuj}{{\cal J}}

\newcommand{\half}{\mbox{\scriptsize $\frac{1}{2}$} }  
\newcommand{\tres}{\mbox{\scriptsize  $\frac{3}{2}$} }  
\newcommand{\funL}{\mbox{ \bf L}}  
\newcommand{\fclas}{{\bf \cal C }} 
 \newcommand{\flineal}{{\bf \cal  L}} 
\newcommand{\degree}{\mbox{degree\,}}
\newcommand{\ener}{ {\varepsilon} }



\newcommand{\RR}{{\sf I}\hspace*{-1.5pt}{\sf R}}
\newcommand{\NN}{{\sf I}\hspace*{-1.5pt}{\sf N}}
\newcommand{\PP}{{\sf I}\hspace*{-1.5pt}{\sf P}}
\newcommand{\QQ}{\hspace*{2pt}{\sf I}\hspace*{-5pt}{\sf Q}}
\newcommand{\ZZ}{{\sf Z}\hspace*{-5pt}{\sf Z}}
\newcommand{\CC}{\hspace*{2pt}{\sf I}\hspace*{-5pt}{\sf C}}
\newcommand{\II}{{\rm I}\hspace*{-1.5pt}{\rm I}}
\newcommand{\KK}{{\sf I}\hspace*{-1.5pt}{\sf K}}
\newcommand{\YY}{\hspace*{-1.5pt}{\sf Y}}
\newcommand{\LL}{{\rm I}\hspace*{-1.5pt}{\rm L}}
\def\bkC{{\rm \kern.24em
             \vrule width.05em height1.4ex depth-.05ex
             \kern-.26em C}}
\def\cyr{\fontencoding{OT2}\fontfamily{wncyr}\selectfont}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% THIS IS LATEX DOCUMENT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\newcommand{\nat}{\NN}
\def\BBox{\vrule height 0.5em width 0.6em depth 0em}
\newcommand{\ga}{\gamma }
\newcommand{\al}{\alpha }
\newcommand{\be}{\beta }
\newcommand{\Ga}{\Gamma }
\newcommand{\e}{\varepsilon}
\newcommand{\z}{\zeta }
\newcommand{\E}{\eta }
\newcommand{\vt}{\vartheta }
\newcommand{\ch}{\chi }
\newcommand{\Oa}{\omega }
\newcommand{\De}{\bigtriangleup }
\newcommand{\vpn}{\mbox{{\bf P}}_{\hspace{-1mm}\mbox{{\scriptsize\bf n}}}}
\newcommand{\gcherta}{(2n+\alpha+\beta+1)\gamma_{n}}
\newcommand{\lrla}{\leftrightarrow }
\newcommand{\vpna}{\mbox{{\bf P}}_{\hspace{-1mm}\mbox{{\scriptsize\bf n}}}^{\hspace{-1mm}\mbox{ {\scriptsize  \bf A}}}  }


\newtheorem{proposition}{Proposition}
\newtheorem{th}{Theorem}
\newtheorem{lema}{Lemma}
\newtheorem{dff}{Definition}

\begin{document}


\pagestyle{myheadings} 
\markboth{\sc \centerline{Proceedings of the IWOP'96} }{\sc \centerline{Proceedings of the IWOP'96} }


\thispagestyle{plain}

\input referencia

\begin{center}

{\large {\bf LAGUERRE POLYNOMIALS IN A RELATIVISTIC QUANTUM-STATISTICAL MODEL.}}\\


{\bf J. Arves\'u }
\end{center}

\noindent % \begin{tabular}{ll}
Key words and phrases:    Dirac equation, orthogonal polynomials, quantum-statistical \\
\phantom{Key words and phrases:}  mechanics. 
%\end{tabular}                      


\begin{abstract}
The main aim of this work is to find analytical expressions for the eigenvalues and
eigenfunctions corresponding to Dirac equation, for hot and dense matter. It is shown that
the Laguerre polynomials depending on the effective charge are solutions for
self-consistent fields. To determine the screening constant and effective charge we introduce
and minimize a functional. These formulas are accurate for machine calculation of bound-bound,
bound-free and free-free transitions, including large values of principal quantum numbers. Hence
these expressions would be in accordance with quantum-statistical results based on more
sophisticated calculations. Comparisons with the solutions of the Schr\"odinger equation
for different substancies are disscused.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%% section 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\setcounter{section}{0}  % comment these two lines out.
\pagenumbering{arabic}   % Tipo de numeros {roman}
\setcounter{page}{23}     %Leave this line commented out.

\section{ Introduction.}

The Knowledge of the wave functions plays a fundamental role in quantum physics, e.g., 
quantum-statistical models, such as opacity and equation of state calculations \cite{af},
\cite{cms}, \cite{ir}, \cite{irw} and \cite{rwi}.\\
When one comes to compute the wave functions and energy levels, 
two methods are usually employed. The first one, which involves sophisticated calculations,
consists in solving the Schr\"odinger or Dirac equation. This is very expensive because the 
model requires very long computational time, mainly when is necessary to take into account the 
multiconfigurability of the heavy atoms, then the use of supercomputers is justified. The 
second one is based on approximate expressions. In light of this development we propose in 
the present work a simple and fast model, which doesn't require to solve the Dirac or 
Schr\"odinger equation, and determines, by means of an analytical expression, the wave functions 
and energy levels for each electron of high-Z elements at high temperatures and densities, 
mainly in the cases, where the computational time can be very long. This model has been 
developed as kernel of future code for opacity calculations and other physical applications.\\ 
The structure of paper is as follows. In Section 2 we provide the required background about the 
solutions of the Dirac equation via Laguerre polynomials, giving analytical expressions for wave 
functions and energy levels, depending on the effective charge and screening constant. 
In Section 3, we define a functional, depending on the effective charge and screening constant. 
Hence we describe the minimization procedure of this functional. The found extremum allows us 
to obtain the wave functions and energy levels, since it depends on the effective charge 
and screening constant. This method is proposed to avoid solving the Dirac equation numerically, 
and so that, we can save a considerable amount of computer memory part. This method is being 
compared with other one \cite{un}. Finally, in Section 4 we provide a table of numerical results 
and some graphics and compare them with results from other methods, in order to determine the 
accuracy of our approximation.
%%%%%%%%%%%%%%%%%%%%%%%%%%% section 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The wave function dependence on effective charge.}
The wave functions corresponding to the numerical solution of the Dirac equation for a
self-consistent field in high temperature plasmas at high-Z elements, have hydrogenic behaviour.
From this behaviour of the wave function suitable to electrons
with any combination of quantum numbers $nlj$, and from their similarity with 
Laguerre functions, the idea of a generalized representation of Hartree, conserving
the spherical symmetry has arised. In this representation the quantum number $j$ has been added.\\
In our approach each atomic electron moves in the individual potential 
$\widetilde{U}(r)=-\frac{Z_{nlj}}r+A_{nlj}$, depending on the nuclear charge and the other
atomic electrons. Here we represent by $Z_{nlj}$ the effective charge and by $A_{nlj}$ the 
screening constant. Although, this assumption is an approximate description for the real hot 
and dense matter. However, it makes possible to understand the behavior of atoms and to find their energy 
levels and wave functions in a very simple way, because it constitutes a very difficult quantum mechanical 
many-body-problem. For the above potential we can solve, analytically, the Dirac equation. 
Therefore the wave functions are expressed in analytical form, depending on $Z_{nlj}$ this  
guarantees to operate with them easily. Thus the efforts shall reduce in order to find the 
appropriate values of the effective charge $Z_{nlj}$ and screening constant $A_{nlj}$, 
respectively.\\
Before proceeding with a description of the method, it is useful to explain which physical 
approximations will be assumed\\
\esp 1. With a good degree of accuracy the atomic nuclei with respect to the atomic electrons
can be consider fixed strenght centers.\\
\esp 2. The electrons are continuosly distributed in the phase space according to Fermi-Dirac statistics.\\
\esp 3. The existence in the atomic nuclei of electric quadrupole, magnetic dipole and all higher 
multipole orders are neglected.\\
The reason for which more work is dedicated to the attainment of the most exactly possible 
wave functions (and not 
justly in the attainment of the energy levels though both are closely linked by means of the 
Dirac equation) 
is conditioned for the calculations of radiative transition probabilities bound-bound,
bound-free and free-free to be incorporated into opacity calculation.
For the purpose of getting a better quantitative understanding of the dependence of the radial 
wave function on the $Z_{nlj}$ and $A_{nlj}$, we summarize here the essential points in the 
procedure to obtain it in a most adequate form. First, we must obtain the exact solution of 
the Dirac differential equation for an electron in a potential 
$\widetilde{U}(r)=-\frac{Z_{nlj}}r+A_{nlj}$, for stationary states with discrete spectra of 
total energy $\widetilde{E}$ it is given by \cite{bs}
$$
\widetilde{\hat H}\widetilde{\Psi}=\left[c(\vec\alpha\vec{\hat p})+mc^2\beta+\widetilde{
U}(r)\right] \widetilde\Psi=\widetilde E\widetilde\Psi,
$$
where $\widetilde\Psi$ is a 1 column-4 row matrix, i.e., the spinor of four components
$$
\widetilde\Psi=\left(\matrix{\widetilde\psi_1\cr\widetilde
	\psi_2\cr\widetilde\psi_3\cr\widetilde\psi_4\cr}\right),\quad(\vec\alpha
	\vec{\hat p})=\sum^3_{k=1}\alpha_k\hat p_k,\quad
\left(\hat p_k=\frac\hbar i\frac\partial{\partial x_k}\right),\quad
\alpha_k=\left(\matrix{0&\sigma_k\cr\sigma_k&0\cr}\right),\quad
\beta=\left(\matrix{I&0\cr0&I\cr}\right),
$$
being the components of $\sigma_k$ and $\beta$ the following two-by-two matrices

$$
\sigma_1=\left(\matrix{0&1\cr1&0\cr}\right),\quad
\sigma_2=\left(\matrix{0&-i\cr i&0\cr}\right),\quad
\sigma_3=\left(\matrix{1&0\cr0&-1\cr}\right),\quad \mbox{\,and\,}\quad
I=\left(\matrix{1&0\cr0&1\cr}\right),\quad
0=\left(\matrix{0&0\cr0&0\cr}\right).
$$
The symbols with tilde are used in this work to distinguish both the hydrogenic-like
magnitudes and wave functions from those wave functions corresponding to the 
exact solution of the Dirac equation for the hydrogen.\\
Writing explicitly the Dirac equation for the four components of the wave 
function, taking into account relativistic units, i.e., the mass of the electron, Planck's
constant $\hbar$, and the speed of light are all equal to 1, it holds 
$$
i\left(\widetilde \varepsilon+\frac{\widetilde\mu}r+1\right)\widetilde\psi_1+
\frac{\partial\widetilde\psi_3}{\partial z}+
\left(\frac\partial{\partial x}-i\frac\partial{\partial y}\right)
\widetilde\psi_4=0,\quad
$$

$$
i\left(\widetilde \varepsilon+\frac{\widetilde\mu}r+1\right)\widetilde\psi_2-
\frac{\partial\widetilde\psi_4}{\partial z}+
\left(\frac\partial{\partial x}+i\frac\partial{\partial y}\right)
\widetilde\psi_3=0,\quad
$$
\begin{equation}
i\left(\widetilde \varepsilon+\frac{\widetilde\mu}r-1\right)\widetilde\psi_3+
\frac{\partial\widetilde\psi_1}{\partial z}+
\left(\frac\partial{\partial x}-i\frac\partial{\partial y}\right)
\widetilde\psi_2=0,\quad
\label{1}
\end{equation}
$$
i\left(\widetilde \varepsilon+\frac{\widetilde\mu}r-1\right)\widetilde\psi_4-
\frac{\partial\widetilde\psi_2}{\partial z}+
\left(\frac\partial{\partial x}+i\frac\partial{\partial y}\right)
\widetilde\psi_1=0,\quad
$$
where
$$
\ba{c}
\dst\widetilde\mu=\widetilde\mu_{nlj}=\frac{Z_{nlj}e^2}{\hbar c}
\approx\frac{Z_{nlj}}{137},\quad
\widetilde E_{nlj}=\widetilde\varepsilon+A_{nlj},
\quad (0<\widetilde E_{nlj} <1)\quad.
\ea
$$
In spherical coordinates $(r,\ \theta,\ \varphi)$ the variables in these equations can be
separated --if one makes use of the fact that the components and the absolute value of 
the total angular momentum of an electron in any central force field are constants of the motion, 
then we can construct the eigenfunction of the eigenvalues which correspond to operators
$\ \hat H,\ \hat J^2$ and $\hat J_{Z}$, i.e, we base on the spherical symmetry the field 
of the nucleus-- by looking for a solution that has the form \cite{nsu}, \cite{nu}
$$
\ba{c}
\dst\left(\matrix{\widetilde\psi_1\cr\widetilde
	\psi_2\cr}\right)=\widetilde f_{nlj}(r)\Theta_{nlm}(\theta,\varphi),\quad
\left(\matrix{\widetilde\psi_3\cr\widetilde
	\psi_4\cr}\right)=(-1)^{\left(l-l'+1\right)/2}\widetilde g_{nlj}(r)
	\Theta_{nl'm}(\theta,\varphi),\quad
\ea
$$
whence
\begin{equation}
\ba{l}
%\displaystyle
\widetilde\Psi=\left(\matrix{\widetilde\psi_1\cr\widetilde
	\psi_2\cr\widetilde\psi_3\cr\widetilde\psi_4\cr}\right)=
	(-1)^{\left(m-\frac12\right)}i^l
\left(\matrix{\sqrt{\frac{j+m}{2l+1}}&\widetilde f_{nlj}(r)&
\Theta_{l,m-\frac12}(\theta,\varphi)\cr
-\sqrt{\frac{j-m}{2l+1}}&\widetilde f_{nlj}(r)&
\Theta_{l,m+\frac12}(\theta,\varphi)\cr
-\sqrt{\frac{j-m+1}{2l+3}}&i\widetilde g_{nlj}(r)&
\Theta_{l+1,m-\frac12}(\theta,\varphi)\cr
-\sqrt{\frac{j+m+1}{2l+3}}&i\widetilde g_{nlj}(r)&
\Theta_{l+1,m+\frac12}(\theta,\varphi)\cr}\right)\quad
\mbox{\,for\,}\quad l=j-\frac12,\quad\\
\\
\widetilde\Psi=\left(\matrix{\widetilde\psi_1\cr\widetilde
	\psi_2\cr\widetilde\psi_3\cr\widetilde\psi_4\cr}\right)=
	(-1)^{\left(m+\frac12\right)}i^l
\left(\matrix{\sqrt{\frac{j-m+1}{2l+1}}&\widetilde f_{nlj}(r)&
\Theta_{l,m-\frac12}(\theta,\varphi)\cr
\sqrt{\frac{j+m+1}{2l+1}}&\widetilde f_{nlj}(r)&
\Theta_{l,m+\frac12}(\theta,\varphi)\cr
\sqrt{\frac{j+m}{2l-1}}&i\widetilde g_{nlj}(r)&
\Theta_{l-1,m-\frac12}(\theta,\varphi)\cr
-\sqrt{\frac{j-m}{2l-1}}&i\widetilde g_{nlj}(r)&
\Theta_{l-1,m+\frac12}(\theta,\varphi)\cr}\right)\quad
\mbox{\,for\,}\quad l=j+\frac12.\quad
\ea
\label{2}
\end{equation}
	By substituting the explicit formulas (\ref{2}) in (\ref{1}) one then obtains \cite{nu} 
\begin{equation}
{\cases{\widetilde f^\prime_{nlj}+\frac{1+\kappa}r\widetilde f_{nlj}-
\left(\widetilde \varepsilon_{nlj}+1+
\frac{\widetilde\mu_{nlj}}r\right)\widetilde g_{nlj}=0,\cr
{\widetilde g^\prime_{nlj}+\frac{1-\kappa}r\widetilde g_{nlj}-
\left(\widetilde \varepsilon_{nlj}-1+
\frac{\widetilde\mu_{nlj}}r\right)\widetilde f_{nlj}=0,\quad}\cr}}
\label{6}
\end{equation}
where
$$
\kappa={\cases{-(l+1)&$\mbox{\,if\,}\quad l=j-\frac12$,\cr
\quad \hspace{0.3cm} l&$\mbox{\,if\,}\quad l=j+\frac12$.\cr}},\quad l=0,1,...,n-1\quad.
$$
For bound states the wave functions $\widetilde f_{nlj}(r)$ and $\widetilde g_{nlj}(r)$ should 
satisfy the normalization condition
\bq
\ba{c}
\dst\int^{\infty}_0\left[\widetilde f^2_{nlj}(r)+\widetilde g^2_{nlj}(r)\right]r^2\,dr=1,\quad
\mbox{\,since\,}\quad
\int^{\pi}_0\int^{2\pi}_0\left|\Theta_{l,m}(\theta,\varphi)\right|^2\,d\theta d\varphi=1.\quad
\ea
\label{normalization}
\eq
The system (\ref{6}) can be written in matrix form
\begin{equation}
\dst\left(\matrix{\widetilde f'_{nlj}(r)\cr \widetilde g'_{nlj}(r)\cr}\right)=\KK \cdot 
\left(\matrix{\widetilde f_{nlj}(r)\cr\widetilde g_{nlj}(r)\cr}\right),\quad
\mbox{\,where\,}\quad
\KK=\left(\matrix{\dst -\frac{1+k}{r}& \dst 1+\widetilde
\varepsilon +\dst \frac{\widetilde\mu}{r}\cr\dst  1-\widetilde\varepsilon - 
\frac{\widetilde\mu}{r}& -\dst \frac{1-k}{r}\cr}\right).\quad
\label{u'through u}
\eq
This leads to the following system of second order differential equations
\begin{equation}
\ba{l}
\dst \widetilde f_{nlj}''(r)-\tau_{1}(r)\widetilde f'_{nlj}(r)+\lambda_{1}(r)\widetilde f_{nlj}(r)=0,\\
\\
%Similary, eliminating $u_1(r)$, we obtain an equation for $u_2(r)$:
\dst \widetilde g''_{nlj}(r)-\tau_{2}(r)\widetilde g'_{nlj}(r)+\lambda_{2}(r)\widetilde g_{nlj}(r)=0,\quad
\ea
\label{eq. of second order for u2}
\end{equation}
where
\begin{equation}
\ba{l}
\tau_{1}(r)=-\left[\frac{2}{r}+
\frac{\widetilde\mu}{r^2(1+\widetilde\varepsilon)+r\widetilde\mu}\right],\,
\lambda_{1}(r)=\frac{1-k^2}{r}-\left[1-(\widetilde\varepsilon+
\frac{\widetilde\mu}{r})^2\right]-\frac{1+k}{r^2}+\frac{\widetilde\mu(1+k)}{r^3(1+
\widetilde\varepsilon)+r^2\widetilde\mu},\\
\\
\tau_{2}(r)=-\left[\frac{2}{r}-
\frac{\widetilde\mu}{r^2(1-\widetilde\varepsilon)-r\widetilde\mu}\right],\,
\lambda_{2}(r)=\frac{1-k^2}{r}-\left[1-(\widetilde\varepsilon+
\frac{\widetilde\mu}{r})^2\right]-\frac{1-k}{r^2}-\frac{\widetilde\mu(1-k)}{r^3(1+
\widetilde\varepsilon)+r^2\widetilde\mu}.
\ea
\label{tau-lambda}
\end{equation}
The equations (\ref{eq. of second order for u2}) are not
generalized equations of hypergeometric type, since $1+\widetilde\varepsilon\not=0$ and 
$\widetilde\mu\not=0$.\\
By means of a linear transformation the system (\ref{eq. of second order for u2}) becomes 
in a similar one, which under certain conditions can be transformed in a generalized system 
of equations of hypergeometric type. Let be the system
$$
\LL\cdot\left(\matrix{\widetilde f_{nlj}(r)\cr\widetilde g_{nlj}(r)\cr}\right)=
\widetilde{\YY}=\left(\matrix{\widetilde y_1\cr \widetilde y_2\cr}\right),\quad
$$
where $\LL$ is a nonsingular matrix independent of $r$. Therefore
\begin{equation}
\ba{c}
\dst\widetilde{\YY}'=\left(\matrix{\widetilde y'_1\cr \widetilde y'_2\cr}\right)=
\widetilde{\KK}\cdot\widetilde{\YY},\mbox{\,where\,\,}
\widetilde{\KK}=\LL\cdot\KK\cdot\LL^{-1}=\left(\matrix{\widetilde k_{11}& \widetilde k_{12}\cr 
\widetilde k_{21}& \widetilde k_{22}\cr}\right),\mbox{\,and\,\,}
\LL=\left(\matrix{\widetilde l_{11}& \widetilde l_{12}\cr 
\widetilde l_{21}& \widetilde l_{22}\cr}\right).\quad
\ea
\label{v'through v}
\eq
The elements of $\widetilde{\KK}$ are
\begin{equation}
\ba{c}
\dst \widetilde{k}_{11}=\dst (\det\LL)^{-1}\left[-\frac{1+k}{r}l_{11}l_{22}-(1+\widetilde\varepsilon+
\frac{\widetilde\mu}{r})l_{11}l_{21}+(1-\widetilde\varepsilon-\frac{\widetilde\mu}{r})l_{12}l_{22}+
\frac{1-k}{r}l_{12}l_{21}\right],\\
\\
\dst \widetilde{k}_{12}=\dst (\det\LL)^{-1}\left[(1+\widetilde\varepsilon+\frac{\widetilde\mu}{r})l^{2}_{11}-
(1-\widetilde\varepsilon-\frac{\widetilde\mu}{r})l^{2}_{12}+\frac{2k}{r}l_{11}l_{12}\right],\\
\\
\dst \widetilde{k}_{21}=\dst (\det\LL)^{-1}\left[(1-\widetilde\varepsilon-\frac{\widetilde\mu}{r})l^{2}_{22}
-(1+\widetilde\varepsilon+\frac{\widetilde\mu}{r})l^{2}_{21}-\frac{2k}{r}l_{21}l_{22}\right],\\
\\
\dst \widetilde{k}_{22}=\dst (\det\LL)^{-1}\left[\frac{1+k}{r}l_{12}l_{21}+(1+\widetilde\varepsilon+\frac{\widetilde\mu}{r})l_{11}l_{21}-
(1-\widetilde\varepsilon-\frac{\widetilde\mu}{r})l_{12}l_{22}-\frac{1-k}{r}l_{11}l_{22}\right].\quad
\ea
\label{matriz-sistema}
\eq
From the equations for $\widetilde{y}_1(r)$ and $\widetilde{y}_2(r)$ we get
\begin{equation}
\dst\widetilde{y}''_1(r)-\widetilde{\tau}_{1}(r)\widetilde{y}'_1(r)+\widetilde{\lambda}_{1}(r)
\widetilde{y}_{1}(r)=0,\quad
\label{eq. of second order for v1}
\end{equation}
\begin{equation}
\dst\widetilde{y}''_2(r)-\widetilde{\tau}_{2}(r)\widetilde{y}'_2(r)+\widetilde{\lambda}_{2}(r)
\widetilde{y}_{2}(r)=0,\quad
\label{eq. of second order for v2}
\end{equation}
where
\bq
\ba{c}
\dst\widetilde{\tau}_{1}(r)=\dst\left(\widetilde k_{11}+\widetilde k_{22}+
\frac{\widetilde k'_{12}}{\widetilde k_{12}}\right),\quad
\widetilde{\lambda}_{1}(r)=\dst\left(\widetilde{k}_{11}\widetilde{k}_{22}-\widetilde{k}_{12}\widetilde{k}_{21}-\widetilde{k}'_{11}+
\frac{\widetilde{k}'_{12}}{\widetilde{k}_{12}}\widetilde{k}_{11}\right),\\
\\
\dst\widetilde{\tau}_{2}(r)=\dst\left(\widetilde k_{11}+\widetilde k_{22}+
\frac{\widetilde k'_{21}}{\widetilde k_{21}}\right),\quad
\widetilde{\lambda}_{1}(r)=\dst\left(\widetilde k_{11}\widetilde k_{22}-\widetilde k_{12}\widetilde k_{21}-\widetilde k'_{22}+
\frac{\widetilde k'_{21}}{\widetilde k_{21}}\widetilde k_{22}\right).
\ea
\eq
The equations (\ref{eq. of second order for v1}) and (\ref{eq. of second order for v2}) 
become into equations of generalized hypergeometric type if we impose the conditions
$$
2kl_{11}l_{12}+\widetilde\mu(l^2_{12}+l^2_{12})=0\quad\mbox{\,and\,}\quad
2kl_{21}l_{22}+\widetilde\mu(l^2_{21}+l^2_{22})=0.\quad
$$
Notice that we select two among four possibilities to transform (\ref{eq. of second order for v1})
and (\ref{eq. of second order for v2}) into equations of generalized hypergeometric type.
Then, $\LL$ takes the form
\bq
\LL=\left(\matrix{\widetilde\mu& \widetilde\nu-k\cr \widetilde\nu-k& \widetilde\mu\cr}\right),
\mbox{\,where\,\,} \widetilde\nu=\sqrt{\kappa^2-\widetilde\mu^2}.
\label{siete}
\eq
In the case when $1+\frac{\varepsilon\kappa}{\widetilde\nu}\not=0$, we obtain the following 
differential equation for $\widetilde y_1(r)$
\begin{equation}
\widetilde y''_1(r)+\frac{2}{r}\widetilde y'_1(r)+\frac{(\widetilde\varepsilon^2-1)r^2+
2\widetilde\varepsilon\mu r-\widetilde\nu(\widetilde\nu+1)}{r^2}\widetilde y_1(r)=0.\quad
\label{eq. diff. of second order for v1}
\end{equation}
In the opposite cases $1+\frac{\varepsilon\kappa}{\widetilde\nu}=0$, then we deduce that
that $\widetilde y_1(r)=0$ and $\widetilde y_2(r)=0$. This kind of functions don't satisfy 
(\ref{normalization}). In order to determine the behavior of $\widetilde y_1(r)$ in the neighborhood 
of $r=0$ we require
$$
|(\widetilde\varepsilon^2-1)r^2+2\widetilde\varepsilon\widetilde\mu r|\ll\widetilde\nu(\widetilde\nu+1),\quad
$$
and then $\widetilde y_1(r)\approx C_1r^{\widetilde\nu}$ for small enough $r$.\\
Equation (\ref{eq. diff. of second order for v1}) is a generalized equation of hypergeometric type. 
In order to (\ref{eq. diff. of second order for v1}) will be an equation of hypergeometric type
$$
\sigma (r)y''+\tau (r)y'+\lambda y=0,\quad
$$
where, $\sigma (r)$ and $\tau (r)$ are polynomials of degree at most 2 and 1, respectively; we substitute 
$\widetilde y_1$ by $\phi (r)y$, where $\phi (r)$ satisfies
$$
\phi^{\prime}/\phi=\pi (r)/\sigma (r).\quad
$$
We choose $\tau (r)=\widetilde\tau (r)+2\pi (r)$ with negative derivative and zero on ($0$,$+\infty$). 
These conditions impose restrictions on the choice of $\pi (r)$, hence we get 
$\pi (r)=-\widetilde ar+\widetilde\nu$, with $\widetilde a=\sqrt{1-\widetilde\varepsilon^{2}}$.\\
Whence
$$
\begin{array}{cccc}
\sigma (r)=r,\quad\tau (r)=2(\widetilde\nu +1-\widetilde ar),\quad\phi (r)=r^{\widetilde\nu}e^{-\widetilde ar},\quad
\lambda=2[\widetilde\varepsilon\widetilde\mu-(\widetilde\nu+1)\widetilde a],\quad\rho (r)=r^{2\widetilde\nu+1}e^{-2\widetilde ar}.
\end{array}
$$
The eigenvalues are determined by \cite{nu}
$$
\lambda+n\tau '+\frac12n(n-1)\sigma ''=0,\quad n=0,1,2,...
$$
and the eigenfunctions follow from the Rodrigues formula \cite{nu}
\begin{equation}
y_n(r)=\frac{C_n}{\rho (r)}\frac{d^n}{dr^n}[\sigma^{n}(r)\rho (r)]=C_nr^{-2\widetilde\nu-1}e^{2\widetilde ar}
\frac{d^n}{dr^n}(r^{n+2\widetilde\nu+1}e^{-2\widetilde ar}).\quad
\label{Rodrigues formula}
\end{equation}
The functions $y_n(r)$ are, up to a constant factor, the Laguerre polynomials $L_{n}^{2\widetilde\nu+1}(x)$,
with $x=2\widetilde ar$.\\
In order to obtain the equivalence with quantum mechanic interpretation we replace $n$ by $n-1$, then
\begin{equation}
\widetilde\mu\widetilde\varepsilon-(n+\widetilde\nu+1)\widetilde a=0.\quad
\label{ eq. for tree values}
\end{equation}
As a consequence the eigenfuctions $\widetilde y_1(r)$ take the form
\begin{equation}
\widetilde y_1(r)={\cases{A_nx^{\widetilde\nu}e^{-x/2}L_{n-1}^{2\widetilde\nu+1}(x) & {\mbox {\em n}=1,2,...},\cr
\hspace{1cm} 0 & \quad {\mbox {\em n}=0}.\quad\cr}}
\label{ form of v1}
\end{equation}
Using (\ref{eq. of second order for v2}), (\ref{ eq. for tree values}) and (\ref{ form of v1}) we find a similar 
differential equation for $\widetilde y_2$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{equation}
\widetilde y''_{2}(r)+\frac{2}{r}\widetilde y'_2+\frac{(\widetilde\varepsilon^2-1)r^2+2\widetilde\varepsilon\mu r+
\widetilde\nu(\widetilde\nu-1)}{r^2}\widetilde y_2=0.\quad
\label{eq. diff. of second order for v2}
\end{equation}
From (\ref{eq. diff. of second order for v2}) it follows that the eigenfunctions are
\begin{equation}
\widetilde y_2(r)=B_nx^{\widetilde\nu-1}e^{-x/2}L_{n}^{2\widetilde\nu-1}(x).\quad
\label{expression for v2}
\end{equation}
Using (\ref{v'through v}), (\ref{matriz-sistema}) and (\ref{siete}) we find 
\bq
\widetilde y'_1=\left(-\frac{\widetilde\nu +1}{r}+\frac{\widetilde\varepsilon\kappa}{\widetilde\nu}\right)
\widetilde y_1+\left(1+\frac{\widetilde\varepsilon\kappa}{\widetilde\nu}\right)\widetilde y_2,\quad
\label{relation-y1-y2}
\eq
putting the explicit formula for $\widetilde y_1$ as well as $r=0$ in (\ref{relation-y1-y2}) and taking into 
account that \cite{nu}
$$
L_n^{\alpha}(0)=\frac{\Gamma (n+\alpha+1)}{n!\Gamma (\alpha+1)}.\quad
$$
We have
$$
A_n=\frac{\widetilde\nu+\widetilde\varepsilon\kappa}{\widetilde an(n+
2\widetilde\nu)}B_n \ \ \ \ (n=1,2,...).\quad
$$
Since we know $\widetilde y_1(r)$ and $\widetilde y_2(r)$, we can give the analytic expression for the radial
components of the wave functions $\widetilde f_{nlj}(r)$ and $\widetilde g_{nlj}(r)$
\begin{equation}
\left(\matrix{\widetilde f_{nlj}\cr \widetilde g_{nlj}\cr}\right)=\frac{1}{2\widetilde\nu(\kappa-\widetilde\nu)}
\left(\matrix{\widetilde\mu & \widetilde\nu-\kappa\cr \widetilde\nu-\kappa& 
\widetilde\mu\cr}\right)\cdot
\left(\matrix{\widetilde y_1\cr \widetilde y_2\cr}\right),\quad
\end{equation}
from which
\bq
\ba{c}
\dst\widetilde f_{nlj}(r)\biggm|_{r=\frac x{2 \widetilde a}}=
\frac{ B_n}{2 \widetilde\nu(\kappa- \widetilde\nu)}
x^{ \widetilde\nu-1}e^{-\frac x2}
\left( f_1xL^{2 \widetilde\nu+1}_{n-1}(x)+
f_2L^{2 \widetilde\nu-1}_n(x)\right),\\
\\
\dst\widetilde g_{nlj}(r)\biggm|_{r=\frac x{2 \widetilde a}}=
\frac{ B_n}{2 \widetilde\nu(\kappa- \widetilde\nu)}
x^{ \widetilde\nu-1}e^{-\frac x2}
\left( g_1xL^{2 \widetilde\nu+1}_{n-1}(x)+
g_2L^{2 \widetilde\nu-1}_n(x)\right),
\ea
\label{f-g}
\eq
where
$$
\ba{l}
\dst B_n=2 \widetilde a\sqrt{\frac{(\kappa-
\widetilde\nu)( \widetilde\varepsilon\kappa-
\widetilde\nu)n!}{ \widetilde\mu\Gamma(n+2 \widetilde\nu)}},\quad
\dst \widetilde\varepsilon=
\frac1{\sqrt{1+\left(\frac{ \widetilde\mu}{n+ \widetilde\nu}\right)^2}},\quad
\dst \widetilde\nu=\sqrt{(j+\frac12)^2- \widetilde\mu^2},\\
\\
\widetilde a=\frac{\frac{ \widetilde\mu}{n+ \widetilde\nu}}
{\sqrt{1+\left(\frac{ \widetilde\mu}{n+ \widetilde\nu}\right)^2}},\quad
\dst f_1=\frac{ \widetilde a \widetilde\mu}{ \widetilde\varepsilon\kappa- \widetilde\nu},\quad
\dst f_2=\kappa- \widetilde\nu,\quad
\dst g_1=
\frac{ \widetilde a(\kappa- \widetilde\nu)}{ \widetilde\varepsilon\kappa- \widetilde\nu},\quad
\dst g_2= \widetilde\mu.\quad
\ea
$$
To calculate $B_n$ we have used the orthonormality condition, the three-term recurrence relation and 
the following expansion formula for the Laguerre polynomials \cite{nu}
$$
L^{\al-2}_{n}(x)=L^{\al}_{n}(x)-2L^{\al}_{n-1}(x)+...
$$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The solutions (\ref{f-g}) in atomic units take the form
$$
\dst \widetilde F_{nlj}(r)\biggm|_{r=\frac x{137\cdot2\widetilde a}}=r\widetilde f_{nlj}(r)=
(137)^{\frac12}\frac{\widetilde B_{nlj}}{4\widetilde a\widetilde\nu(\kappa-\widetilde\nu)}
x^{\widetilde\nu}e^{-\frac x2}
\left(\widetilde f_1xL^{2\widetilde\nu+1}_{n-j-\frac32}(x)+
\widetilde f_2L^{2\widetilde\nu-1}_{n-j-\frac12}(x)\right),
$$
$$
\dst \widetilde G_{nlj}(r)\biggm|_{r=\frac x{137\cdot2\widetilde a}}=r\widetilde g_{nlj}(r)=
(137)^{\frac12}\frac{\widetilde B_{nlj}}{4\widetilde a\widetilde\nu(\kappa-\widetilde\nu)}
x^{\widetilde\nu}e^{-\frac x2}
\left(\widetilde g_1xL^{2\widetilde\nu+1}_{n-j-\frac32}(x)+
\widetilde g_2L^{2\widetilde\nu-1}_{n-j-\frac12}(x)\right),
$$
$$
\ba{c}
\dst \widetilde B_{nlj}=\frac{2\widetilde a^2}{\widetilde\mu}
\sqrt{\frac{(\kappa-\widetilde\nu)(n-j-\frac12)!
[\widetilde a\kappa(n-j-\frac12+\widetilde\nu)-\widetilde\mu\widetilde\nu]}
{\Gamma(n-j-\frac12+2\widetilde\nu)}},\quad
\widetilde a=\frac{\frac{\widetilde \mu}{n-j-\frac12+\widetilde\nu}}{\sqrt{1+
\left(\frac{\widetilde \mu}{n-j-\frac12+\widetilde\nu}\right)^2}},\\
\\
\dst \widetilde f_1=
\frac{\widetilde a\widetilde\mu^2}{\widetilde a\kappa(n-j-\frac12+\widetilde\nu)-
\widetilde\mu\widetilde\nu},\quad
\widetilde f_2=\kappa-\widetilde\nu,\quad
\widetilde g_1=\frac{\kappa-\widetilde\nu}{\widetilde\mu}\widetilde f_1,\quad
\widetilde g_2=\widetilde\mu,\\
\\
\dst \widetilde n=
n-\frac{\widetilde\mu^2}{(j+\frac12)+\sqrt{(j+\frac12)^2-\widetilde\mu^2}},\quad
\widetilde\varepsilon_{nlj}=
-\frac{\frac{Z^2}{\widetilde n^2}}{1+\frac{\widetilde\mu^2}{\widetilde n^2}+
\sqrt{1+\frac{\widetilde\mu^2}{\widetilde n^2}}},
\widetilde\nu=\sqrt{(j+\frac12)^2-\widetilde\mu^2},\quad\\
\\
\dst n=1,2,...\ ,\quad j=\frac12,\frac32,...,n-\frac12\quad.
\ea
$$
Notice that we have replaced $\widetilde f_{nlj}(r)$ and $\widetilde g_{nlj}(r)$ 
by $\widetilde F_{nlj}(r)$ and $\widetilde G_{nlj}(r)$, respectively to remark the 
difference between the wave function in relativistic units and the wave function 
in atomic units.\\
From the above expressions it is clear that whether we are able to find
the effective charge we can construct the atomic wave functions. In the next section we
propose a simple method to determine it.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Calculation of the effective charge.}
At this point it is important  to construct and minimize a 
functional, whose dependence is given in terms of the variables $Z_{nlj}$ and $A_{nlj}$,
so the wave functions and energy levels can be expressed analytically for a
self-consistent field \cite{au}. This procedure constitutes an alternative method to avoid solving 
numerically the Dirac equation for each state at given potential.\\
When one solves a quantum-statistical problem we calculate the Thomas-Fermi 
potential $U_{T-F}(r)$ or the Hartree-Fock-Slater $U_{H-F-S}(r)$ one. Both potentials depend 
on the temperature ($\theta$) and density ($\rho$). The product $rU(r)$
is approximated by a linear function on a restricted interval of $r$ in which the absolute 
maximum or minimum of the wave function is located, i.e., in the region of $r$ where the 
wave function certainly is different from zero (the most probable position for the 
particle). Let be take the following linear function,\\
$$
F(r)\equiv Z_{nlj}-rA_{nlj}=-r\widetilde U(r),\quad r\in [0,+\infty[.\quad
$$
All this can be expressed mathematically by means of a functional $\Phi(Z_{nlj},A_{nlj})$,
\begin{equation}
\Phi(Z_{nlj},A_{nlj})=\int^{\infty}_0\left[rU(r)-
r\widetilde U(r,Z_{nlj},A_{nlj})\right]^2
\left[\widetilde F^2_{nlj}(r)+\widetilde G^2_{nlj}(r)\right]\,dr.\quad
\label{0}
\end{equation}
Whereby the radial components of the wave functions corresponding to the self-consistent 
field $U(r)$ will be almost equal to the radial components of the analytic wave functions 
$\widetilde F_{nlj}(r)$ and $\widetilde G_{nlj}(r)$ for the parametric potential $\widetilde U(r)$.\\
Notice that the behavior of
the curve $rU(r)$ depends on the choice 
range of temperature and density. Although this model is applicable on the restricted range of 
temperature and density our numerical simulations shown that it cover a wide region of densities 
and temperatures: $10^{-3} g/cm^{3}\leq\rho\leq10^{2} g/cm^{3}$, $5 ev\leq\theta\leq15 Kev$, 
depending on the substance and quantum number combinations.\\
Without a preliminary analysis the functional $\Phi(Z_{nlj},A_{nlj})$, can be minimized
under conditions $\frac{\partial\Phi(Z_{nlj},A_{nlj})}{\partial Z_{nlj}}=0$ and 
$\frac{\partial\Phi(Z_{nlj},A_{nlj})}{\partial A_{nlj}}=0$. Once values 
$Z_{nlj}$ and $A_{nlj}$ are found to construct the wave function, but we have prefered to use 
other way, consisting in the addition of a boundary condition.
This procedure guarantees a good agreement between all energy levels corresponding to $U(r)$ and 
$\widetilde U(r)$, respectively.\\
The boundary condition is
\begin{equation}
\int^{\infty}_0\left[U(r)-\widetilde U(r,Z_{nlj},A_{nlj})\right]
\left[\widetilde F^2_{nlj}(r)+\widetilde G^2_{nlj}(r)\right]\,dr=0.
%\footnote{This boundary condition 
%is equivalente to consider in terms of perturbation theory that: $\varepsilon=
%\widetilde\varepsilon+\triangle\varepsilon+..., \triangle\varepsilon \rightarrow 0$.}
\label{3}
\end{equation}
If this boundary condition is assumed, then $A_{nlj}$ may be expressed
through $Z_{nlj}$ as
$$
A_{nlj}=2\widetilde\varepsilon_{nlj}+\int_0^\infty U(r)\left[\widetilde F_{nlj}(r)+
\widetilde G_{nlj}(r)\right]\,dr.
$$
Here the fact that the potential energy $-\frac{Z_{nlj}}{r}$ is two times the 
total energy $\widetilde\varepsilon_{nlj}$ is used.\\
Substituting $A_{nlj}$ into the expression (\ref{0}), it is observed that $\Phi(Z_{nlj},A_{nlj})$
will depend only on $Z_{nlj}$
\begin{equation}
\Phi(Z_{nlj})=\int^\infty_0\left[rU(r)+Z_{nlj}-2r\widetilde\varepsilon_{nlj}-
rI\right]^2\left[\widetilde F^2_{nlj}(r)+\widetilde G^2_{nlj}(r)\right]\,dr,\quad
\label{17}
\end{equation} 
where
$$
I\equiv\ I(Z_{nlj})=
\int_0^\infty U(r)\left[\widetilde F_{nlj}(r)+\widetilde G_{nlj}(r)\right]\,dr.
$$
To determine the value $Z_{nlj}$ for which $\Phi(Z_{nlj})$, given by (\ref{17}), has a minimum 
we use a iterative minimization method \cite{dy}. This method quickly converges, with the
prescribed accuracy (in few iterative steps), if we start the iteration near an unknown $Z_{nlj}$.
In such a way we can choose as initial value of iteration the effective charge for 
nonrelativistic case $Z_{nl}$, \cite{un}, as well as the Thomas-Fermi potential \cite{lat}, 
for $U(r)$, i.e., the functional (\ref{17}) has a minimum in the neighborhood of $Z_{nl}$. The
reason for such a choice is that the properties of matter for the considered range of 
temperatures and densities is approximately described by nonrelativistic quantum physics. We want
to remark here that, in order to do this procedure self-consistent, it is necessary to use, in 
the next iterations, the hydrogenic-like potential (instead of the Thomas-Fermi one).\\
For the minimization procedure it is assumed that we are able to find the initial value 
$Z_{nl}$. Let consider that the fine structure constant $\alpha=\frac{e^2}{\hbar c}$ is approximately equal 
to zero (extreme non relativistic limit). Then we have the nonrelativistic case \cite{au},  
\cite{un}
\begin{equation}
\dst-\frac12\widetilde R''(r)+\left[\widetilde U(r)+\frac{l(l+1)}{2r^2}\right]\widetilde R=
\widetilde\varepsilon\tilde R,
\mbox{\,where\,\,}\widetilde U(r)=-\frac{Z_{nl}}r+A_{nl}.\quad
\label{schrodinger}
\end{equation}
The functional analogue to (\ref{17}) is, in this situation
\bq
\dst\widetilde\Phi(Z_{nl},A_{nl})=\int^{\infty}_0\left[rU(r)-r\tilde
U(r,Z_{nl},A_{nl})\right]^2\widetilde R_{nl}^2(r,Z_{nl})\,dr.
\label{func-nl}
\eq
Moreover, the boundary condition have the form
\begin{equation}
\dst\int^{\infty}_0\left[U(r)-\widetilde U(r)\right]\widetilde R^2_{nl}(r)\,dr=0.
\label{frontera-nl}
\end{equation}
The energy levels and radial components for Schr\"odinger equation are well known \cite{bs}, \cite{un}. 
In fact we get
\begin{equation}
\dst\widetilde\varepsilon_{nl}=-\frac{Z^2_{nl}}{2n^2}+A_{nl},
\quad\mbox{\,and\,}\quad\widetilde
R_{nl}(r)=\frac{\sqrt{Z_{nl}}}n\sqrt{\frac{(n-l-1)!}{(n+l)!}}x^{l+1}
e^{-\frac x2}L^{2l+1}_{n-l-1}(x),\,\mbox{\,where\,}\,
x=\frac{2Z_{nl}r}n.
\label{enl}
\end{equation}
Hence (\ref{func-nl}) becomes the following functional \cite{au}, \cite{un}
\begin{equation}
\ba{l}
\dst\widetilde\Phi(Z_{nl},A_{nl})=\frac1{2n}\int^{\infty}_0e^{-x}
\biggl[\frac{nx}{2Z_{nl}}U\left(\frac{nx}{2Z_{nl}}\right)-Z_{nl}+
A_{nl}\frac{nx}{2Z_{nl}}\biggr]^2xL_{nl}\,dx,
\ea
\label{f_nl}
\end{equation}
where $L_{nl}(x)=\frac{x^{2l+1}}{(n+l)!}(n-l-1)![L^{2l+1}_{n-l-1}(x)]^2$.
Using (\ref{frontera-nl}) the functional will depend only on $Z_{nl}$, then the minimum condition
($\frac{d\,\widetilde\Phi(Z_{nl})}{d\,Z_{nl}}=0$) leads us to the equation
$$
\dst aZ_{nl}^2-b(Z_{nl})Z_{nl}+c(Z_{nl})=0;
$$
where
$$\ba{c}\dst a=\frac{n^2-l(l+1)+1}n,\hspace{.3cm}
\dst b(Z_{nl})=\frac{2n^2-2l(l+1)+1}n(\alpha_0+\beta_0)+(\alpha_1+\beta_1)-
\frac1{2n}(\alpha_2+\beta_2),\\
\\
\dst c(Z_{nl})=\gamma+\frac1n[5n^2+1-3l(l+1)]\alpha_0\beta_0-
\frac1{2n}(\alpha_0\beta_2+\alpha_2\beta_0),\quad
\alpha_k=\int^{\infty}_0e^{-x}x^kL_{nl}(x)[rU(r)]\,dx
\ea
$$
with $k=0,\ 1,\ 2$, and $r=\frac{nx}{2Z_{nl}}$,
$$\ba{c}
\dst \beta_k=\int^{\infty}_0e^{-x}x^kL_{nl}(x)
\left[-r\,\frac d{dr}(rU(r))\right]\,dx,\quad
\gamma=\int^{\infty}_0e^{-x}x^kL_{nl}(x)
\left[-r\,\frac d{dr}(rU(r))rU(r)\right]\,dx.
\ea
$$
Taking into account that
$$
\dst \int^{\infty}_0e^{-x}x^kL_{nl}(x)\,dx=2^kZ_{nl}^{k-1}n^{2-k}\int^{\infty}_0
r^{k-1}\tilde R^2_{nl}(r)\,dr,
$$
the following equalities
$$
\ba{c}
\dst \int^{\infty}_0e^{-x}\frac1xL_{nl}(x)\,dx=\frac1{2l+1},\quad
\int^{\infty}_0e^{-x}L_{nl}(x)\,dx=1,\quad
\int^{\infty}_0e^{-x}xL_{nl}(x)\,dx=2n,\\
\\
\dst\int^{\infty}_0e^{-x}x^2L_{nl}(x)\,dx=2\left[3n^2-l(l+1)\right],\,
\int^{\infty}_0e^{-x}x^3L_{nl}(x)\,dx=4n\left[5n^2+1-3l(l+1)\right].
\ea
$$
hold.\\
Thus our strategy is quite simple: To find $Z_{nl}$ using the iterative method for
the following second order equation
$$ 
  a{Z_{nl}^{(s+1)}}^2-b(Z_{nl}^{(s)})Z_{nl}^{(s+1)}+c(Z_{nl}^{(s)})=0,\quad
  $$
in the form
$$
  Z_{nl}^{(s+1)}=\frac1{2a}\left[b(Z_{nl}^{(s)})+\sqrt{b^2(Z_{nl}^{(s)})-
4ac(Z_{nl}^{(s)})}\,\right],
$$
putting
$$
\dst Z^{(0)}_{10}=Z,\quad Z^{(0)}_{n0}=Z_{n-1\ n-2},\quad
Z^{(0)}_{nl}=Z_{n\ l-1}\quad(l\not=0).
$$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

In Figure 1 we provide the scheme of the used algorithm.

\vskip 1cm

\begin{figure}[h]
{\footnotesize %\scriptsize
\begin{picture}(100,450)(20,20)
%\begin{picture}(100,480)(20,20)
%\put(0,510){{\large\bf 3~\quad Diagram of calculation.}}
\put(140,20){\line(0,1){20}}\put(140,20){\line(1,0){40}}
	\put(180,40){\line(0,-1){20}}\put(180,40){\line(-1,0){40}}
\put(95,60){\line(1,1){20}}\put(95,60){\line(1,0){100}}
	\put(215,80){\line(-1,-1){20}}\put(215,80){\line(-1,0){100}}
\put(55,140){\line(0,1){20}}\put(55,140){\line(1,0){70}}
	\put(125,160){\line(0,-1){20}}\put(125,160){\line(-1,0){70}}
	\put(59,140){\line(0,1){20}}\put(121,160){\line(0,-1){20}}
\put(60,200){\line(1,2){5}}\put(60,200){\line(1,-2){5}}
	\put(120,200){\line(-1,2){5}}\put(120,200){\line(-1,-2){5}}
	\put(65,190){\line(1,0){50}}\put(65,210){\line(1,0){50}}
\put(60,240){\line(0,1){20}}\put(60,240){\line(1,0){60}}
	\put(120,260){\line(0,-1){20}}\put(120,260){\line(-1,0){60}}
	\put(64,240){\line(0,1){20}}\put(116,260){\line(0,-1){20}}
\put(60,300){\line(1,2){10}}\put(60,300){\line(1,-2){10}}
	\put(120,300){\line(-1,2){10}}\put(120,300){\line(-1,-2){10}}
	\put(70,280){\line(1,0){40}}\put(70,320){\line(1,0){40}}
\put(60,340){\line(0,1){20}}\put(60,340){\line(1,0){60}}
	\put(120,360){\line(0,-1){20}}\put(120,360){\line(-1,0){60}}
	\put(64,340){\line(0,1){20}}\put(116,360){\line(0,-1){20}}
\put(60,380){\line(0,1){20}}\put(60,380){\line(1,0){60}}
	\put(120,400){\line(0,-1){20}}\put(120,400){\line(-1,0){60}}
	\put(64,380){\line(0,1){20}}\put(116,400){\line(0,-1){20}}
\put(60,420){\line(0,1){20}}\put(60,420){\line(1,0){60}}
	\put(120,440){\line(0,-1){20}}\put(120,440){\line(-1,0){60}}
	\put(64,420){\line(0,1){20}}\put(116,440){\line(0,-1){20}}
\put(60,470){\line(1,1){20}}\put(60,470){\line(1,0){60}}
	\put(140,490){\line(-1,-1){20}}\put(140,490){\line(-1,0){60}}
\put(90,470){\vector(0,-1){30}}
\put(90,420){\vector(0,-1){20}}
\put(90,380){\vector(0,-1){20}}
\put(90,340){\line(0,-1){20}}
\put(90,280){\vector(0,-1){20}}
\put(90,240){\vector(0,-1){30}}
\put(90,190){\vector(0,-1){30}}
\put(90,140){\line(0,-1){20}}
\put(90,120){\line(-1,0){50}}
\put(40,120){\line(0,1){80}}
\put(40,200){\vector(1,0){20}}
\put(120,200){\line(1,0){20}}
\put(140,200){\line(0,-1){78}}
\put(140,118){\line(0,-1){18}}
\put(140,100){\line(-1,0){120}}
\put(20,100){\line(0,1){200}}
\put(20,300){\vector(1,0){40}}
\put(120,300){\line(1,0){40}}
\put(160,300){\vector(0,-1){220}}
\put(160,60){\vector(0,-1){20}}
\put(160,120){\vector(-1,0){70}}
\put(92,481){Input}
\put(71,473){$Z,\,A,\,\theta,\,\rho$}
\put(76,427){PARAM}
\put(72,387){VTHOM}
\put(78,347){POT(1)}
\put(69,305){$n=1,n_{\rm max}$}
\put(70,290){$l=0,h-1$}
\put(79,251){aprel $\phi$}
\put(70,244){$(n,\,l,\,za,\,e)$}
\put(80,198){j=ja,1}
\put(80,151){aprel $r$}
\put(65,144){$(n,\,l,\,j,\,zr,\,\varepsilon r)$}
\put(120,71){results {\tt output.res}}
\put(105,63){$n,l,j,Z_{nl},\varepsilon_{nl},Z_{nlj},\varepsilon_{nlj}$}
\put(154,27){end}
\put(225,485){File {\tt input.dat}}
\put(225,477){$Z$ --- nuclear charge}
\put(225,469){$A$ --- atomic weight}
\put(225,461){$\theta$ --- temperature in Kev}
\put(225,453){$\rho$ ---density in g/$cm^{3}$}
\put(225,445){{\tt out.res} --- file of results}
\put(225,430){parametrs of}
\put(225,420){calculation}
\put(225,390){calculation of Thomas- Fermi}
\put(225,380){potential for fixed}
\put(225,370){values of $\theta$ and $\rho$}
\put(225,350){approximation of Amaldi}
\put(225,340){for potential}
\put(225,310){$n_{max}$ --- maximum value of principal}
\put(225,300){quantum number n}
%\put(225,290){(max number of shell)}
\put(225,270){(nonrelativistic method)}
\put(225,250){resolution of the Schr\"odinger equation} 
\put(225,240){by nonrelativistic method}
\put(225,230){out put $Z_{nl}$ and $\varepsilon_{nl}$}
\put(225,200){$j+l=\widetilde j+\frac12$, $\widetilde j$ --- �����}
\put(225,190){inner quantum number j}
\put(225,180){$ja=0$, if $l=0$, then $ja=1$}
\put(225,150){solution of relativistic}
\put(225,140){method}
\put(225,130){$Z_r\leftarrow Z_{nlj}\quad
\varepsilon_r\leftarrow\varepsilon_{nlj}$}
\put(205,485){$\longleftarrow$}
\put(205,430){$\longleftarrow$}
\put(205,390){$\longleftarrow$}
\put(205,350){$\longleftarrow$}
\put(205,310){$\longleftarrow$}
\put(205,270){$\longleftarrow$}
\put(205,250){$\longleftarrow$}
\put(205,200){$\longleftarrow$}
\put(205,150){$\longleftarrow$}
\put(100,325){$Zn^{(0)}=Z$}
\put(100,265){$Za^{(0)}=Zb$}
\put(100,165){$Zr=Za$}
\put(75,105){$Zb=Za^{(0)}$}
\end{picture}
}
\caption{Algorithm for computation of effective charge by relativistic method.}
\end{figure}

\section{Summary of numerical results and conclusions.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5


\begin{table} 
{\scriptsize
\renewcommand{\arraystretch}{.3}
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}
\hline
& & & & & & & & & &\\
$n l j$  & $Z_{nl}P$ & $Z_{nl}S$ & $Z_{nlj}D$ & $Z_{nlj}J$ & $r_{\rm max}S$ & 
$r_{\rm max}D$ & $\varepsilon_{nl}P$ & $\varepsilon_{nl}S$ &
$\varepsilon_{nlj}D$&$\varepsilon_{nlj}J$  \\
& & & & & & & & & & \\\hline
& & & & & & & & & & \\\hline
& & & & & & & & & &\\
 $1 0 \frac12$& $12.21$& $12.50$& $12.50$& $12.21$& $8.02\cdot 10^{-2}$& $7.99\cdot 10^{-2}$& $6.37\cdot 10^1$& 
 $6.37\cdot 10^1$& $6.39\cdot 10^1$& $6.38\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $2 0 \frac12$& $11.18$& $11.39$& $11.40$& $11.18$& $4.60\cdot 10^{-1}$& $4.59\cdot 10^{-1}$& $1.14\cdot 10^1$& 
 $1.14\cdot 10^1$& $1.14\cdot 10^1$& $1.14\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $2 1 \frac32$& $10.44$& $10.70$& $10.70$& $10.44$& $3.74\cdot 10^{-1}$& $3.74\cdot 10^{-1}$& $1.04\cdot 10^1$& 
 $1.04\cdot 10^1$& $1.04\cdot 10^1$& $1.04\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $3 2 \frac32$& $9.70$& $9.81$& $9.81$& $9.70$& $9.17\cdot 10^{-1}$& $9.17\cdot 10^{-1}$& $3.05\cdot 10^0 $& 
 $3.05\cdot 10^0 $& $3.05\cdot 10^0$& $3.05\cdot 10^0$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $3 2 \frac52$& $9.70$& $9.81$& $9.81$& $9.70$& $9.17\cdot 10^{-1}$& $9.17\cdot 10^{-1}$& $3.05\cdot 10^0$& 
 $3.05\cdot 10^0$& $3.05\cdot 10^0$& $3.05\cdot 10^0$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $4 0 \frac12$& $10.09$& $10.12$& $10.12$& $10.09$& $2.43\cdot 10^0$& $2.43\cdot 10^0$& $1.10\cdot 10^0$& 
 $1.09\cdot 10^0$& $1.09\cdot 10^0$& $1.10\cdot 10^0$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $4 2 \frac52$& $9.37$& $9.42$& $9.42$& $9.37$& $2.25\cdot 10^0$& $2.25\cdot 10^0$& $8.82\cdot 10^{-1}$& 
 $8.84\cdot 10^{-1}$& $8.84\cdot 10^{-1}$& $8.82\cdot 10^{-1}$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $4 3 \frac52$& $9.12$& $9.30$& $9.30$& $9.12$& $1.72\cdot 10^0$& $1.72\cdot 10^0$& $8.12\cdot 10^{-1}$& 
 $8.12\cdot 10^{-1}$& $8.13\cdot 10^{-1}$& $8.11\cdot 10^{-1}$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $4 3 \frac72$& $9.12$& $9.30$& $9.30$& $9.12$& $1.72\cdot 10^0$& $1.72\cdot 10^0$& $8.12\cdot 10^{-1}$& 
 $8.12\cdot 10^{-1}$& $8.13\cdot 10^{-1}$& $8.11\cdot 10^{-1}$\\
& & & & & & & & & &\\\hline
\end{tabular} }
\caption{$Z$ = 13. $\theta$=$1.0000\cdot 10^{-1}$ $kev$ and $\rho$ = $2.7000\cdot 10^{-1}$ $g/cm^{3}$ given 
Thomas-Fermi potential, as initial approximation. The energy values and $r_{max}$ are given in 
atomic units.}
%%\end{table} 

\vskip 1cm 

%%%\input rue2.tex

%%% \begin{table}
{\scriptsize
\renewcommand{\arraystretch}{.3}
\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}
\hline
& & & & & & & & & &\\
$n l j$  & $Z_{nl}P$ & $Z_{nl}S$ & $Z_{nlj}D$ & $Z_{nlj}J$ & $r_{\rm max}S$ & 
$r_{\rm max}D$ & $\varepsilon_{nl}P$ & $\varepsilon_{nl}S$ &
$\varepsilon_{nlj}D$&$\varepsilon_{nlj}J$  \\
& & & & & & & & & & \\\hline
& & & & & & & & & & \\\hline
& & & & & & & & & &\\
$1 0 \frac12$& $77.89$& $78.36$& $78.48$& $78.12$& $1.28\cdot 10^{-2}$& $1.05\cdot 10^{-2}$& 
$2.77\cdot 10^3$& $2.77\cdot 10^3$& $3.08\cdot 10^3$& $3.08\cdot 10^3$ \\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\
$2 0 \frac12$& $74.49$& $75.49$& $76.42$& $75.64$& $6.94\cdot 10^{-2}$& $6.02\cdot 10^{-2}$& 
$5.53\cdot 10^2$& $5.53\cdot 10^2$& $6.37\cdot 10^2$& $6.37\cdot 10^2$ \\
& & & & & & & & & & \\\hline
& & & & & & & & & &\\
$2 1 \frac12$& $72.82$& $74.38$& $75.67$& $74.27$& $5.38\cdot 10^{-2}$& $4.48\cdot 10^{-2}$& 
$5.35\cdot 10^2$& $5.36\cdot 10^2$& $6.18\cdot 10^2$& $6.17\cdot 10^2$\\

& & & & & & & & & &\\\hline
& & & & & & & & & &\\\hline
& & & & & & & & & &\\
$4 3 \frac72$& $58.39$& $59.25$& $59.31$& $58.43$& $2.70\cdot 10^{-1}$& $2.68\cdot 10^{-1}$& 
$8.53\cdot 10^1$& $8.54\cdot 10^1$& $8.57\cdot 10^1$& $8.57\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\
$5 0 \frac12$& $63.64$& $63.65$& $65.13$& $65.15$& $6.27\cdot\cdot 10^{-1}$& $5.86\cdot 10^{-1}$& 
$5.57\cdot 10^1$& $5.55\cdot 10^1$& $6.00\cdot 10^1$& $6.03\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\
$5 1 \frac12$& $62.12$& $62.28$& $63.82$& $63.42$& $6.25\cdot 10^{-1}$& $5.83\cdot 10^{-1}$&
 $5.39\cdot 10^1$& $5.39\cdot 10^1$& $5.82\cdot 10^1$& $5.80\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $6 4 \frac92$& $55.81$& $55.84$& $55.86$& $55.82$& $8.37\cdot 10^{-1}$& $8.35\cdot 10^{-1}$& 
 $2.87\cdot 10^1$& $2.87\cdot 10^1$& $2.88\cdot 10^1$& $2.88\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $6 5 \frac92$& $55.24$& $55.43$& $55.44$& $55.25$& $6.50\cdot 10^{-1}$& $6.47\cdot 10^{-1}$& 
 $2.84\cdot 10^1$& $2.84\cdot 10^1$& $2.85\cdot 10^1$& $2.85\cdot 10^1$\\
& & & & & & & & & &\\\hline
& & & & & & & & & &\\ 
 $6 5 \frac{11}2$& $55.24$& $55.43$& $55.43$& $55.25$& $6.50\cdot 10^{-1}$& $6.48\cdot 10^{-1}$& 
 $2.84\cdot 10^1$& $2.84\cdot 10^1$& $2.84\cdot 10^1$& $2.84\cdot 10^1$\\
& & & & & & & & & &\\\hline
\end{tabular} }
\caption{$Z$ = 79. $\theta$=$1000\cdot 10^{-1}$ $kev$ and $\rho$ = 1.9300  $g/cm^{3}$ 
given  Thomas-Fermi potential, as first approximation. The energy values and $r_{max}$ are given in 
atomic units.}

\vskip .6cm 

In the next tables are shown the following result corresponding to the gold ($Au$) and aluminium ($Al$) for fixed 
$\rho$ and $\theta$ and for some combination of quantum numbers $nlj$, although we computated for each 
combination of quantum numbers:  

\vskip .6cm 

$Z_{nlj}J$ and $\varepsilon_{nlj}J$ --- effective charge and energy levels found from our model.\\
\esp $Z_{nl}S$ and $\varepsilon_{nl}S$ --- effective charge and energy levels corresponding to the 
numerical solutions for Schr\"odinger equation.\\
\esp $Z_{nlj}D$ and $\varepsilon_{nlj}D$ --- effective charge and energy levels corresponding to the 
numerical solutions for Dirac equation.\\
\esp $Z_{nl}P$ and $\varepsilon_{nl}P$ --- effective charge and energy levels corresponding to the 
nonrelativistic cases.\\
\esp $r_{max}S$ --- position of absolute maximum of wave function corresponding to the numerical 
solution of Schr\"odinger 
equation.\\
\esp $r_{max}D$ --- position of absolute maximum of wave function corresponding to the numerical 
solution of Dirac equation.


\end{table} 



To check the quantum-statistical model describing the state of matter in the region of
high temperatures, two elements, one of low-Z (aluminium) and other of high-Z (gold), of interest
in several physical applications have been employed. The range on temperature and density,
for which the model is applicable, covers the region $10^{-3} g/cm^{3}\leq\rho\leq10^{2} g/cm^{3}$, 
$5 ev\leq\theta\leq15 Kev$, depending on the substance and quantum number combinations. 
However, in this work fixed representative temperatures and densities has been selected. Results 
from the model are compared with other numerical results. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%\input rue1.tex

\normalsize

Notice the small difference of results compared with several methods. In general, the 
tendency for high principal quantum number $n$ is the convergence between all methods.\\
A good agreement between all radial wave functions calculating by several methods for gold are 
illustrated in the next figures. The large componets are quite similar in all cases. Here our 
method has called Z-method. It is important to
remark that the graphic 3 gives results comparing calculated the small 
components for Z-method with other one and large radial component of the numerical solution of the
Dirac equation. Therefore, we can neglect the influence of the small components to 
generate opacity calculation code. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[t]
\vspace*{3cm}
 \begin{center}
  \setlength{\unitlength}{1cm} 
 \hspace{-1cm}
   \begin{picture}(12,9)
    \epsfxsize=13cm
    \epsffile{gra.eps}
   \end{picture}
 \end{center}
  \vspace*{-1cm}
 \caption{Comparison of numerical results. Here $r_{0}$ is the  mean value of atomic radius ($r_{0}=1.388
(\frac{A}{\rho})^{\frac{1}{3}}$).}
\label{fig1}
\end{figure}



In summary, we demostrated that simple formulas for atomic wave functions and energy levels can 
be used to determine other properties of the matter in very short computational time for any 
substance at several densities and temperatures. The numerical results show small deviations 
from numerical solutions for high-Z elements, which is balanced by the simplicity and 
rapidity of the computations, and for low atoms is regarded complete coincidence of the wave 
functions, calculated with help of different methods. In this way the model can be used in the 
opacity and equation of state calculations scheme, because it has a wide flexibility depending 
on plasma conditions. 
 
\vskip .4cm 

\underline{ACKNOWLEDGEMENTS}: 
 
\vskip .4cm 

I want to express my special thanks to three people from M. V. Keldysh Institute of 
Applied Mathematics, Russian Academy of Sciences, 
who made this report possible. Professor A. F. Nikiforov introduced me in these areas, orthogonal
polynomials and mathematical physics. He was the advisor of my master thesis on which is based 
this work and also answered all my questions, and was always willing to explain the extreme difficult 
quantum-statistical models. Professor V. G. Novikov for his useful disscusion and remarks. 
Dr. A. D. Solomiannaya for her help with Fortran program. I also would like to thank the Professor 
F. Marcell\'an and Dr. R. Alvar\'ez-Nodarse for their suggestions, remarks and helping me 
prepare the english version.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      BIBLIOGRAFIA    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{10}

\bibitem{af} {\sf D. R. Alexander and J. W. Ferguson}: {\it Low-Temperature 
Rosseland Opacities}. {\sc ApJS.} {\bf 437} (1994), 879-891.

\bibitem{au} {\sf J. Arves\'u}: {\sc Attainment of relativistic hydrogenic-like
wave functions for self-consistent fields in high-temperature plasmas}. Master
thesis in mathematical physics. (Moscow State University, {\em M. V. Lomonosov}, 1995.)
[in Russian].

\bibitem{bs} {\sf H. A. Bethe, E. S. Salpeter}: {\sc Quantum Mechanics of 
One-and Two-Electron Atoms}. (Springer, Berlin, 1957.)

\bibitem{cms} {\sf T. R. Carson, D. F. Mayers and D. W. N. Stibbs}: {\it The 
calculation of Stellar Radiative Opacity}. {\sc Mon. Not. R. astr. Soc.} 
{\bf 140} (1968), 483-536.

\bibitem{dy} {\sf V. P. Dyakonov}: {\sc Handbook of algorithms and programs in basic
program for PC}. (Nauka, Moscow, 1987.)[in Russian].


\bibitem{ir} {\sf C. A. Iglesias and F. J. Rogers}: {\it Opacities for the 
Solar Radiative Interior}. {\sc ApJS.} {\bf 371} (1991), 408-417.

\bibitem{irw} {\sf C. A. Iglesias, F. J. Rogers and B. G. Wilson}: {\it Spin-Orbit 
Interaction effects on the Rosseland mean Opacity}. {\sc ApJS.} 
{\bf 397} (1992), 717-728.

\bibitem{lat} {\sf R. Latter}: {\it Temperature behavior of the Thomas-Fermi 
statistical model for atoms}. {\sc Phys. Rev.} {\bf 99} (6) (1955), 1854-1870.

\bibitem{nsu} {\sf A.F. Nikiforov, S. K. Suslov and  V. B. Uvarov}: 
{\sc Orthogonal Polynomials in Discrete Variables}. Springer Series in 
Computational Physics. (Springer Verlag, Berlin, 1991.) 
 
\bibitem{nu} {\sf A. F. Nikiforov and V. B. Uvarov}: {\sc Special Functions of 
Mathematical Physics}. (Birkh\"auser Verlag, Basel, 1988.)

\bibitem{rwi} {\sf F. J. Rogers, B. G. Wilson and C. A. Iglesias}: {\it Parametric 
potential method for generating atomic data}. {\sc Phys. Rev. A} {\bf 38} (10) (1988),
5007-5020.

\bibitem{un} {\sf V. B. Uvarov and A. F. Nikiforov}: {\it About one approximation 
method of solution for Schr\"odinger equation}. {\sc ZhVM i MF} {\bf 1} (1) (1961)
[in Russian].

\end{thebibliography}

\begin{flushleft}
{  Jorge Arves\'u } \\
Departamento de Matem\'aticas.\\
 Escuela Polit\'ecnica Superior.\\
Universidad Carlos III de Madrid. \\
c/ Butarque 15, 28911, Legan\'es, Madrid. Spain. \\
{\sf E-Mail}{\small \tt : arvesu@dulcinea.uc3m.es }
%%{\sf URL}{\small \tt: http:\backslash\backslash dulcinea.uc3m.es\backslash users\arvesu\arvesu.html}
\end{flushleft}

\newpage
$$ \ $$

\end{document}



