Self-Consistent Thomas-Fermi Approximation

• Self-Consistent Thomas-Fermi approximation is a semiclassical method that determines local densities and potentials together, replacing discrete orbitals with local Fermi-gas relations.
• It is applied across atomic, nuclear, and low-dimensional models to efficiently capture density-functional behavior and correctly impose global constraints like charge neutrality and thermal equilibrium.
• Numerical strategies such as Fourier methods, radial grid integrations, and imaginary-time propagation illustrate its adaptability while highlighting trade-offs in accuracy and computational complexity.

The self-consistent Thomas–Fermi approximation is a semiclassical mean-field framework in which densities, local occupation functions, or auxiliary fields are determined together with the potentials they generate, subject to global constraints such as fixed particle number, charge neutrality, or thermodynamic equilibrium. In its classical form it replaces discrete orbitals by local Fermi-gas relations; in later variants it appears in orbital-free density-functional theory, relativistic mean-field treatments of nuclei and Wigner–Seitz cells, trapped-gas functionals with gradient corrections, and magnetic or phase-coupled generalizations in lower dimensions (Zhang et al., 2014, Li et al., 31 Jul 2025, Lehtomäki et al., 2017, Levitt et al., 18 Apr 2025).

1. Variational structure and defining features

Across its major realizations, the approximation is organized around a free-energy or energy functional whose stationarity conditions couple local matter variables to self-generated fields. In relativistic non-uniform nuclear matter, nucleon densities and scalar densities are computed from local Fermi–Dirac distributions, while meson and Coulomb fields are obtained from Poisson-type equations sourced by those same densities; the Wigner–Seitz cell radius is then selected by minimizing the free energy per baryon at fixed $T$, $Y_p$, and $\rho_B$ (Zhang et al., 2014). In finite nuclei, the same logic appears in zero-temperature form: local Fermi momenta $k_i^F(r)$ are determined from constant chemical potentials and position-dependent $\sigma$, $\omega$, $\rho$, and Coulomb fields, and the fields are iterated until they reproduce the densities that generated them (Li et al., 31 Jul 2025). In a fully local orbital-free electronic variant, the Euler equation becomes algebraic in $\rho^{1/3}$, so the density is obtained pointwise from the external potential and a chemical potential fixed by normalization (Rasanen et al., 2013). In trapped 2D dipolar gases, the same closure is written as a nonlinear Schrödinger-like equation for $\psi=\sqrt{n}$, with the effective potential depending on $\psi^2$ and $Y_p$0 (Zyl et al., 2012).

Context	Primary variables	Closure mechanism
Atomic orbital-free electrons	$Y_p$1, $Y_p$2, $Y_p$3	Euler equation and normalization
Relativistic nuclei / Wigner–Seitz cells	$Y_p$4, $Y_p$5, $Y_p$6, $Y_p$7	Local Fermi distributions plus field equations
Trapped gases / gauge-coupled systems	$Y_p$8, $Y_p$9, or $\rho_B$0	TFvW or magnetic-TF minimization

The self-consistent label therefore refers less to a specific functional form than to a closure principle: the fields entering the local kinetic or thermodynamic relation are the same fields produced by the resulting matter distribution.

2. Atomic and orbital-free realizations

In the atomic setting, the canonical Thomas–Fermi boundary-value problem is the dimensionless equation

$\rho_B$1

with boundary conditions $\rho_B$2 and $\rho_B$3. A useful reformulation introduces $\rho_B$4, which converts the problem into a quasi-linear equation and makes it possible to construct rational-power approximants that match the Majorana small-$\rho_B$5 expansion and the Sommerfeld asymptotic law $\rho_B$6 within a single analytic expression (Esposito et al., 2020). This mathematical line belongs to the older statistical-atom tradition in which self-consistency means solving Poisson screening and local Fermi filling simultaneously.

