Multi-Conﬁguration Calculation of Ionization Potential Depression

: The modelling of ionization potential depression in warm and hot dense plasmas constitutes a real theoretical challenge due to ionic coupling and electron degeneracy effects. In this work, we present a quantum statistical model based on a multi-conﬁguration description of the electronic structure in the framework of Density Functional Theory. We discuss different conceptual issues inherent to the deﬁnition of ionization potential depression and compare our results with the famous and widely-used Ecker-Kröll and Stewart-Pyatt models.


Introduction
Ionization potentials are fundamental quantities in atomic physics. Their experimental and theoretical studies have aroused the interest of scientists for many years, either for accurate determination [1] or looking for regularities along isoelectronic sequences [2], etc. Ionization potential depression (IPD) or continuum lowering in dense plasmas has been the subject of many investigations over the past years [3][4][5][6]. This many-body effect significantly alters the ionization balance and the corresponding charge state distribution, which have a strong influence on thermodynamic [7], transport [8][9][10] and radiative properties of the system [11,12]. Such a topic is of crucial importance for astrophysics (stellar physics, planetary science, ...) [13][14][15][16], solid state physics and inertial confinement fusion [17,18]. It has an impact on the cross-sections of the different atomic processes, such as collisional excitation or ionization [19,20]. However, an accurate description of IPD in an interactive many-body environment is extremely difficult [21][22][23][24][25][26], because strong correlations and quantum degeneracy have to be taken into account in a consistent manner. Different semiempirical models have been proposed for the IPD. In particular, the Ecker-Kröll (EK) [27] and the Stewart-Pyatt (SP) model [28] have been widely applied in plasma physics for simulating physical properties and the analysis of experimental measurements. Thanks to the recently developed experimental facilities, it has become possible to explore warm/hot and high-energy density matter, whereby the plasmas are found to be strongly coupled and nearly degenerate. The experimental outcomes can be used to benchmark physical assumptions and theoretical models [19,20]. Over the last few years, new experiments related to the IPD have been performed with high-intense laser beams at LCLS (Linac Coherent Light Source) [29] at SLAC National Accelerator Laboratory [30,31], at ORION in the United Kingdom [32], and at NIF (National Ignition Facility) at Lawrence Livermore National Laboratory [33,34]. Such experiments have revealed the lack of a consistent understanding of IPD, since no satisfying explanation can be drawn for the corresponding measurements using the commonly accepted IPD models (EK and SP). Therefore, a deeper investigation of the IPD in warm/hot dense regime is required. Attempts for a better understanding of the recent experimental data have been made using different models and simulation methods. In particular, two-step Hartree-Fock calculations [35], simulations based on the finite-temperature density functional theory [36,37], classical molecular dynamics simulations [38], and Monte-Carlo simulations [39] have been worked out. Although such outstanding approaches are very successful to explain properties of warm dense matter, they are restricted due to their high computational cost. Besides these simulation methods, other analytic improvements have also been proposed, such as fluctuation model [40], analytic self-consistent approach valid in a wide range of temperatures and densities [41] and atomic-solid-plasma model [42]. The problem of line-shifts in LTE (Local Thermodynamic Equilibrium) plasmas was addressed by Alexiou et al. in the framework of static ions and "active spectator" model. The authors showed that lines can have unshifted (or slightly shifted) peaks with asymmetric tails if the free electron wave functions are strongly localized by the random field of the plasma or by dynamic effects [43] (see in particular the subsection 3.2 entitled "The last bound level"). To emphasize the likely importance of dynamic effects, it is worth quoting the last paragraph of Ref. [44]: "Reading Dick Lee's comments on the first draft of this paper I realized that, although the work deals with behavior of electrons in a static plasma, the most powerful localization mechanism is dynamic: the finite life-time of a "free" electron. A bound electron becomes "free" when it is kicked out of an ion by a photon, electron, or a collective process. At that time its wave function is still localized at, approximately, the same volume which it had occupied before. Then the wave function begins its time-dependent spread out, and, eventually would have become macroscopic, but, before that happens, the electron is recaptured by an ion, and ends up in another (local) bound state. So, between the time it was "kicked out" and the time it was recaptured it had only a finite period for spreading. This mechanism is unique to a hot matter, where there are enough holes in which the electron can be recaptured". Very recently, Ren et al. explored the evolution of systems driven to high energy densities using a non-LTE, Fokker-Planck collisional-radiative code. The authors investigated the evolution dynamics of a solid-density plasma driven by an XFEL, and explored the relaxation of the plasma to local thermodynamic equilibrium on femtosecond timescales in terms of the charge state distribution, electron density, and temperature [45]. Motivated by the theoretical as well as experimental challenges related to the understanding of IPD, we present, starting from a quantum statistical theory, a self-consistent multi-configuration approach to IPD. The main idea consists in grouping together (into a so-called superconfiguration) all the configurations representing the excited states of an ion. Therefore, each ion is represented by one single superconfiguration. Then, a self-consistent calculation is carried out for each superconfiguration within the Density Functional Theory (DFT) in the Local Density Approximation (LDA). In such a way, the energies of the subshells and therefore the total energies of the ions, take into account the so-called "density effects" or plasma environment effects. The IPD is obtained as the ionization energy of a given ion from which the "unscreened" value, corresponding to the isolated ion, is subtracted.
Our DFT approach is different from the one used in Ref. [36]. In the latter work, the authors have employed the finite-temperature plane-wave DFT code ABINIT to calculate electron states and energies within a plasma, combined with the projector augmented-wave (PAW) potential formalism used to model a range of different aluminum ions, constructed by creating excited-electron configurations in inner atomic shells and freezing them in the core of the potential. Such a technique belongs to the so-called "ab initio" or moleculardynamics methods. Our approach is inspired from a different class of models, often referred to as "average-atom models". The main difference between our approach and usual average-atom models is that we impose an integer number of electrons, in order to deal with a specific ion, in the superconfiguration formalism (the average-atom model provides a mean electronic structure, and therefore fractional numbers of bound electrons, in the different subshells and for the total).
Ecker-Kröll and Stewart-Pyatt analytical models are recalled in Section 2. The selfconsistent calculation of electronic structure including screening is presented in Section 3, and the unscreened case in Section 4. The calculation of Ionization Potential Depression is described in Section 5. The results are analyzed in Section 6 and compared to an alternative theoretical approach proposed by Massacrier [46][47][48]. The influence of the temperature on the isolated-atom energy values is discussed, together with the validity of Rayleigh-Schrödinger perturbation theory.

