few fermions in a harmonic trap — position density matrix from diagrammatics

Question

I'm trying to calculate the momentum distribution of a 1D system of non-interacting identical fermions in a harmonic trap.

Given Feynman's answer (from his Statistical Mechanics book) for the position density matrix of a single trapped particle at $T>0$, $ \rho_1 (x, x'; \beta) = \sqrt{\cfrac{m \omega}{2 \pi \hbar \sinh (\beta \hbar \omega) }} \exp \left\{ \cfrac{-m \omega}{ 2 \hbar \sinh (\beta \hbar \omega) } \left[ (x^2 + x' ^2) \cosh (\beta \hbar \omega) - 2x x' \right] \right\} $ ,

the translationally invariant distribution of $ (x'-x) $ is $ \tilde{\rho} (s; \beta) = \int_{-\infty} ^\infty \mathrm{d}x \mathrm{d}x' \delta (s - (x' - x) ) \rho_1 (x, x'; \beta) = \cfrac{\mathrm{e}^{\frac{-m \omega s^2}{4\hbar} \coth \frac{\beta \hbar \omega}{2}}}{2 \sinh \frac{\beta \hbar \omega}{2}} $ .

Here, we define $\beta\equiv 1/(k_\text{B}T)$. The Fourier transform of $\tilde{\rho} (s ; \beta)$ is the momentum distribution of the system, which is also a Gaussian.

How would you construct a two-fermion version out of this? What about 3 fermions?

My first attempt was $ \tilde{\rho}_{2\text{F}} (s;\beta) = \frac{1}{2!} \left( 2 \tilde{\rho} (s; \beta) \tilde{\rho} (0; \beta) - 2 \tilde{\rho} (s; 2\beta) \right) $

using Feynman diagrams, keeping in mind the anti-periodicity $\rho_1 (x, x'; \beta + \beta) = -\rho_1 (x, x'; \beta)$ for fermions.

Visualizing the path along the imaginary time on a hypercylinder of circumference $\beta$ would look something like

Note here that diagrams with fermion paths intersecting each other have zero statistical weight.

I've written the full statement of my problem at http://mathb.in/1393 . Apparently my answer doesn't agree with two other numerical calculations... perhaps missing a normalizing factor in one of the terms.

Any comments appreciated.

Answer 1

Remark: This answer will include a series expansion of the result with a description of a method for its summation into a closed form (without doing the full calculation).

The position density function which is more known as the "Heat kernel" can be constructed from the position representation of the energy eigenfunctions as follows:

$\rho(x, x^{\prime}, \beta) = \sum_n \Psi_n(x) \overline{\Psi_n(x^{\prime})} e^{-\beta E_n}$

For the single harmonic oscillator

$\Psi_n(x) = \frac{1}{\sqrt{2^n\sqrt{\pi} n!}}e^{-\frac{x^2}{2}}H_n(x)$, and $H_n(x)$ are the Hermite polynomials and $E_n = \hbar \omega (n+\frac{1}{2})$.

The summation into a closed form can be performed (for example) using the generating function:

$H_n(x) = \frac{n!}{2\pi i}\oint \frac{e^{(2tx-t^2)}}{t^{n+1}} dt$

Now in the case of two identical fermionic oscillators, the wave functions are given by:

$\Psi_{m,n}(x_1, x_2) = \frac{1}{\sqrt{2}} (\Psi_m(x_1) \Psi_n(x_2) -\Psi_m(x_2) \Psi_n(x_1))$, corresponding to the energies

$E_{m,n} = \hbar \omega (m+n+1)$

Which give the Heat kernel of the two identical fermionic oscillators upon the insertion in the Heat kernel expansion. Here also, the series can be summed using the generating functions although with more work.