Later orbital-free atom models altered the basic variables rather than abandoning the semiclassical structure. The Englert–Schwinger model is a self-consistent statistical atom model written as a potential functional $\rho_B$7, with density obtained from $\rho_B$8, exchange treated at Dirac-LDA level, and strongly bound electrons treated separately from the semiclassical continuum. In the benchmark reported there, the model generally offers the same accuracy as TFD-$\rho_B$9vW, but corrects the failure in the Pauli potential near the nucleus and is limited chiefly by its inability to describe low-$k_i^F(r)$0 atoms (Lehtomäki et al., 2017). A different local orbital-free construction replaces the nonlocal electron–electron term by

$k_i^F(r)$1

so that the resulting Euler equation is algebraic in $k_i^F(r)$2; the chemical potential is then fixed by $k_i^F(r)$3, and reasonable total energies are obtained for Hooke’s atoms, jellium, and model atoms up to $k_i^F(r)$4 electrons (Rasanen et al., 2013).

A more radical atomic variant recasts the Pauli principle through polymer self-consistent field theory. In that construction, electrons are modeled as Gaussian threads in a four-dimensional thermal-space, and the Pauli potential is postulated as a Flory–Huggins excluded-volume interaction between different electron pairs. Polymer scaling arguments show that the excluded-volume term reproduces the Thomas–Fermi $k_i^F(r)$5 scaling in the uniform limit, while the self-consistent solution for the first eighteen elements yields radial densities with correct shell structure and total binding-energy errors below $k_i^F(r)$6 for the lightest atoms and $k_i^F(r)$7 or less for atoms heavier than nitrogen (Thompson, 2022). This goes beyond ordinary TF precisely because the Pauli term is not inserted as a fixed local power of the density but generated self-consistently from the densities of different contour species.

Finite-temperature atomic TF theory also admits a systematic leading correction from strongly bound electrons. In the ion-sphere setting, the TF free energy scales as

$k_i^F(r)$8

while inclusion of quantum K-shell electrons produces

$k_i^F(r)$9

The $\sigma$0 Scott term is therefore a $\sigma$1 correction relative to TF, more dominant than the $\sigma$2 relative suppression of gradient, exchange, and correlation corrections. The construction remains valid only while the L shell stays non-ionized (Segev et al., 2014).

3. Relativistic nuclear matter and finite nuclei

In nuclear physics, the self-consistent Thomas–Fermi approximation is most closely associated with relativistic mean-field models of non-uniform matter. At subnuclear density, the medium is modeled as a body-centered cubic lattice of heavy nuclei surrounded by dripped nucleons and electrons, then reduced to a spherical Wigner–Seitz cell of radius $\sigma$3. For each $\sigma$4, $\sigma$5, and $\sigma$6, the favored state is obtained by minimizing the free energy with respect to $\sigma$7, while local proton and neutron densities follow from Fermi–Dirac occupations in spatially varying $\sigma$8, $\sigma$9, $\omega$0, and Coulomb fields (Zhang et al., 2014). Surface and Coulomb energies then emerge directly from the gradient terms in the meson and electromagnetic fields, rather than from an imposed liquid-drop parameterization.

This STF scheme was explicitly compared with the parameterized Thomas–Fermi procedure used in the Shen equation of state. In the latter, nucleon profiles are constrained to a specific functional form and the surface term is represented phenomenologically by

$\omega$1

For $\omega$2, the PTF free energy per baryon is slightly lower than STF by $\omega$3–$\omega$4 MeV around $\omega$5–$\omega$6, largely because surface energy is underestimated; increasing $\omega$7 to $\omega$8 brings the gradient and Coulomb energies much closer to STF values. The same calculations find a spherical bubble phase near $\omega$9, with free-energy differences below $\rho$0 relative to droplets but noticeable changes in chemical potentials (Zhang et al., 2014).

The same relativistic STF formalism can be applied to isolated nuclei. There the unknowns are the radial meson and Coulomb fields and the local Fermi momenta $\rho$1, $\rho$2, with proton and neutron numbers fixed by

$\rho$3

