Recursive determination of phase shifts for screened Coulomb potentials

In the calculation of hot-plasma atomic structure, the continuum wavefunctions are characterized by phase shifts, which therefore determine the scattering cross-sections. In this short paper, we propose a recurrence relation for the phase shifts in the case of a particular type of parametric potentials widely used in atomic-structure codes. These potentials have to be linear combinations of static screened Coulomb potentials (Yukawa-type potentials) multiplied by polynomial functions.


Introduction
The potential felt by an electron in the atom is not coulombic at intermediate distances, and the difference between the effective potential and the Coulomb potential manifests itself through the quantum defect associated with bound states and the phase shift associated with continuum states. For a given angular momentum, the quantum defect g ℓ n [1] is defined, in atomic units, by where  ℓ n is the binding energy of the excited electron in the ℓ n orbital and Z the charge of the ionic core. When the electron energy is above the continuum limit, its wavefunction is characterized by a phase shift d ℓ for each orbital quantum number ℓ. The asymptotic solution of the radial Schrödinger equation where C ℓ is a constant,  = k 2 and + ℓ J 1 2 andℓ J 1 2 are Bessel functions of the first kind. When  ¥ r , y ℓ has the asymptotic form The difference between the actual and the coulombic phase shift is the signature of the non-coulombic part of the electron-ion interaction. Phase shifts are important for the photo-ionization cross-sections and a reliable calculation of the impact ionization cross-section of an atom requires an accurate determination of the continuum wavefunctions in the incident and in the exit channels. In the theory of collisions, phase shifts determine the scattering cross-sections [2,3]. For instance, the elastic scattering cross-section of particles is The phase shift plays a major role in phase-amplitude methods [4][5][6][7][8]; but whatever the technique chosen for its computation (solving Schrödinger equation or using semi-classical methods [9]), the difficult point is the determination of the potential V(r). Due to configuration mixing, the number of involved radial integrals can be very large, and analytical potentials can be an alternative to self-consistent field methods. Tannous et al [10] tried to determine potentials which incorporate the effect of exchange while keeping a local character. Parametric potentials are often used [11], and are the key ingredient of a number of atomic-structure codes, such as Hebrew University Lawrence Livermore Atomic Code [12][13][14][15], OPAL [16], super transition arrays [17][18][19] or flexible atomic code [20]. In this work, we propose, following Tietz [21,22], to evaluate the difference between phase shifts associated to consecutive values of ℓ by the relation which is obtained by replacing the wavefunction by its asymptotic expression (4). Assuming in addition that the difference d d Two parametric potentials used in widely used atomic-structure codes are presented in sections 2 and 3, and the new recursion relation is explained in section 4.

Klapisch's parametric potentials
The electronic shell structure in an atomic potential was introduced by Klapisch [23]. The parametric potential method consists in optimizing an analytical central-field potential according to a chosen quality criterion (variational, spectroscopic, etc). The method is particularly interesting, compared to Hartree-Fock for example, when the wavefunctions are nearly solutions of a central field (this is the case for alkali-like spectra or for highly ionized atoms) or if the spin-orbit is strong (rare gases, medium and heavy atoms) [24,25]. It was implemented in the RELAC code [26,27] and was widely used, for instance in the calculation of photo-ionization crosssections [28][29][30][31]. Let us consider a closed electronic subshell of q electrons, having spherical symmetry, described by the normalized radial density of charge: where α is a positive real constant (α plays the role of the inverse of a screening length). The potential of another electron at radius r in the field of the latter density of charge and a nucleus of charge Z reads [32] a a = - To each subshell i corresponds a parameter a i . The parameters a i , Î [ ] i N 1, s , N s being the number of subshells, are adjusted by a least-square procedure to reproduce experimental data or ab initio calculations. The function a ℓ ( ) f r , , of equation (10) can be expressed in terms of incomplete Gamma functions (see appendix A), which can be expressed in terms of confluent hypergeometric functions (or 'Kummer functions') F 1 1 . The first term in the summation ( j=0) corresponds to the Yukawa potential [33], for which it is difficult to obtain an analytical treatment of bound states, although significant progress was made in the last few years (see for instance [37]).
Historically, Yukawa showed in the 1930s that such a potential arises from the exchange of a massive scalar field such as the field of a massive boson. Since the field mediator is massive, the corresponding force has a certain range, which is inversely proportional to the mass of the mediator particle [34].
The screening (or shielding) effect describes the attraction between an electron and the nucleus in any atom with more than one electron. It can be defined as a reduction in the effective nuclear charge felt by an electron, due to its interactions with the other electrons and to the interactions of these electrons with the nucleus. In lowdensity, high-temperature plasmas, electric-field screening can be taken into account using the Debye-Hückel model [35].
In solid-state physics, the screened potential is used to calculate the electronic band structure of a large variety of materials, often in combination with pseudopotential models.

