Thomas–Fermi equation

Numerical solutions of the Thomas–Fermi equation

In mathematics, the Thomas–Fermi equation for the neutral atom is a second order non-linear ordinary differential equation, named after Llewellyn Thomas and Enrico Fermi, which can be derived by applying the Thomas–Fermi model to atoms. The equation reads

$\frac{d^{2}y}{d{x}^2}=\frac{1}{\sqrt{x}}{y}^{3/2}$

subject to the boundary conditions

$y(0)=1,\{\begin{array}{c}y(\mathrm{\infty })=0\text{for neutral atoms}\\ y(x_0)=0\text{for positive ions}\\ y(x_1)-x_1{y}^{\prime }(x_1)=0\text{for compressed neutral atoms}\end{array}$

If $y$ approaches zero as $x$ becomes large, this equation models the charge distribution of a neutral atom as a function of radius $x$. Solutions where $y$ becomes zero at finite $x$ model positive ions. For solutions where $y$ becomes large and positive as $x$ becomes large, it can be interpreted as a model of a compressed atom, where the charge is squeezed into a smaller space. In this case the atom ends at the value of $x$ for which $dy/dx=y/x$.

Transformations

Introducing the transformation $z=y/x$ converts the equation to

$\frac{1}{x^2}\frac{d}{dx}\left(x^2\frac{dz}{dx}\right)-z^{3/2}=0$

This equation is similar to Lane–Emden equation with polytropic index $3/2$ except the sign difference. The original equation is invariant under the transformation $x\to cx,y\to c^{-3}y$. Hence, the equation can be made equidimensional by introducing $y=x^{-3}u$ into the equation, leading to

$x^2\frac{d^{2}u}{d{x}^2}-6x\frac{du}{dx}+12u=u^{3/2}$

so that the substitution $x=e^t$ reduces the equation to

$\frac{d^{2}u}{d{t}^2}-7\frac{du}{dt}+12u=u^{3/2}.$

Treating $w=du/dt$ as the dependent variable and $u$ as the independent variable, we can reduce the above equation to

$w\frac{dw}{du}-7w=u^{3/2}-12u.$

But this first order equation has no known explicit solution, hence, the approach turns to either numerical or approximate methods.

Sommerfeld's approximation

The equation has a particular solution $y_p(x)$, which satisfies the boundary condition that $y\to 0$ as $x\to \mathrm{\infty }$, but not the boundary 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 which can satisfy the other boundary condition in 1932. If the transformation $x=1/t,w=yt$ is introduced, the equation becomes

$t^4\frac{d^{2}w}{d{t}^2}=w^{3/2},w(0)=0,w(\mathrm{\infty })\sim t.$

The particular solution in the transformed variable is then $w_p(t)=144t^4$. So one assumes a solution of the form $w=w_p(1+\alpha t^{\lambda })$ and if this is substituted in the above equation and the coefficients of $\alpha$ are equated, one obtains the value for $\lambda$, which is given by the roots of the equation ${\lambda }^2+7\lambda -6=0$. The two roots are ${\lambda }_1=0.772,{\lambda }_2=-7.772$, where we need to take the positive root to avoid the singularity at the origin. This solution already satisfies the first boundary condition ($w(0)=0$), so, to satisfy the second boundary condition, one writes to the same level of accuracy for an arbitrary $n$

$W=w_p(1+\beta t^{\lambda }{)}^n=[144t^3(1+\beta t^{\lambda }{)}^n]t.$

The second boundary condition will be satisfied if $144t^3(1+\beta t^{\lambda }{)}^n=144t^3{\beta }^n{t}^{\lambda n}(1+{\beta }^{-1}{t}^{-\lambda }{)}^n\sim 1$ as $t\to \mathrm{\infty }$. This condition is satisfied if $\lambda n+3=0,144{\beta }^n=1$ and since ${\lambda }_1{\lambda }_2=-6$, Sommerfeld found the approximation as $\lambda ={\lambda }_1,n=-3/{\lambda }_1={\lambda }_2/2$. Therefore, the approximate solution is

$y(x)=y_p(x)\{1+[y_p(x){]}^{{\lambda }_1/3}{\}}^{{\lambda }_2/2}.$

This solution predicts the correct solution accurately for large $x$, but still fails near the origin.

Solution near origin

Enrico Fermi provided the solution for $x\ll 1$ and later extended by Edward B. Baker. Hence for $x\ll 1$,

$\begin{array}{cc}y(x)=& 1-Bx+\frac{1}{3}{x}^3-\frac{2B}{15}{x}^4+\cdots \\ & \cdots +x^{3/2}[\frac{4}{3}-\frac{2B}{5}x+\frac{3B^2}{70}{x}^2+(\frac{2}{27}+\frac{B^3}{252})x^3+\cdots ]\end{array}$

where $B\approx 1.588071$.

It has been reported by Salvatore Esposito that the Italian physicist Ettore Majorana found in 1928 a semi-analytical series solution to the Thomas–Fermi equation for the neutral atom, which however remained unpublished until 2001. Using this approach it is possible to compute the constant B mentioned above to practically arbitrarily high accuracy; for example, its value to 100 digits is $B=1.588071022611375312718684509423950109452746621674825616765677418166551961154309262332033970138428665$.