Approximate Solution of the Thomas–Fermi Equation for Free Positive Ions

: The approximate solution of the nonlinear Thomas–Fermi (TF) equation for ions is found by the Fermi method. The solution is based on the new asymptotic representation of the TF ion size valid for any ionization degree. The two universal functions and their derivatives, introduced by Fermi, are calculated by recent effective algorithms for the Emden–Fowler type equations with the accuracy sufﬁcient for majority of applications. The comparison of our results with those obtained previously shows high accuracy and validity for arbitrary values of ionization degree. This study could potentially be of interest for the statistical TF method applications in physics and chemistry.

In the TF statistical model [8], the electrostatic potential V within the ion with N bound electrons and nuclear charge Z is defined by the expression where V 0 = (Z − N)e/r 0 is the potential at the TF ion boundary r 0 . The φ(x) is the TF function that satisfies the TF equation, where x = r/r TF is measured in TF unit r TF = 1 4 ( 9π 2 2Z ) 1/3 a 0 , a 0 is the Bohr radius. This nonlinear ordinary differential equation (ODE), which is the foundation of TF method, has the following dimensionless form: together with boundary conditions φ(0) = 1 and φ(x 0 ) = 0. Here, x 0 = r 0 /r TF -is the dimensionless ion size (boundary radius) in the TF model, which has to satisfy the condition The ionization degree q is determined by the number of bound electrons N and nuclear charge Z: q = (Z − N)/Z. It could be formally shown that if N = Z, and hence q = 0 (i.e., for neutral atom), the dimensionless boundary radius x 0 stems to infinity. In this case, the Equation (1) will determine the TF function for the neutral atom. Conventionally, this solution is designated by the function φ 0 , and it is evidently universal for all Z.
For the free positive ions, the TF Equation (1) solution φ(x) depends only on the ionization degree q. According to Fermi [4], the approximate solution of (1) could be represented as a sum of two functions where δ(x) is assumed small in comparison with the function φ 0 (x). It is evident, that δ(0) = 0. Expanding the right hand side part of (1) in series up to the first order over δ/φ 0 , with account of (3), one can get the linear ODE for δ: Next, it is assumed [4] that where k(q) -is the coefficient depending on the ionization degree q, while η 0 (x) -is the universal function independent from q. Then it follows from (4) with the initial conditions η 0 (0) = 0 and η 0 (0) = 1. The solution of Equation (6), satisfying the fitting boundary conditions, could be put in the form [4] The values of η 0 (x) and η 0 (x) could be determined as from (7), as by the direct numerical integration of (6) [4]. The latter way was realized in [5,6,[8][9][10][11][12], and the results were tabulated. The interpolation formulas for the coefficient k(q) < 0 were proposed as well. So, in [6] k = −0.083q 3 , and in [12] k = −0.0542q 2.86 . However, such interpolation expressions as shown in the present work give satisfactory results only in the range of values q ≤ 0.3.
A relatively simple analytical approximate solution of Equation (4), applicable for the determination of φ(x) in the interval from 0 to x 0 , was obtained by Sommerfeld in [7]: here z = x 3 √ 144 λ 2 , λ 1 ≈ 7.772 and λ 2 ≈ 0.772. The equation for determination of z 0 and the boundary ion radius x 0 follows from (2) and (8): The expression in front of the square brackets in (8) is the Sommerfeld approximate solution of TFE for the neutral atom φ 0 [7]. The solution (8) of Equation (1) could be significantly improved if the exact numerical values are used for the function φ 0 [8]. Then, if to perform the corresponding replacement, Equation (9) for the determination of z 0 (x 0 ) is changed However, as was established in the present work, the x 0 (q) dependence determination from Equations (9) and (10) is possible only for values of q less than certain q max . The representation (3) for φ is formally invalidated if the correction δ becomes comparable with the value of φ 0 , i.e., for values of x close to x 0 . In order to avoid this and improve the approximate solution of Equation (1) according to the Fermi method, Gombas [8] proposed to define the coefficient k(q) from the equality φ( using the exact values of boundary radius x 0 . However, the exact values of x 0 are not known beforehand and should be determined during the procedure of the solution of Equation (1) taking into account (2). This complexity could be overcome using the asymptotic representation for x 0 (q). The aim of the present paper is an improvement of the Fermi method for solution of Equation (1) for the positive ions in order to get results applicable for practically any values of q with adequate accuracy for applications. That is why the functions φ 0 (x), η 0 (x) and their derivatives were calculated with high accuracy by the recently developed algorithms [27,28] and tabulated (see Appendix A). The coefficients k(x 0 ) are determined from (11) based on tabulated values of φ 0 (x) and η 0 (x), presented in the Appendix A (see Section 2). The asymptotic representation is obtained for the dependence x 0 (q) and hence for k(q), applicable for the wide range of the q variation and convenient in practical use (see Section 3). The comparison of k(q) and x 0 (q) dependencies with those obtained from the direct numerical solution of TFE and other known data is performed (see Section 3). The accuracy estimations of the derived approximate solution are also presented in the Section 3. The obtained results are summarized in Section 4.

Calculations of k(q) Coefficients
The coefficient k (see expression (5)) is the necessary element of the approximate solution of Equation (1) with the help of the representation (3). As we pointed out in the introduction, this coefficient could be determined from the equality k = −φ 0 (x 0 )/η 0 (x 0 ). In addition, x 0 and k depend on the ionization degree q. The exact value of boundary radius x 0 for an arbitrary ion is not known beforehand, as it determined from the additional condition (2) during the solution of Equation (1). There are two ways to avoid this difficulty. The first version is to find the numerical solution of Equation (1) together with x 0 for some values of the ionization degree, and after k, and next with the help of interpolation to find those quantities for any q. The second one is to use for x 0 (q) the asymptotic representation. Both versions are considered below.
Following [10], Equation (1) could be represented via the equivalent system of two ODEs of the first order: In (12) ; the initial conditions are: p(0) = 1, . The numerical solutions of the system (12) for different φ (0) were found by the Runge-Kutta method of the eighth order [29]. The values of x 0 (condition (2)) and k (the equality (11) and data from Tables A1 and A2) were determined correspondingly and presented in the Table 1 (compare with [12,16]). The dependence of k coefficient on the ionization degree q is shown in Figure 1. The solid line corresponds to k values shown in Table 1. The dotted curve corresponds to results obtained in [12] for the k values in the course of numerical solution of Equation (1). The dashed curve and the dash-dotted one are related to the two interpolation formulas, k = −0.083q 3 [6] and k = −0.0542q 2.86 [12], respectively. It is seen clearly in the figure that these interpolation formulas suit well only for q less than 0.3. The k values obtained in this work and shown in the Table 1 correspond well to data from [12]. However, it is worth noting that the data of Table 1 encircle the ionization degree values 0.01 ≤ q ≤ 0.99, in contrast to [12], where the k values are calculated only for q ≤ 0.65. Table 1. The boundary ion radius x 0 and k coefficient (Equation (11)), corresponding to the different values of the ionization degree q, are calculated using the numerical solution of the ODE system (12). q x 0 log 10 (|k|) q x 0 log 10 (|k|) q x 0 log 10 (|k|) Now consider the determination of the x 0 (q) asymptotic representation, based on the consideration of the TFE (1) analytical solution for ions in the two limiting cases of low and high ionization degree q.
The construction of a solution of Equation (1) for ions with the low ionization degree (q 1) is based on the expansion in powers of small parameter [17,19]. Skipping the details of calculations it is possible to show that, in the present case, the analytical dependence of x 0 (q) could be put in the form (see [17,19]): where Λ 0 ≈ 32.7294, λ 2 ≈ 0.772. The first two coefficients of expansion d k in (13) are d 1 = −0.91670 and d 2 = 0.02608. The expansion of solution of Equation (1) in the functional series for ions with the high ionization degree (q → 1) [14,15] over powers of ratio N/Z = 1 − q 1 is also known. This approach allows us to obtain the expansion for the boundary radius x 0 in the form: where x 0m are the certain coefficients of the expansion (14).  Table 1. The dotted curve corresponds to obtained in [12] k values in the course of numerical solution of Equation (1). The dashed and dash-dotted curves correspond to the two interpolation formulas, |k| = 0.083q 3 and |k| = 0.0542q 2.86 obtained in [6,12], respectively.
The dependence of the boundary radius x 0 on the ionization degree was also studied in the Thomas-Fermi-Dirac model in [20]. By numerical integration of corresponding equations the set of dependencies x 0 (q) were obtained for the different values of Z − N. It is shown there that for the increasing ionization degree (Z − N) → ∞ (Z N) the set of obtained curves converge to the curve x 0 (q), obtained in the TF model without the account of exchange. In so doing, the function x 0 (q) has the asymptotics Joining expressions (13) and (15) it is possible to derive the simple asymptotic representation for x 0 (q) suitable for practical calculations The comparison with the data of Table 1 shows that (16) allows us to calculate the x 0 (q) values with the root-mean-square error about 0.6% and maximal relative error not exceeding 1.2%. The dependence of x 0 value on the ionization degree q is presented in the Figure 2. The x 0 values from the Table 1 form the solid line. The dotted curve corresponds to x 0 values, calculated along with the asymptotic representation (16). The dashed and dashdotted curves are related to solutions of Equations (9) and (10), correspondingly. For solution of Equation (10) the data of Table A1 for the function φ 0 (x) is used. The analysis of solutions of Equation (9) (the Sommerfeld's solution) and (10) (the solution of Sommerfeld-Fermi) shows that the determination of self-consistent x 0 (q) dependence is possible only until the certain limiting value q max . For Equation (9) q max is around 0.65 and for Equation (10) q max is about 0.72. However, the Sommerfeld solution corresponds well to the values x 0 (q), obtained by the numerical solutions of Equation (1) only for the low ionization degrees q ≤ 0.15. The modified solution of Sommerfeld-Fermi fits to the numerical x 0 (q) values (Table 1) Table 1. The dotted line is related to the x 0 values calculated using the asymptotic representation (16). The dashed and dash-dotted curves correspond to the solution of Equations (9) and (10), accordingly. In the course of the solution of Equation (10), the data of Table A1 for φ 0 were used.
Thus, the results obtained in this section show that application of the asymptotic representation for x 0 (q) (16) or the usage of the numerically calculated values of the k(q) coefficients together with the values of functions φ 0 and η 0 from Tables A1 and A2 allows us to find approximate solution of Equation (1) for the arbitrary ionization degree q.

Approximations of Functions φ 0 (x), η 0 (x) and Their Derivatives
In cases when the lower accuracy is quite sufficient, it is possible to use the good approximate formulas for functions φ 0 (x) and η 0 (x) instead of the interpolation of known tabulated data (Tables A1 and A2). For example, the representation of φ 0 (x) in terms of the rational functions from [13] is well known The comparison with the data of Table A1 shows that such approximate formula allows to calculate the values of function φ 0 (x) with the root-mean-square error about 0.2% and the maximal relative error not exceeding 0.5%.
For the function η 0 (x), we obtained the new approximate formula with the help of the least-square method The estimate of accuracy of the Formula (18) for the values of functions η 0 (x) and η 0 (x) along with the data of Table A2 gives the root-mean-square error about 4% and maximal relative error not exceeding 8%. For example, the graphics of functions η 0 (x) and η 0 (x), based on the data of Table A2 are shown in the Figure 3 in solid lines, and using the approximate Formula (18)  Thus, it is possible to construct the approximate solution of Equation (1) for any value of the ionization degree q (with the accuracy sufficient for the majority of applications). In order to do this, by setting the q value, it is necessary to determine the value of the boundary radius x 0 with the help of the representation (16). After that, the value of coefficient k is determined from the relation (11) and the approximate Formulas (17) and (18). Finally, the solution of Equation (1) for given q is formed from the expressions (5), (3) with the account of (17) and (18).
The graphics of functions φ(x) and φ (x) using the numerical solution of the ODE system (12) (solid lines) and using the approximate solutions (17) and (18) Table A1. It could be concluded that, in total, this approximate method of the solution construction provides quite good results. The noticeable deviation of the approximate solution from the exact numerical one is seen only for the function φ (x) near the ion boundary only for q > 0.8. Also, the study demonstrates the consistency of the boundary conditions (2) and (11) Table A1 data.
The new asymptotic representation of the boundary ion radius x 0 (q) is obtained. ' The values of k coefficient are determined as from the numerical solution of TFE for '! ion as using the asymptotic representation for x 0 (q). The values x 0 (q) and k(q) are '" presented in Table 1. Note that data of the Table 1  Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest. Here we used results [27,28], based on the parametric representation of TFE so-# lution, introduced and applied in the even unpublished paper of E. Majorana, recon-$ structed along with his notes in [21]. The main peculiarities of this approach are de-% scribed in [21]. It was recently rigorously mathematically developed and generalized & on the Emden-Fowler type of equations in [27,28], where the detailed step by step de-' scription of the solution construction is given.  Table A1 data.

Conclusions
The approximate φ solution for the free positive ions of the TF equation was found by the E. Fermi method [4] in the assumption of the φ representation in the form: The values of functions φ 0 and η 0 are calculated by the new effective analyticalnumerical algorithm [27,28] with high accuracy. The obtained values of φ 0 , η 0 and their derivatives are given in Tables A1 and A2.
The new asymptotic representation of the boundary ion radius x 0 (q) is obtained. The values of k coefficient are determined from the numerical solution of TFE for ion using the asymptotic representation for x 0 (q). The values x 0 (q) and k(q) are presented in Table 1. Note that data of the Table 1 most fully cover the range of the ionization degree values 0.01 ≤ q ≤ 0.99, in distinction from other known papers.
The results of the direct k(q) numerical calculations of this work show that the validity of known earlier interpolation formulas for k are limited by values: q < 0.3.
Moreover, it follows from the analysis of the other approximate solutions of Equation (1) performed here (Sommerfeld solution (9) and Sommerfeld-Fermi solution (10)) that the determination of the x 0 (q) dependencies in these solutions is possible only up to the certain limit value q max . For Equation (9) q max ≈ 0.65, and for Equation (10) q max ≈ 0.72.
The simple analytical approximation of η 0 (x) is obtained, which along with the known analytical approximation of φ 0 (x) and the obtained here asymptotics of x 0 (q) dependence opens possibility for the approximate solution of the Equation (1) for any value of ionization degree q with an accuracy sufficient for major applications. The estimations of the accuracy of the obtained approximate solution confirm that in total this approach ensures very good results.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Calculation of Functions φ 0 (x), η 0 (x) and Their Derivatives
Here we used results [27,28], based on the parametric representation of TFE solution, introduced and applied in the even unpublished paper of E. Majorana, reconstructed along with his notes in [21]. The main peculiarities of this approach are described in [21]. It was recently rigorously mathematically developed and generalized on the Emden-Fowler type of equations in [27,28], where the detailed step by step description of the solution construction is given.
In [27,28] the algorithm of computing the TFE solution and its derivative of the many-electron neutral atom and for the given in advance any accuracy value in arbitrary point of ray is described in detail. We realized this algorithm for computing φ 0 (x) and φ 0 (x) for the interval x ∈ [0, 50], but don't repeat its rather lengthy description, referring to original papers [27,28]. We present in Table A1 only the exact values of functions φ 0 (x) and φ 0 (x) rounded to the tenth decimal sign (compare with [10,24]). It is worth to indicate only that in this algorithm the corresponding subsidiary functions of the parametric representation [27,28] are expanded in the Taylor series that in our case contained 120 terms. Table A1. Values of functions φ 0 (x) and φ 0 (x) calculated on the basis of the method from [27,28] and rounded to the tenth decimal sign. In our work the functions η 0 (x), η 0 (x) are obtained by the direct numerical integration of the system of two ODE of the first order, equivalent to (6): Here u(t) = η 0 (x), v(t) = 2η 0 (x) and t = √ x. Note that Fermi [4] used namely this Equation (6) for the determination of the functions η 0 (x), η 0 (x). The system (A1) should be supplemented by the initial conditions u(0) = 0 and v(0) = 2. The integration of system (A1) was performed by the explicit Runge-Kutta method of the eighth order [29]. The calculated in this way functions η 0 (x) and η 0 (x) are presented in Table A2 (compare with [11]). Table A2. The values of functions η 0 (x) and η 0 (x), obtained in the course of solution of the system of Equations (A1). For the function φ 0 (x) the data from Table A1 is used. The minimal value of the ionization degree q of the ion with the nuclear charge Z is determined by the value q = 1/Z. The consideration of ions with Z ≈ 100 is of interest from the practical point of view that corresponds to q = 0.01 and the value of boundary ion radius x 0 ≈ 34 (see Table 1). Thus the represented in Tables A1 and A2 values of functions φ 0 (x), φ 0 (x), η 0 (x) and η 0 (x) for x ∈ [0, 50], cover in full all x 0 , corresponding to q ≥ 0.01.