Parametric potential of Rogers et al
Rogers et al [16,38] defined the electron configuration as having two components: the first one is a parent configuration consisting of all the electrons in a given configuration except one. The excluded electron defines the second component or 'running' electron. The authors introduced a parametric potential for each parent configuration. This potential consists of a long-range Coulomb part and a screened function represented by a sum of static screened Coulomb (Yukawa-type) potentials [33]: is the number of electrons in the parent ion, N n the number of electrons in the shell with principal quantum number n, n * the maximum value of n for the parent configuration and a n the screening constant for electrons in shell n. For a given ion, it is possible to define as many parent configurations as occupied subshells. The screening parameters are obtained by an iterative procedure consisting in the resolution of spin-averaged Dirac equation and matching the eigenvalues to the ionization potentials. Such an approach enables one to account for effects (electron correlations, spin-other-orbit interactions) which are not necessarily included in ab initio calculations. The authors provided accurate fits of the screening parameters along an iso-electronic sequence. The estimation of the screening constants was improved by Mabong et al [39] in order to include relativistic effects and later by Mendoza et al [40] and others. Rogers' potential encountered a great success for the generation of atomic quantities (levels, oscillator strengths, photo-ionization cross-sections, etc) required in astrophysics. The potential can be applied to any excitation and ionization state. Moreover, since its Fourier transform is simple, it is a good candidate for the determination of cross-sections of elementary processes.
In the present work we only consider the static screening effects. It is however important to keep in mind that the dynamic screening effects are also important not only in weakly coupled plasmas but also in strongly coupled plasmas such as a dusty plasmas [41][42][43]. Many interesting studies were (and are being) carried out on the subject (time-dependent density functional theory, linear response theory, etc), but including the effect of dynamic screening properly in our model would be a difficult task and is beyond the scope of the present article. In the framework of the linear response theory, Zangwill and Soven [44] applied the random-phase-approximation to calculate the photo-absorption cross-section of rare gases. This was the starting point for the development of the time-dependent density functional theory. Extending the scope of this work to plasmas leads to some difficulties because the atoms cannot be considered without delocalized electrons. The latter may contribute to particlehole transitions but may also be responsible for collective effects. A proper description of dynamic screening would then require a full quantum-mechanical description of both bound and free electronic states. The effect of dynamic screening on the energies of the electronic subshells would be difficult to quantify but we expect it to be rather small. Dynamical screening is known to have an impact on the radiative spectra (see for instance [45], through channel mixing (photo-excitation/photo-ionization) and configuration interaction.

Recursion relation for the phase shift
The recursion evaluation of the phase shift following the requirements of Tietz [21] (see equation (6)) involves only integrals of the type where m J and n J are Bessel functions of the first kind. Although efficient algorithms do exist (for instance Lucas [46] proposed a method which makes use of extrapolation on a sequence of partial sums, and requires rewriting the product of Bessel functions as the sum of two more well-behaved functions), the evaluation of infinite integrals involving products of Bessel functions is tedious. The integral (13) can be expressed, if l m n + + > 0 and a > 0, as (see equation where Γ is the usual Gamma function and F p q is a generalized hypergeometric function defined by In particular, F 1 1 is the confluent hypergeometric function, which includes Bessel functions, Laguerre polynomials, incomplete Gamma functions, etc as special cases. The function F 2 1 , commonly referred to as the Gauss hypergeometric series, includes Chebyshev, Legendre and Jacobi polynomials. In the particular case where b g = = k, we find, after some algebraic manipulation, the following identity:  This result was also obtained by Mahajan [48]. A similar integral, involved in the difference between d -( ) The numerical evaluation of hypergeometric functions F p q was investigated by many authors. For instance, the program of Perger et al [49] relies on a direct numerical evaluation of the series. The only transformation formula employed is the use of the linear transformation when | | z is approaching 1 from below. Recently, Willis [50] proposed an acceleration procedure through precise remainder asymptotics. He expressed the asymptotics of the remainders of the partial sums of the generalized hypergeometric function F p q through an inverse power series where the exponent λ and the asymptotic coefficients c k may be recursively computed to any desired order from the hypergeometric parameters and argument. The hypergeometric series F 3 2 can also be evaluated numerically using the following integral representation [51,52]: with  d a 0 and  e b 0. Such a representation is interesting since the integral over the whole space (from 0 to ¥) of the left-hand side of equation (18) is replaced by a finite-range integral between 0 and 1. Other numerical techniques exist. For example, Wills [53] and a few years later Bretz [54], arranged the series expansion into a nested form. The ease of computing numerically the F p q using Horner's rule [55] for polynomial evaluation has been found to yield fortran programs which are very accurate (see for instance the program by Srinivasa Rao in [56]). The generalized hypergeometric function is evaluated as a a However, we do not pretend that the numerical evaluation of hypergeometric functions is more efficient from a numerical point of view (speed, accuracy) than direct integration. Such considerations depend on the numerical integration technique used, on the machine, etc. The purpose of the present work is to provide analytical expressions, which interest, beyond the purely mathematical aspect, relies in the fact that they can open the way to algebraic manipulations, through the use of properties, recurrence relations, etc of hypergeometric functions.
Let us consider the static screened Coulomb potential 2 7 r for which the calculation of the phase shifts was investigated by many authors. Rogers [57] proposed a direct integration near the origin and firstand second-order WKB approximations at larger distance. Grandjouan and Deutsch [58] used Numerov method and variable-phase approach, applied to the Sturm-Liouville form of the Schrödinger equation. Bechler and Pratt [59] performed the calculation in the framework of perturbation theory. Following the prescription of Tietz [21], we get In order to test the validity of the approximation, let us take the following potential (see figure 1) from [60]: with the parameters given in table 1 for the case of argon, in table 3 for mercury and in table 5 for uranium. The integral (18) must be evaluated for the three different pairs g l ( ) , i i :   Table 1. Numerical values of the potential parameters g i and l i for argon (from [60] for the Potential given by equation (30)).   Table 3. Numerical values of the potential parameters g i and l i for mercury (from [61,62] with the potential given by equation (30)). Num. int. This work Num. int. .

2
The numerical integration technique is described in [63]; it consists in a modification of the Gauss-Jackson method [64]. For mercury and uranium, the relativistic effects cannot be neglected. Therefore, due to the high values of the atomic numbers (80 and 92 respectively), a proper description of the atomic structure would require solving Dirac's equation. However, in that case, a generalization of the present approach would be more complicated (in particular because the wavefunction has two components). When relativity is taken into account in the calculation of atomic structure, the energies of s and p orbitals become lower than in the non-relativistic case, on the contrary to d and f orbitals. The spin-orbit splitting is a very important relativistic effect. It concerns only non-s orbitals. For s orbitals the mass-velocity correction dominates over the Darwin correction, while for non-s orbitals the Darwin correction is small. For inner-core s and p orbitals, the contraction of the orbital is easily explained by the relativistically increased mass. For instance, the energy of s 1 orbital of uranium is (in atomic units): −3690.78 in the non-relativistic (Schrödinger) case, −4255.56 in the relativistic (Dirac) case, and −4114.71 in the semi-relativistic Pauli approximation. However, in the present work, we are dealing with continuum (free-electron) wavefunctions which are much less sensitive to relativistic effects than the bound states. Moreover, the relativistic effects become less and less important as the energy increases, and all our examples correspond to very high values of continuum energy. We are presenting an approximate method, in which the wavefunctions are represented by Bessel functions and we are only interested in phase shifts. The potential error resulting from the fact that we neglect relativistic effects is much smaller than the one resulting from the approximations mentioned above.
Using the function 'Timing' of the software Mathematica [65], we find that our approach (see formula (31)) is much more efficient than the direct numerical integration of the expression given in equation (7). Different time evaluations are given in table 7 for different values of atomic number Z and orbital quantum number ℓ. As expected, the numerical cost of the direct integration increases with ℓ.
As explained above, Tietz approximation relies on two assumptions: the wavefunctions are replaced by their approximate form and the phase shift differences must be small. The later approximation is not necessary; indeed, instead of equation (6) we can take  Table 5. Numerical values of the potential parameters g i and l i for uranium (from [63] for the potential given by equation (30)).

Conclusion
Starting from Tietz formula, in which the wavefunctions are replaced by their asymptotic form, we proposed an analytic expression for the difference between two consecutive phase shifts (in terms of orbital quantum number ℓ) for a class of screened Coulomb potentials. This recursion relation, which involves a F 3 2 generalized hypergeometric function, enables one to obtain approximate phase shift in atomic structure codes using parametric potentials with respect to k using equation (31) and the potential used in [60] (the parameters are provided in table 1). for several values of Z and ℓ: comparison between the present approach (equation (31)) and the numerical integration (equation (7)). consisting of the summation of a pure Coulomb part and static screened Coulomb potentials multiplied by polynomial functions.

Acknowledgments
The author would like to thank D Teychenné and A Decoster for useful discussions about the numerical evaluation of the generalized hypergeometric functions.