Thomas-Fermi-Dirac Theory

There are two ways to derive equations for Thomas-Fermi-Dirac theory. One way is to start from grand potential and derive all equations from it. The other way is to start with low level equations and build our way up. Will start with the former.

Top Down Approach¶

The full interacting grand partition function is:

We are approximating it with a trial grand potential that is non-interacting:

And the relation between the two is:

Proof¶

The given relation is an approximation to the exact grand potential Ω[β, μ] = −(1/β) ln Ξ, derived using a mean-field-like approach based on a first-order cumulant expansion (or equivalently, approximating the average of an exponential by the exponential of the average). This is a standard technique in statistical mechanics for approximating partition functions when using a trial Hamiltonian, particularly in theories like Thomas-Fermi-Dirac (TFD) or mean-field approximations where fluctuations beyond the trial system are neglected.

Step 1: Recall the exact grand partition function and grand potential¶

The exact grand partition function for the interacting system is:

where the trace is over all particle numbers N and states, H is the full interacting Hamiltonian (including kinetic energy, external potential V_en, and electron-electron interactions), μ is the chemical potential, β = 1/(kT), and N is the particle number operator.

The exact grand potential is then:

This is intractable to compute exactly for interacting systems, so we introduce a non-interacting trial system.

Step 2: Introduce the trial Hamiltonian and its grand partition function¶

Define a trial Hamiltonian H_trial that is non-interacting (i.e., a sum of single-particle terms, as given in the notes: H_trial = ∑_i (p_i²/2 + V(x_i)), where V(x) = V_en(x) + V_ee(x) + V_xc(x) is an effective mean-field potential). The trial system is solvable exactly, as it reduces to independent fermions in the potential V(x).

The trial grand partition function is:

and the trial grand potential is:

As shown in the notes, this can be explicitly computed in the semiclassical (Thomas-Fermi) approximation:

(Here, the factor of 2 accounts for spin degeneracy.)

Step 3: Rewrite the exact Ξ in terms of the trial system¶

To connect the exact and trial systems, rewrite the exact grand partition function using H = H_trial + (H - H_trial):

This can be expressed as:

where ⟨⋅⟩_trial denotes the expectation value in the trial ensemble:

Let ΔH = H - H_trial for brevity (ΔH captures the difference between the true interactions and the mean-field approximation in H_trial). Then:

Step 4: Take the logarithm to get the grand potential¶

The grand potential is:

This is still exact, but the term ln ⟨e^{-β ΔH}⟩_trial is difficult to compute.

Step 5: Apply the first-order cumulant approximation (mean-field approximation)¶

To approximate ln ⟨e^X⟩_trial where X = -β ΔH, use the cumulant expansion:

In the mean-field approximation (common in TFD and similar theories), we truncate at first order, neglecting higher-order cumulants (which represent fluctuations/correlations beyond the mean field):

This is equivalent to the Jensen’s inequality-based approximation ⟨e^X⟩ ≈ e^{⟨X⟩} (since ln is concave, this gives an upper bound, but in practice, it’s used as an approximation here).

Substituting X = -β ΔH:

Plug this back into the expression for Ω:

This is the desired relation. The equality in the notes is understood as this approximation (not exact equality), which becomes accurate in the limit where correlations (higher cumulants) are small, such as in high-density or semiclassical regimes relevant to TFD.

Additional Notes on the Approximation¶

• Validity: This first-order truncation is exact if ΔH is constant (no fluctuations), but in general, it’s a mean-field approximation. In TFD, it’s justified because the trial H_trial already includes mean-field effects via V(x), and the semiclassical kinetic energy approximation further simplifies things.

• Connection to variational principle: This relates to the Gibbs-Bogoliubov-Feynman inequality, which states Ω ≤ Ω_trial + ⟨H - H_trial⟩_trial (an upper bound on Ω). In the notes, the approximation treats it as equality for the purpose of deriving the TFD equations, with minimization over the trial parameters (e.g., V(x) or the density n_e(x)) to find the best approximation.

Evaluate $\langle H - H_{\text{trial}} \rangle_{\text{trial}}$ explicitly¶

Where the full Hamiltonian contains:

And the non-interacting trial Hamiltonian is for each particle (the total Hamiltonian is a sum of these):

where: $V({\bf x}) = V_{en}({\bf x}) + V_{ee}({\bf x}) + V_{xc}({\bf x})$.

Now, evaluate the correction $\langle H - H_{\text{trial}} \rangle_{\text{trial}}$. The true $H$ includes Hartree $U_H = \frac{1}{2} \sum_{i \neq j} \frac{1}{|\mathbf{r}_i - \mathbf{r}_j|}$ and exchange-correlation (approximated as LDA for simplicity). In the mean-field:

