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

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.

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.

• Sorry for the long delay in thanking you. I have just got around to working plotting your Ψ with Mathematica ("Slater") and found it to agree exactly with the "diagrammatic" version. The 2 curves lie right on top of each other. So it looks like there was a bug in my path integral Monte Carlo code ... Thanks again! – DCY May 11 '13 at 23:12
• I actually tried the "Slater" version with N=3 too. Everything except the "diagrammatic" version seem to agree fairly closely. I guess for N=3, my diagrammatic reasoning is somehow wrong. – DCY May 12 '13 at 3:59