A recent comparison with full relativistic mean-field calculations, performed with the same nonlinear $\rho$4–$\rho$5–$\rho$6 interaction for nuclei from $\rho$7 to $\rho$8, shows the expected pattern: STF reproduces average bulk trends but yields smoother density profiles, smaller neutron skins than RMF for the same parameter set, and no shell or spin–orbit structure because nucleons are treated as a local Fermi gas rather than by Dirac orbitals (Li et al., 31 Jul 2025).

Hot $\rho$9 hypernuclei introduce an additional refinement: the subtraction procedure. Two solutions are computed in a large box, one for liquid plus gas ($\rho^{1/3}$0) and one for gas alone ($\rho^{1/3}$1); the physical hypernuclear density is then

$\rho^{1/3}$2

This construction makes bulk observables independent of box size and permits direct study of liquid–gas coexistence in strange finite systems. Within that framework, the $\rho^{1/3}$3 central density is very sensitive to temperature, the $\rho^{1/3}$4 radius becomes very large at high temperature, and the single-$\rho^{1/3}$5 binding energy decreases strongly with increasing temperature (Hu et al., 2016).

4. Low-dimensional trapped gases and semiclassical pairing

In two dimensions, the form of the TF kinetic functional changes because the density of states is constant. For a spin-polarized 2D Fermi gas,

$\rho^{1/3}$6

and for a harmonically trapped dipolar gas the full Thomas–Fermi–von Weizsäcker functional becomes

$\rho^{1/3}$7

Writing $\rho^{1/3}$8 converts the Euler equation into

$\rho^{1/3}$9

which is then solved by imaginary-time propagation. The von Weizsäcker term is not generated by the standard 2D gradient expansion in the usual way, but it restores the edge correction missing from pure TF and smooths the unphysical sharp cutoff of the LDA density (Zyl et al., 2012).

A different 2D realization occurs in “flat” artificial atoms. There the nontrivial point is the combination of a 2D density of states with 3D electrostatics: electrons move in a plane, but the Hartree field is generated by the ordinary 3D Coulomb law. The resulting integral TF equation is linear but involves a singular Coulomb kernel. Diagonalization of that kernel yields a fully self-consistent analytic solution for general radial confinement, and for harmonic confinement the entire problem collapses to the two universal scaling variables

$\psi=\sqrt{n}$0

In the small-dot limit $\psi=\sqrt{n}$1, $\psi=\sqrt{n}$2; in the large-dot limit $\psi=\sqrt{n}$3, $\psi=\sqrt{n}$4 and the density approaches the classical profile proportional to $\psi=\sqrt{n}$5 (Ovchinnikov et al., 2013).

Self-consistent TF ideas also extend to superfluid pairing. In that case the $\psi=\sqrt{n}$6 limit is taken not of the static density functional alone but of the weak-coupling gap equation written in the basis of the mean field. The result is an energy-space integral equation for $\psi=\sqrt{n}$7 involving the TF level density $\psi=\sqrt{n}$8 and semiclassical pairing matrix elements $\psi=\sqrt{n}$9. Along isotopic chains, the TF gaps smooth out the quantal arch structure, but that structure can be almost recovered if shell fluctuations are reintroduced into the level density. The same formalism predicts strong suppression of pairing at the drip line and was used to analyze pairing in Wigner–Seitz cells of the neutron-star inner crust (Viñas et al., 2011).

5. Generalizations to gauge fields, non-Hermitian trapping, and self-gravity

Recent work has generalized self-consistent TF theory to systems with density-generated gauge fields. For 2D abelian anyons in the magnetic-gauge picture, restricting to Slater determinants and replacing point fluxes by a density-sourced vector potential

$\psi^2$0

produces a self-consistent magnetic field

$\psi^2$1

The corresponding magnetic Thomas–Fermi functional is

$\psi^2$2

with

$\psi^2$3

and minimization gives

$\psi^2$4

This captures the qualitative energy trends of the Hartree/Chern–Simons–Schrödinger model, including the subtle dependence on the statistics parameter $\psi^2$5 (Levitt et al., 18 Apr 2025).