• Hartree: $\langle U_H \rangle_{\text{trial}} \approx E_{ee} = \frac{1}{2} \int n_e(\mathbf{r}) V_{ee}(\mathbf{r}) \, d^3 r$, and $\langle \sum_i V_{ee}(\mathbf{r}_i) \rangle_{\text{trial}} = \int n_e V_{ee} \, d^3 r = 2 E_{ee}$, so contribution: $E_{ee} - 2 E_{ee} = -E_{ee}$.

• Dirac exchange: $E_{xc} = \int \epsilon_{xc}(n_e(\mathbf{r})) n_e(\mathbf{r}) \, d^3 r$, with $\epsilon_{xc} \propto -n^{1/3}$. Then $V_{xc} = \frac{\delta E_{xc}}{\delta n_e} = \frac{4}{3} \epsilon_{xc}$, so $\langle \sum_i V_{xc}(\mathbf{r}_i) \rangle_{\text{trial}} = \int n_e V_{xc} \, d^3 r = \frac{4}{3} E_{xc}$, and contribution: $E_{xc} - \frac{4}{3} E_{xc} = -\frac{1}{3} E_{xc}$.

• External: Cancels exactly (included in both $H$ and $H_{\text{trial}}$).

In the variational approximation: $\Omega[\rho] = \mathrm{Tr} \rho (H - \mu N + \frac{1}{\beta} \ln \rho)$, using a non-interacting trial $\rho_{\text{trial}}$:

where T_s is exact in KS (orbital kinetic) or approximate in TFD (semiclassical).

For U (true Coulomb), the trial average $\langle U \rangle_{\text{trial}} = \frac{1}{2} \mathrm{Tr} \rho_{\text{trial}} \sum_{i \neq j} \frac{1}{r_{ij}}$. Since $\rho_{\text{trial}}$ is non-interacting (product of single-particle states), it factorizes to the Hartree term E_H, but to capture quantum exchange-correlation, we approximate:

where E_xc is added as a functional of the trial density n(r) (e.g., Dirac: $E_x = -\frac{3}{4} \left( \frac{3}{\pi} \right)^{1/3} \int n^{4/3} d^3 r$).

More rigorously, here is how it works using DFT (both KS and OF) for the full system at zero temperature, but the analogous derivation works for a finite temperature also:

In other words, both for the full Hamiltonian and for the trial Hamiltonian, we make some approximation to the kinetic term (TF in our case), and the difference is taken care of using the XC potential.

In the correction:

With $\int n v_H = 2 E_H$ and $\int n v_{xc} = \frac{4}{3} E_{xc}$ (for Dirac), it becomes $-E_H - (1/3) E_{xc}$.

Thus the correction that must be applied is:

We start with a grand potential for non-interacting fermions (they interact in terms of the mean field potential $V({\bf x})$, but not directly), and add to it the corrections above:

The potential $V({\bf x}) = V_{en}({\bf x}) + V_{ee}({\bf x}) + V_{xc}({\bf x})$ is the total potential that the electrons experience (it contains nuclear, Hartree, and XC terms) and $E_{ee}$ is the Hartree energy:

For simplicity, we assume here that $V_{xc}$ only contains the exchange of the homogeneous electron gas. For a general XC functional, the relation is nonlinear and one must simply numerically calculate the XC energy density $e_{xc}({\bf x})$ and calculate the XC energy using:

In our case here, we have $e_{xc} = {3\over4}V_{xc}({\bf x})$, which is only true for the exchange in homogeneous electron gas. Otherwise the relation is nonlinear.

The density is a functional derivative with respect to $\mu$:

By defining the function $\Phi(n_e({\bf x}))$:

we can express the grand potential using $n_e$ as follows:

Now we can calculate the free energy:

where we used the fact that $\mu = {1\over \beta} \Phi(n_e({\bf x})) + V({\bf x})$, i.e. the left hand side $\mu$ is a constant, thus the sum of the terms on the right hand side is also constant (even though the individual terms are not).

We can calculate the entropy $S=-\left(\partial\Omega\over\partial T\right)_{V,\mu}$ as follows:

The total energy is then equal to:

From which we can see that the kinetic energy $E_{kin}$ is equal to:

The relation between the total energy and free energy can be also written as:

But it gives the same result as we obtained above.

To determine the kinetic part of the free energy, we set all potentials equal to zero ($V({\bf x}) = V_{en}({\bf x}) = V_{ee}({\bf x}) = V_{xc}({\bf x}) = 0$) and obtain:

If the potentials are zero, then the pressure can be calculated from:

If the potentials are not zero, then one can calculate the pressure using:

Summary:

where:

and $\mu N$ is calculated as follows:

So $F_e$ can also be expressed as: