Generalized Thomas-Fermi equation: existence, uniqueness, and analytic approximation solutions

1.   Introduction

The classical Thomas-Fermi problem for the neutral atom is a second-order non-linear ordinary differential equation, named after Llewellyn Thomas and Enrico Fermi ^[1,2,3,4,5] which can be derived by applying the Thomas-Fermi model to atoms. The Thomas-Fermi model assumes that all electrons are subject to the same conditions and energy conservation law, and has potential energy $e\Phi$^[6] so when assuming that the potential is spherically symmetric, then the charge density $\rho$ and the potential energy are related through the Poisson's equation

where $\hbar$ is the Planck's constant, $r$ is the distance from the nucleus, and $\rho$ is given by

where $e$ is the electronic charge and $m$ is the mass. Substituting in the above equation yields

with the corresponding boundary conditions:

where $Z$ is the atomic number, and

By introducing the following transformation $r = \mu x$ for some appropriate parameter $\mu$ and $y = \dfrac{r\Phi(r)}{eZ}$, we arrive at the so-called differential equation of the Thomas-Fermi equation

This equation models the charge distribution of a neutral atom as a function of the radius $x$. It should be noted here that the basic Thomas-Fermi (TF) model for ions is subject to the boundary conditions $y(0) = 1$ and $y(x_{0}) = 0$, where $x_{0} > 0$ is the dimensionless ion size which measures the boundary radius and satisfies the relation $-x_{0}y^{\prime}(x_{0}) = q$ for some ionization factor $q$. When the nuclear charge equals the number of (bound) electrons, then $q = 0$ which occurs when $x_{0}\longrightarrow\infty$, and the the problem, in this case, describes the neutral atom model case ^[7]. The TF equation (1.6) has connections to other important partial differential equations, for example, it is considered a special case of the well-known Poisson equation, and can also be viewed as an Euler-Lagrange equation associated with the Fermi energy ^[8]. The Thomas-Fermi model has deep connections to the quantum gravity theory where it is reformulated at the Planck scale ^[9,10,11].

The equation has a particular solution $y_{p}(x)$, which satisfies the boundary condition $y\rightarrow 0$ as $x\rightarrow \infty$, but not the initial condition $y(0) = 1$. This particular solution is $y_{p}(x) = {\frac {144}{x^{3}}}$.

Arnold Sommerfeld used this particular solution and provided an approximate solution that can satisfy the other boundary condition ^[4]:

where $\lambda_1 = 0.772$ and $\lambda_2 = -7.772$. This solution predicts the correct solution accurately for large $x$ but still fails near the origin. A considerable amount of literature is devoted to the numerical solutions for the classical version of this problem ^{[12,13,14,15,16]}. In ^[12], numerical solution was obtained using the variational principle. J. Boyd ^[13] obtained a numerical solution using rational Chebyshev functions. The authors in ^[14,15] obtained numerical solutions using a spectral method based on the fractional order of rational Bessel functions. Pikulin ^[16] used a semi-analytical numerical method to compute the solution. Furthermore, different methods such as homotopy analysis and iterative methods are used to investigate the approximate solutions to this problem, see, e.g., ^{[17,18,19,20,21]}.

However, a free boundary value issue is also implemented to approach the initial Thomas-Fermi equation. As a result, the free boundary value issue is changed into a nonlinear boundary value problem that is defined on a closed interval. An adaptive approach is used to tackle the issue utilizing the moving mesh finite element method ^[22].

The present paper investigates the generalized boundary value problem of the Thomas-Fermi equation

where

and we assume that $0\leq y(x)\leq 1$.

The Thomas-Fermi equation is a special case of this equation when $p = 1$. As pointed out in ^[23], this generalized TF equation is related to non-integrable Abel equations, and therefore no closed solutions are possible for any case of this type of equation.

In this paper, we aim to provide an analytic approximate solution in explicit form for problem (1.9) with (1.10). In Section 2, we establish a theorem that provides the lower and upper bounds of the solution $y$ and guarantees the existence of the solution to this problem as well as a theorem on the uniqueness of the solution. In Section 3, we present analytic approximate solutions in different explicit forms to this problem. Also, an interesting variation of the Adomian decomposition method (ADM) ^{[24,25,26,27,28,29,30,31,32,33]} is presented, which allows the determination of the solution in an easily-computed series. In Section 4, we carry out an analysis of the solution and compare it with other numerical solutions.

2.   An existence and uniqueness theorem

2.1.   Existence theorem

We first prove a result on the double inequalities for the lower and upper bounds of the solution $y$, which is an important tool in the proof of the existence of the solution to problem (1.9) with (1.10).