Analytical Ionization Potential Depression Models
Let us consider an atom (atomic number Z) which can be partially or totally ionized and k the ionization state (0 ≤ k ≤ Z), having a fraction F k . The ionic density n i is given by where N A is the Avogadro number, ρ the matter density (in g/cm 3 ) and A the molar atomic mass (in g). The ionic populations are then N k = F k .n i . We define were n e represents the electron density. Let us introduce (setting the Boltzmann constant

In particular
is the Wigner-Seitz radius (i.e., the radius of the average spherical cell representing the average ion) and R(n i + n e ) = 3 4π(n i + n e )

Debye-Hückel Model
The expression of the IPD in the Debye-Hückel (DH) theory [49] is where i is the charge state and z = i + 1.

Stewart-Pyatt Model
The IPD has been formulated in the sixties following two different ways. First, Stewart and Pyatt proposed a model using the finite-temperature potential for the average electrostatic potential near nuclei of the plasma particles [28]. The bound electron are considered as part of the unperturbed ion. The plasma free electrons are described by Fermi-Dirac statistics and ions by Maxwell-Boltzmann statistics. In such a model, the bound electrons do not contribute to the reduction of the ionization energy. The Stewart-Pyatt ionization potential depression is given by [50]: There are some differences between the understanding and implementation of Stewart-Pyatt model in the literature. Such a point is clarified in Appendix A.

Ecker-Kröll Model
Second, Ecker and Kröll formulated a generalized Saha equation as a function of the chemical potential of the plasma [27]. This model assumes two different forms for the IPD depending on the particle (ions + electrons) density. Let us introduce the critical density the radius R EK = R(n i + n e ) and the constant If n i + n e ≤ n c , we take (i still represents the charge state): and if n i + n e ≤ n c , the EK IPD is provided by and the continuity between Equations (1) and (2) is ensured at n c . The dependence of C EK on atomic number Z is discussed in Appendix B.

Modified Ecker-Kröll Model
Recent experimental results have shown discrepancies with the SP model which is the most widely used among the IPD models and shown good agreement with the EK model version in which C EK has been set to 1 (referred to as "modified EK" model) according to experimental considerations [30]. In contrast, other experimental results obtained independently have corroborated the SP model [32] (see Section 6.1). In the following, we will compare our simulation results with both models SP and EK with C EK = 1.

Self-Consistent Calculations-General Case (Screening of Ions)
A superconfiguration [51,52] is an ensemble of configurations. The superconfiguration approximation consists in partitioning the ensemble of configurations, i.e., in gathering them into groups. A superconfiguration represents a number of configurations of the plasma and is composed by an ensemble of supershells (i.e., groups of subshells) and their respective populations [51,52]. For example, superconfiguration represents all the configurations with shells n = 1, n = 2 and n = 3 full, a total of 19 electrons distributed in all possible ways (consistent with the Pauli exclusion principle) in the subshells 4s, 4p, 4d, 4 f , 5s, 5p and 5d (the third supershell consists of 7 orbitals 4s to 5d), and 3 electrons in the last supershell made of subshells 5 f and 5g only.
In the present work, we have one superconfiguration par ion (which means that all the existing subshells are gathered in the same supershell), and perform therefore one self-consistent calculation of electronic structure per ion.
For a superconfiguration of ionization Z * , the electron densities can be adjusted in order to ensure the constraints: and where n e (r) = n b (r) + n f (r) is the total electron density (bound + free electrons), or, for an average atom The electron density and the mono-electronic potential are calculated in a self-consistent way in the DFT/LDA theory by an iterative scheme involving Kohn-Sham orbitals and/or Fermi functions. Let k be a bound state, the associated Kohn-Sham mono-electronic orbital obeys the Schrödinger equation (throughout the paper we intentionally omit the "hat" characteristic of quantum-mechanical operators): H being the Hamiltonian of the system. The mono-electronic potential involved in the definition of H are where −Z/r is the potential of the nucleus, V e (r) the potential due to the electrons and V xc [n e (r)] the exchange-correlation potential.
In the nonrelativistic self-consistent spherical ion-cell model, we use the quantum Schrödinger states for free electrons as well as for bound electrons, i.e., the electron number density is obtained by summing over a complete set of bound and free states: where f represents the Fermi-Dirac distribution (µ is the chemical potential, Lagrange multiplier fixing the number of bound electrons of the ion): and where Y m (θ, φ) are the usual spherical harmonics. As the potential is constant outside the ion sphere, the radial wave function of a bound state outside the ion sphere can be written as y n (r) = A n r κ (−ikr) (10) where κ is the modified spherical wave function of the third kind, and the wave vector k is defined as k = √ 2 . The constant A n is determined by the normalization and the boundary conditions at R ws . The condition that the logarithmic derivative of the wave function is continuous at R ws determines the one-electron eigenvalues . The radial wave function of a free state outside the ion-sphere can be written as where j and n denote spherical Bessel functions of the first and second kinds and δ , is the phase shift. Shape resonances accompanying the smearing of a bound state in the continuum are taken into account in the density of states [11,53,54]. The bound-electron density then reads and the free-electron density One has to solve, at each iteration of the self-consistent calculation, the Poisson equation, for the determination of the electronic potential V e (r) = 1 r r 0 4πr 2 n e (r )dr + ∞ r 4πr n e (r )dr .
The total Coulomb potential of the atom (nucleus + electrons) is We assume that beyond the Wigner-Seitz sphere, the electron charge density is equal to zero, the ionic charges cancelling locally the electron density (n e (r) = 0 for r ≥ R ws ) so that, for r ≥ R ws , V e (r) = R/r and V c (r) = 0, because of the conditions (3)- (5).
It is therefore possible to limit ourselves to the calculation of the potential inside the Wigner-Seitz sphere: where Practically, we solve Equation (14) by an inward/backward algorithm (Froese method [55]): (i) From r = 0 to r = R ws : calculation of Q e (r) by a four-point quadrature formula, the electron density behaving as r α (the value of α is specified in Appendix D) when r → 0.
(ii) From r = R ws to r = 0: calculation of the second term of relation (14) by a four-point quadrature formula.
The mono-electronic energies k , solutions of Schrödinger equation, rely on the assumption that the potential V scf is equal to zero outside the sphere, which is not the case if the exchange-correlation is non zero, and must therefore be translated from The self-consistent calculation is initialized by a simple Coulomb potential, and we have seen that the one-electron eigenenergies and wave functions are obtained solving Schrödinger equation. The electron density is built using the latter, the chemical potential is determined by ensuring the integer number of electron in the superconfiguration, the new radial potential is obtained through Poisson's equation, and the next iteration starts, etc. until convergence. We use a convergence criterion on the central potential, considering that convergence is achieved at iteration n if max 0≤r≤R ws

Self-Consistent Calculations-Unscreened Case
If we impose, for a (super) configuration of ionization Z * , the cancellation of the free-electron density n f (r), one calculates an isolated ion of charge Z * , with the occupation constraint where n e (r) = n b (r) is the bound-electron density, equal to the total density. The monoelectronic potential involved in the definition of the Hamiltonian H is in that case where −Z/r is the potential of the nucleus, V b (r) the potential due to the bound electrons, and V xc [n b (r)] the exchange-correlation potential. As in the screened case, at each iteration, one has to solve the Poisson equation (see Appendix C for the calculation of the electronic potential the total Coulomb potential of the isolated ion (nucleus + electrons) still being We assume that, outside the sphere of radius R ∞ (our "numerical infinity), the bound electron density is zero (typically smaller than 10 −7 atomic units), the density of the plasma being sufficiently low so that the wave functions (or their derivatives) can be considered as In practice, R ∞ ≈ 10 R ws , but we always check that the results are not affected if we take a higher value (in the limit of the numerical capabilities of the code, of course). The mono-electronic energies k , solutions of Schrödinger Equation (6), are dependent on the assumption that the potential V scf is zero outside the sphere, which is not the case, and must be shifted from the quantity

Difference between Screened and Unscreened Energies
Let us consider an isolated (unscreened) ion of average ionization Z * . The energy of the mono-electronic Kohn-Sham orbital is k + V scf (R ∞ ) associated to the wave function ψ k of state k, and solution of Equation (6). The mono-electronic potential entering the definition of the Hamiltonian H is Let us now consider the screened ion. The energy of the mono-electronic Kohn-Sham orbital is˜ k +Ṽ scf (R ∞ ) associated to the wave functionψ k which obeysHψ k =˜ kψk . The mono-electronic potential entering the HamiltonianH is The ionization potential lowering is, for this bound state i.e., where δ k =˜ k − . This value can be obtained by subtracting two simulations, one screened and one unscreened.

First-Order Perturbation Theory
The introduction of a free-electron distribution screening the isolated ion can be understood as a perturbation δH of the Hamiltonian of the isolated ion: H = H 0 + δH. The first-order correction to energy reads and δV b (r) =Ṽ b (r) − V b (r) denotes the variation of V b due to the free-electron density. The ionization potential depression reads then, at first order 6. Discussion of the Results and Comparison to an Alternative Approach 6.1.

Analysis of the Results
As mentioned in the introduction, ten years ago, Ciricosta et al. [30] used the Linac Coherent Light Source to generate solid-density aluminum plasmas at temperatures up to 180 eV. By varying the photon energy of the X rays that create and probe the plasma, they observed the K-fluorescence and could directly measure the position of the K edge of the ions in the plasma. The values they found disagree with the predictions of the SP model, but are consistent with the (modified) EK model, which predicts significantly larger depression of the ionization potential.
The Orion laser was used by Hoarty et al. [32] to study dense plasmas created by a combination of short pulse laser heating and compression by laser driven shocks. In that way, the plasma density was systematically varied between 1 and 10 g/cm 3 using aluminum samples buried in plastic foils or diamond sheets. The aluminum plasma was heated to electron temperatures between 500 and 700 eV allowing the plasma conditions to be diagnosed by K-shell emission spectroscopy. The K-shell spectra show the effect of the ionization potential depression as a function of density. The measured data were compared to simulated spectra accounting for the change in the ionization potential by the SP and EK simple models and were found to be in closer agreement with simulations using the SP model. The IPD is obtained by a simple procedure, consisting in subtracting the isolatedatom (infinite R ws radius) subshell energy from the subshell energy subject to the plasma environment. In our modeling, the isolated-atom case is obtained using our electronicstructure code for a very large Wigner-Seitz radius (representing the limit where the electron density tends to zero). Figures 1-5 display the IPD calculated using our approach and different analytical models (DH, SP, EK, modified EK) at temperatures consistent with Ciricosta's experiment, i.e., T = 50, 70, 90, 120 and 180 eV.
The general trend is that our approach is always closer to SP for the lower charge states (keeping in mind that modified EK and SP are rather close in that case). For the highest charge states, at the lowest considered temperature (T = 50 eV), our results are closer to modified EK than SP, but the situation changes progressively as the temperature increases, and for the highest temperature (T = 180 eV) considered by Ciricosta [30], our results are closer to SP. Figures 1-5 also show the values obtained when the unscreened ionization potential, subtracted to the screened one subject to plasma effects, is taken at the plasma temperature. We can see that the quantity obtained that way (black dashed curve) increases with charge state. It is not an IPD stricto sensu, i.e., as commonly understood in the past and current works on the subject. Indeed, at plasma temperature T, excited states of a given ion are populated, and the ion "feels" the plasma environment differently as compared to the zerotemperature case (which corresponds to the "true" isolated atom). This has automatically an impact on the ionization potential. Figure 6 represents the IPD of the different charge states for an aluminum plasma with different models in conditions related to Hoarty's experiment with different modeling options. The figure shows a comparison between our work (with exchange-correlation and subtracted isolated-atom contribution taken at T = 0.001 eV), the values obtained by subtracting unscreened-atom ionization potential taken at the plasma temperature with exchange-correlation, the values obtained by subtracting unscreened-atom contribution taken at the plasma temperature without exchange-correlation, and usual analytical expressions: DH, SP, EK and modified EK.
This shows that accounting for exchange-correlation is important and that our results with T = 0.001 K isolated atom lie between EK and SP, while our results with the isolated atom taken at the real temperature are between SP and modified EK. Our results with unscreened ionization potential at T and no exchange-correlation are very close to modified EK, which is probably fortuitous (or at least difficult to interpret).

Reformulation of the Potential
The model presented in this subsection was presented by Massacrier a few years ago [46][47][48]; we will refer to it as "GM model" in the following. Let us consider the Poisson equation Massacrier suggested to use [46][47][48]: where provided that the latter integral exists. Integrating Equation (17) by parts, one gets If we apply this to the potential created by the free electrons, and with the notations used throughout the paperṼ Taking Equation (16) for the expression of the IPD and neglecting the variations of the exchange-correlation potential and the influence of the free electron on the bound-electron density, one gets i.e.,

Correction to Recover Perturbation Theory
Setting which is precisely the expression of the IPD used by Massacrier. However, perturbation theory teaches us that the IPD is equal, at first order, to Neglecting the variation of exchange-correlation and the influence of the free electrons on the bound-electron density (δB b = 0), one gets and taking Z * /R ∞ , we find that the IPD should be stricto sensu meaning that the GM approach does not take into account the term As can be checked in Tables 1-3, the contribution of such a term is small compared tõ V f (0)) [46][47][48].
The values predicted by the GM model are close to ours and always higher. They coincide exactly with our values for the lowest and highest charge states and exhibit a linear behaviour (see . Figure 12 displays the IPD of the different charge states for an aluminum plasma with different models in conditions related to Ciricosta's experiment (extreme temperatures 50 and 180 eV) compared to our the experimental values. It confirms that our results are close to the ones obtained by Massacrier (GM model), but always significantly smaller than the experimental values, whereas the latter lie between SP and EK (closer to SP for low charge sates, and to EK for high charge states). Table 1. Correction to the IPD due to Equation (18) for an Al plasma at T = 70 eV and ρ = 2.7 g/cm 3 (conditions of Ciricosta's experiment [30]

Conclusions
We have presented a multi-configuration computation of the ionization potential depression in the framework of density functional theory. The IPD is obtained by a simple procedure, consisting in subtracting the isolated-atom (infinite R ws radius) subshell energy from the subshell energy in the presence of the plasma environment. Several theoretical aspects were considered, such as the connection with perturbation theory, and the relative importance of the different contributions. In our case, the isolated atom is obtained using our electronic-structure code for a very large Wigner-Seitz radius (representing the case where the electron density tends to zero). It turns out that the temperature at which the unscreened-atom ionization potential is computed changes drastically the results; the values are different, as well as the variations.
Following the usual definition, in order to compare with analytical modified EK and SP models, the true IPD requires the evaluation of the unscreened ionization potential at the lowest temperature as possible (in the present work T = 0.001 eV, the lowest temperature which could be computed with our code), which corresponds to the isolated-atom value. The IPD, calculated in that way, decreases with charge state. We discussed also the values obtained when the unscreened ionization potential, subtracted to the screened one subject to plasma effects, is taken at the plasma temperature. The quantity obtained that way (which is not a real IPD) increases with charge state.
The general trend is that our approach is always closer to SP for the lowest charge states (keeping in mind that modified EK and SP are rather close in that case). As concerns the higher charge states, at the lowest considered temperature (T = 50 eV), our results are closer to modified EK than SP, but the situation changes progressively as the temperature increases, and for the highest considered temperature (T = 180 eV), our results are in better agreement with SP.
Concerning the Hoarty experiment, our results with T = 0.001 K isolated atom lie between EK and SP, while the values obtained with the unscreened ionization potential taken at the real temperature are between SP and modified EK. The values obtained while taking the unscreened ionization potential at the plasma temperature T and no exchangecorrelation are very close to modified EK.
As concerns the Ciricosta experiment, our results are close to the one obtained by Massacrier (GM model), but always significantly smaller than the experimental values, whereas the latter lie between SP and EK (closer to SP for low charge sates, and to EK for high charge states).
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy and ethical restrictions.

Acknowledgments:
The author is indebted to Gérard Massacrier and Gérard Dejonghe for fruitful discussions on the subject in 2015.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Differences between SP Models in the Literature
In the present work, we have used the formula of Preston et al. [50]: but Calisti et al. use an alternative formulation [38,56]: Let us transform the above relation; we have in particular and therefore (1 + z 1 ) in Ref. [38] plays the role of z(1 + z * ) in Ref. [50].

Appendix B. Value of C EK at Critical Density n c
As we have seen in Section 2, the modified Debye length reads and the radius and since and which means that C EK is proportional to 1/Z.

Appendix C. Solving the Poisson Equation
In atomic units, the Poisson equation reads It can also be put in the integral form V(r) = d 3 r n(r ) |r − r | since the Green function associated to the Laplacian verifies ∇ 2 1 r = −4πδ(r).
The potential V is given by rV(r) = Following the Froese-Fischer technique, we calculate separately z(r) by outward numerical integration (from r = 0 to R ws ), apply the boundary condition, and compute y(r) by inward integration. The following quadrature formulas are used: d r k z dr = r k f as well as d r −(k+1)y dr = − (2k + 1) r k+2 z.

Appendix D. Numerical Scheme Used for the Determination of the Wave Functions
Let us write the Schrödinger equation in the form: with F (r) = −2[ − V(r)] and lim r→0 rV(r) = Z. In order to obtain the function y(r), we resort to Noumerov's scheme [57,58]. In this "shooting" method, y left is first integrated from left (nucleus) to right (outward integration), and in a second step, y right is integrated from right (boundary) to left (inward integration). Both functions y left and y right are set equal at the matching point r t . For y left , one has with y i = y(r i ), F i = −2[ + V(r i )] and h = r i+1 − r i . If the values of y(r) at r = r a and r = r b = r a + h are known, the above numerical scheme enables one to deduce the value of y left (r) at any point of the grid. Using Duhamel's principle, the solution of Equation (A1) reads y(r) = C r +1 + r +1 2 + 1 r 0 1 − r r 2 +1 f (r )(r ) − y(r )dr .
Since r f (r)w(r) has a limit when r tends to zero, in order to calculate the integral for small values of r one can replace r f (r)w(r) by a + br, where a = −2Z. After integration, one has w(r) = 1 + r − Z + 1 + br 2(2 + 3) .
In order to compute y right , we use, still following Ref. [58], the semi-classical approximation at the boundary and integrate inward according to Then, at the semi-classical turning point: y left (r t ) = y right (r t ) provides the eigenvalue, ensuring the normalization of the wave function over the whole space. One has to check that y n (r) has exactly n − − 1 nodes. The above mentioned scheme is rather fast and robust, but another alternative would be to use the (faster) phasefunction method [58,59]. One can see from Equation (A3) that the exponent α characterizing the dependence of the electron density with respect to r close to the origin is equal to 2 , where is the orbital (or azimuthal) quantum number. As concerns the free states, only an outward integration is performed, and the solution is matched to the expression of the wave function outside the sphere given by Equation (11).