In a different direction, the stationary Gross–Pitaevskii equation with a harmonic real potential and a linear imaginary PT-symmetric term also admits a self-consistent TF limit. After amplitude–phase decomposition $\psi^2$6 and the rescaling $\psi^2$7, the leading-order TF system in $\psi^2$8 is

$\psi^2$9

For sufficiently small $Y_p$00, this truncated problem has a unique positive solution, the TF radius remains $Y_p$01, and the edge layer is governed by the Painlevé II equation with the Hastings–McLeod solution (Gallo et al., 2014). Self-consistency there couples density, phase gradient, and effective potential rather than density and Coulomb fields.

Not all such constructions are physically reliable. For self-gravitating Bose–Einstein condensates governed by the Gross–Pitaevskii–Poisson system, the naive TF approximation drops the kinetic term and yields a compact density

$Y_p$02

However, evaluating the full energy on that profile shows that the kinetic energy diverges logarithmically at the surface and that the TF total energy is positive, whereas a bound self-gravitating configuration should have negative total energy. An alternative smooth profile proportional to $Y_p$03 was therefore proposed for numerical work (Toth, 2014).

6. Numerical strategies, accuracy, and limitations

The numerical realization of self-consistent TF theory varies sharply with domain. Atomic polymer-based orbital-free DFT has been solved with Matsen’s bilinear Fourier method using zeroth-order spherical Bessel functions and a $Y_p$04 squared-$Y_p$05-norm convergence criterion; relativistic Wigner–Seitz matter has been treated on radial grids of up to 1201 points with composite Simpson integration; the 2D dipolar TFvW equation has been propagated in imaginary time with split operators, FFT, and Crank–Nicolson updates; and the anyon Hartree problem underlying magnetic TF was solved numerically with plane-wave spectral discretization and a Riemannian L-BFGS method on the Stiefel manifold (Thompson, 2022, Zhang et al., 2014, Zyl et al., 2012, Levitt et al., 18 Apr 2025).

Domain	Numerical strategy	Representative outcome
Atomic orbital-free Pauli model	Bilinear Fourier / spherical Bessel basis	Shell structure for H–Ar
Relativistic Wigner–Seitz matter	Radial ODEs on fine grids	Droplet, bubble, and uniform phases compared
2D dipolar TFvW gas	Imaginary-time split operator	Smooth edge replacing TF cutoff
Magnetic anyons	Plane waves + Riemannian L-BFGS	mTF reproduces qualitative energy trends

Performance is correspondingly domain-specific. The polymer-based atomic Pauli model attains correct shell structure and binding-energy errors below $Y_p$06 for the lightest elements and $Y_p$07 or less from oxygen through argon (Thompson, 2022). In non-uniform nuclear matter, STF and parameterized TF often differ mainly through surface treatment rather than bulk thermodynamics (Zhang et al., 2014). In anyonic systems, the magnetic TF density profile differs only subtly in real space from ordinary TF, but the energy retains the predicted oscillatory dependence on $Y_p$08 through $Y_p$09 (Levitt et al., 18 Apr 2025).

The limitations are equally systematic. Conventional TF and simple TFDW-type functionals do not produce atomic shell structure; the Englert–Schwinger model fails for low-$Y_p$10 atoms; relativistic STF omits shell effects, spin–orbit splittings, pairing, and, in some implementations, non-spherical pasta geometries; the finite-temperature TFS correction breaks when L-shell ionization invalidates the assumption of strongly bound inner electrons; and the self-gravitating BEC example shows that a self-consistent stationary solution can still be globally inconsistent as a bound state (Lehtomäki et al., 2017, Li et al., 31 Jul 2025, Segev et al., 2014, Toth, 2014). This suggests that “self-consistent Thomas–Fermi approximation” is best understood not as a single model, but as a family of semiclassical closure schemes whose accuracy depends on how local kinetic physics, nonlocal fields, and discrete strongly bound states are partitioned.