Theorem 2.1. The generalized boundary value problem of TF equation (1.9) with (1.10) has at least one solution $y \in \mathbb{C}^{2}[0, \infty)$ such that

where $y_{1}(x) = \frac{1}{\left(x+1\right)^{\frac{1+\sqrt{5}}{2}}}$, $y_{2}(x) = -x\Big(\big(K(0, x)-K(1, x)\big)\Big)$, and $K$ is the modified Bessel function of the second kind.

Proof. Using the following inequality:

we obtain

which will be helpful later.

In view of $0\leq y(x)\leq 1$, it follows that

and with $y^{\frac{p}{p+1}}\leq 1$, we have

Consequently,

Hence,

On the other hand, in view of the solution $y$ remains in the interval $[0, 1]$ and since $\frac{p}{p+1} < 1$, we have

Using now $x < x+1$ to obtain $\frac{1}{x} > \frac{1}{x+1}$. Hence, $\left(\frac{ 1}{ x}\right)^{\frac{p}{p+1}} > \frac{ 1}{\left(x+1\right)^{\frac{p}{p+1}}}$. Thus,

or

It can be checked easily that

Hence,

Thus, from (2.7) and (2.12), we obtain

where

For comparison purposes, we have the following linear boundary value problems:

and

Then, suitable comparison problems are

and

To find the solution $y_{1}$ of problem (2.17), we write

Let $\xi = x+1$. Thus this equation becomes

The substitution $\xi = e^{-t}$ leads to a constant coefficient linear equation

Thus,

where $C_{1}$ and $C_{2}$ are two constants. Using the boundary conditions $y_{1}(0) = 1, \ y_{1}(\infty) = 0$ to find $C_{1} = 1$ and $C_{2} = 0$. This gives

To find the solution $y_{2}$ of problem (2.18), we bring back the form of the confluent hypergeometric equation with parameters $a$ and $b$ ^[34,35,36]:

which has a regular singularity at $0$ and an irregular one at infinity; and whose solution is $_1F_1(a; c; x)$.

The equation of problem (2.18) can be simply written as

We introduce the change of variables $\xi = 2x$ and $y_{2}(x) = v(\xi)$. Then

The transformation $v(\xi) = \xi e^{-\frac{\xi}{2}}w(\xi)$ leads to

which is the confluent hypergeometric equation with parameters $a = 3/2$ and $b = 2$.

Thus the general solution of problem (2.18) is given in terms of the modified Bessel functions as

where $I$ and $K$ are the modified Bessel functions of the first and second kind, respectively. $c_1$ and $c_2$ are arbitrary constants, which can be determined from the boundary conditions. Indeed, to satisfy these conditions $y_2(0) = 1$ and $y_2(\infty) = 0$, we get $c_1 = 0$ and $c_2 = -1$, and so the required solution is given by

For small $x$, we have

Hence, the condition $y_{2}(0) = 1$ is satisfied.

We are now able to apply the method of upper and lower solutions. For more details about this technique, we refer the reader to (Chapter 7, ^[37]), which is applicable when $f(x, y)$ has a singularity at $x = 0$ and the Lipschitz constants $L_{1}(x) = -\frac{ 1}{ \left(x+1\right) ^{2} }$ and $L_{2}(x) = -\left(1+\frac{1}{x}\right)$ are functions of the independent variable $x$ and continuous everywhere except for $L_{2}(x)$ at $x = 0$.

It should be noted here that $y_{1}$ and $y_{2}$ are both twice continuously differentiable and satisfy the above differential inequalities functions (2.15) and (2.16) on $(0, \infty)$ with $y_{1} < y_{2}$. Furthermore, the function $f(x, y)$ is continuous and bounded in

This completes the proof.

To show the variation of these two extremum functions, we present in Figure 1 the variation of the upper and lower functions in terms of the independent variable $x$. In addition, the extremum functions are independent of the parameters $p$, which makes them the optimum functions for all kinds of Thomas-Fermi equations.

2.2.   Uniqueness theorem

Theorem 2.2. The generalized boundary value problem of TF equation (1.9) with (1.10) has at most one solution $y \in \mathbb{C}^{2}[0, \infty)$.

Proof. To obtain an important result on the uniqueness, we assume that $\bar{y}_{1}$ and $\bar{y}_{2}$ are two different solutions to problem (1.9) with (1.10). Then,

and

where $g(x, y) = y \left(\frac{y}{x}\right)^{\frac{p}{p+1}}$.

Consider the positive function $h(x) = \frac{1}{2}\left(\bar{y}_{1}-\bar{y}_{2}\right)^{2}$. Thus $h$ vanishes at zero and infinity. Therefore, if is not identically zero it must have a positive maximum at a point $\bar{x}$, where $\bar{x} > 0$. Thus, its graph is concave down at $\bar{x} > 0$, and we have

