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}{dx^{2}}}={\frac {1}{\sqrt {x}}}y^{3/2}$

subject to the boundary conditions

$y(0)=1,\quad \quad {\begin{cases}y(\infty )=0\quad {\text{for neutral atoms}}\\y(x_{0})=0\quad {\text{for positive ions}}\\y(x_{1})-x_{1}y'(x_{1})=0\quad {\text{for compressed neutral atoms}}\end{cases}}$

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\rightarrow cx,\ y\rightarrow 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}{dx^{2}}}-6x{\frac {du}{dx}}+12u=u^{3/2}$

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

${\frac {d^{2}u}{dt^{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\rightarrow 0$ as $x\rightarrow \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}{dt^{2}}}=w^{3/2},\quad w(0)=0,\ w(\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\rightarrow \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{aligned}y(x)={}&1-Bx+{\frac {1}{3}}x^{3}-{\frac {2B}{15}}x^{4}+\cdots {}\\[6pt]&\cdots +x^{3/2}\left[{\frac {4}{3}}-{\frac {2B}{5}}x+{\frac {3B^{2}}{70}}x^{2}+\left({\frac {2}{27}}+{\frac {B^{3}}{252}}\right)x^{3}+\cdots \right]\end{aligned}}$

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$.