Since

or

Hence,

From (2.32) and (2.33), we have

Applying the mean value theorem to the function $g$ with respect to $\bar{y}$, we obtain

where $0\leq\bar{y}_{1} < \bar{y}^{*} < \bar{y}_{2}\leq 1$.

On the other hand, differentiating the function $g(x, y)$ with respect to $y$, we obtain

Hence,

Consequently,

Substituting this into (2.37), we obtain

which contradicts the assumption that $\frac{2p+1}{p+1} > 0$, $\frac{\bar{y}^{*}(\bar{x})}{\bar{x}} > 0, \ \bar{x} > 0$ and $\left(\bar{y}_{1}(\bar{x})-\bar{y}_{2}(\bar{x})\right)^{2} > 0$. So $h(x) = \frac{1}{2}\left(\bar{y}_{1}-\bar{y}_{2}\right)^{2}\equiv 0$. This shows the uniqueness of the solution and completes the proof of the theorem.

3.   Explicit approximate solutions

We conclude here based on Theorem 2.1, which may offer advantages in finding out lower and upper solutions of our problem (1.9) with (1.10) in explicit forms such that $y_{1} \leq y \leq y_{2}$ on $[0, \infty)$, and consequently we should expect $y$ to take similar explicit forms in the whole region with the corresponding boundary conditions.

3.1.   The first approximation

To obtain an approximate solution $y$ to problem (1.9) with (1.10), we first make the following approximation.

A possible linear approximation of a function $f(x)$ at $x = x_{0}$ may be obtained using the equation of the tangent line

If we choose $f(x) = \sqrt{\beta x}, \ \beta > 0$ and $x_{0} = \frac{1}{\beta}$, then

when $x$ is close enough to $x_{0} = \frac{1}{\beta}$. Hence,

Substituting (3.3) into the nonlinear term of the ODE of problem (1.9), we obtain

For the solution $y$ of the approximate equation (3.4), by Theorem 2.1, we expect that the solution $y$ can be obtained in the form

where $\beta > 0$ and $m < 0$ are two parameters to be determined. Inserting the ansatz given by (3.5) into Eq (3.4), we obtain

If we assume that $m-2 = \frac{mp}{p+1}+m-\frac{2p}{p+1}$, that is $m = -\frac{2}{p}$, then, we derive the following relation between the parameters

that is

Thus the first analytic approximate solution to the generalized TF equation is given by

3.2.   The second approximation

The term $\left(\beta x\right)^{\frac{p}{p+1}}$ can be approximated by $\left(1+\beta x\right)^{\frac{p}{p+1}}$ for sufficiently large values of $\beta x$; that is,

Substituting (3.10) into the nonlinear term of the ODE of problem (1.9), we obtain

It follows that $m = 1+\frac{2}{p}$ and $\beta = \left[\left(\frac{p+2}{p}\right)^{2}+\frac{p+2}{p}\right]^{-\frac{p+1}{p+2}}$.

Thus, the second analytic approximate solution to the generalized TF equation for $x$ large is given as

3.3.   The third approximation

For $x$ near $1$, we can substitute $x\approx 1$ in the denominator of the nonlinear term of the ODE of problem (1.9) to find

It follows that $m = 2+\frac{2}{p}$ and $\beta = \left[\left(\frac{2p+2}{p}\right)^{2}+\frac{2p+2}{p}\right]^{-\frac{1}{2}}$.

Thus, the third analytic approximate solution to the generalized TF equation for $x$ near $1$ is given as

Thus, our approximate solutions can be obtained by direct approaches.

3.4.   The fourth approximation using ADM

In this section, we consider an interesting variation of the modified Adomian decomposition method (ADM) ^{[24,25,26,27,28,29,30,31,32,33,34]}, which permits the determination of the solution of nonlinear initial-boundary value problem (1.9) with (1.10).

Rewrite the ODE of problem (1.9) with (1.10) in Adomian's operator-theoretic form

where $L = \frac{d^{2}}{d x^{2}} \ \mbox{and} \ N(y) = y^{\frac{2p+1}{p+1}}$.

Applying $L^{-1}$ to both sides of Eq (3.15) and using the initial condition $y(0) = 1$, we obtain

where $B = y'(0)$ is an unknown constant to be determined by using the boundary condition $y(\infty) = 0$.

According to the Adomian decomposition method ^{[24,25,26,27,28,29,30,31,32,33]}, assuming the decomposition

where $A_{n}$ are the Adomian polynomials ^[24,25,33]. Thus, Eq (3.16) becomes

We identify

Hence, a new recurrence relation for $y_{n}, \ n\geq 0$, is established as

where the Adomian polynomials $A_{n}$ ^[24,25,33] for the $N(y) = y^{\frac{2p+1}{p+1}}$ term are

The first few components of the solution $y_{n}, \ n\geq 0$ are given by

Hence,

3.4.1.   Computation of $B$

It remains now to apply the second boundary condition $y\rightarrow 0$ as $x\rightarrow \infty$ to the function $y(x)$. This boundary condition cannot be applied directly to the series (3.23). Recall that it is customary to combine the series solutions obtained by the decomposition method with the Padé approximants to provide an effective tool to treat boundary value problems on an infinite or semi-infinite interval ^[33]. To illustrate this, we choose $p = 1$. For convenience, we list below, by using (3.21), few terms of the Adomian polynomials $A_{n}$

The first few components of the solution $y_{n}, \ n\geq 0$, are given by

Hence,

Setting $x^{\frac{1}{2}} = \xi$ into (3.26), we obtain

which is indeed the same approximation of $y$ that obtained by Baker in 1930 ^[5] and Wazwaz ^[33]. In applying the boundary condition $y(\infty) = 0$ to the diagonal Padé approximants $P_{10, 10} = [10 / 10]$, we obtain the approximation for the initial slope $B = y'(0) = -1.588077$, which is a very good approximation to accuracy $10^{-5}$ comparing to the value obtained by Parand et al. as $-1.588071$ ^[14]. These values are also in good agreement with the obtained numerical value $y'_{n}(0) = -1.564036$ for $p = 1$.

4.   Analysis of solutions

We are now in the position to explore some mathematical results and investigate the numerical treatment of the boundary value problem (1.9) with (1.10). In Figure 2, we present the different solutions of problem (1.9) with (1.10) with the particular case $p = 1$. The first approximation (solid blue line) is in good agreement with the numerical solution and Sommerfeld's approximation. On the other hand, the third approximation (black dash-dotted line) is in good agreement with the numerical solution for small values of the independent variable $x$. While the second approximation diverges slightly from the other solutions for small and intermediate values of the independent variable $x$. All solutions coincide together for large values of $x$. Due to the potential limits of the numerical volume ^[38,39,40], we chose the maximum value of the independent variable as $x = 14$. The numerical solution is obtained, using the Maple software, and the available mid-rich sub-method, which is a midpoint method with the same enhancement schemes. So, the midpoint sub-methods are capable of handling harmless end-point singularities that the trapezoid sub-methods cannot. For the enhancement schemes, Richardson extrapolation is generally faster, but deferred corrections use less memory on difficult problems ^[41,42].

In addition, we present in Tables 1 and 2 a comparison between the numerical solution and different proposed approximations for the case $p = 1$, for small and large values of the independent variables $x$. These numerical values show clearly that the first and third approximations agree very well with the numerical solution in all ranges of the independent value $x$.

Now, we can explore other interesting cases with $p\ne 1$, to show the efficiency of the suggested approximations and their validity ranges. In Figure 3, we present the different solutions of (1.9) with (1.10) with the particular cases $p = 2, 3$. The first and third approximations (solid blue line, black dash-dotted line) are in good agreement with the numerical solution for small values of the independent variable $x$. On the other hand, the third approximation remains in good agreement with the numerical solutions, while the first approximation diverges from the numerical solution by increasing the parameter $p$. The second approximation is still larger than all approximations over the small and intermediate domains of $x$.

Our overall findings demonstrate that it is possible to acquire a good approximation to the generalized TF equation. The charge distribution of a neutral atom as a function of radius $x$ is also well-known to be described by this equation if and only if $y(x)$ approaches zero as $x$ grows in size. Solutions with $y(x) = 0$ at a finite $x$ are used to mimic positive ions. For solutions where $y(x)$ becomes significant and positive as $x$ increases significantly, it can be viewed as a model of a compressed atom, where the charge is squeezed into a smaller region. These broad comments are adequately supported by our plots. The proposed investigation might be useful in dense media where quantum gravity's effects could be felt strongly.

5.   Conclusions

The goal of this study is to solve the generalized TF equation which governs several physical issues, such as quantum systems, that naturally differ significantly from Fermi or Bose statistics, as well as some astrophysical or cosmological contexts where quantum electrostatics may exhibit more intertwined screening effects. The TF equation is modeled in this investigation as a singular boundary value problem with an upper and lower solution theory. The existence-construction of the aforementioned upper-lower solutions is also explored. Excellent approximations are proposed and the obtained results are in good agreement with those obtained numerically. We anticipate that the approximation solutions we have presented will be useful in assisting with the investigation of the TF model-governed physics issues.

Acknowledgements

The authors extend their appreciation to the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University for funding this work through Research Group no. RG-21-09-14.

Conflict of interest

The authors declare that they have no competing interests.