# density matrices for few identical fermions, harmonically trapped

## iibewegung

### Diagrammatic expansion

According to Feynman, density matrix for a 1D harmonic oscillator for ( T>0 ) is
[
\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}
]
where ( \beta \equiv 1/k*\text{B} T ).

The translationally invariant distribution of $ (x'-x) $ is
[
\tilde{\rho} (s; \beta) = \int*{\Omega*\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}}
]
where the integral is over all space ("( \Omega*\infty )").

The Fourier transform of ( \tilde{\rho} (s ; \beta) ) is the momentum distribution of the system, which is also a Gaussian.

Now we want to extend this to a system of 2 or 3 non-interacting identical fermions in the same harmonic potential. Following the textbooks, the many-particle density matrix is
[
\rho*{N*} (\mathbf{R}, \mathbf{R}' ; \beta ) = \sum*n \mathrm{e}^{-\beta E*n} \Psi*n ^* (\mathbf{R}) \Psi*n (\mathbf{R}')
]
where ( \mathbf{R} \equiv \left{ x*1 , x*2 , \cdots , x*N \right} ) are the many-body coordinates and ( \Psi *n (\mathbf{R}) ) is a generic many-body eigenfunction which works for all statistics (fermions, bosons, boltzmannons). In terms of the single-particle density matrix, the many-boltzmannon (Distinguishable particles) density matrix is
[
\rho*\text{D} (\mathbf{R}, \mathbf{R}' ; \beta ) = \prod*{j=1} ^N \rho*1 (x*j, x'*j; \beta) \; .
]
The many-fermion version is
[
\rho*{N\text{F}} (\mathbf{R}, \mathbf{R}' ; \beta ) = \cfrac{1}{N!} \sum*\mathscr{P} (-1)^\mathscr{P} \rho*\text{D} (\mathbf{R}, \mathscr{P} \mathbf{R}' ; \beta ) \equiv \sum*\mathscr{P} \sum*n \mathrm{e}^{-\beta E*n} \Psi*n ^* (\mathbf{R}) \Psi*n (\mathscr{P} \mathbf{R}')
]
where ( \mathscr{P} ) is a permutation operation (and the number of times it is done) with other identical particles and the sum is over all possible permutations.

To get a single-particle distribution of ( (x'-x) ) for two particles, we apply
[
\int*{\Omega*\infty} \mathrm{d}\mathbf{R} \mathrm{d}\mathbf{R}' \left[ \delta (s - (x'*1 - x*1)) \delta(x'*2 - x*2) * \; + \delta (x'*1 - x*1) \delta(s - (x'*2 - x*2)) * \; \right]
]
to ( \rho*{2\text{F}} ) to get
[
\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)
]
which is a composition of ( \tilde{\rho} (s; \beta) ) calculated above.

The case ( N=3 ) is a bit more complicated. We first apply
[
\int*{\Omega*\infty} \mathrm{d}\mathbf{R} \mathrm{d}\mathbf{R}' \left[ \delta (s - (x'*1 - x*1)) \delta(x'*2 - x*2) \delta(x'*3 - x*3) * \right. \; \
+ \delta (x'*1 - x*1) \delta(s - (x'*2 - x*2)) \delta(x'*3 - x*3) * \; \
\left. + \delta (x'*1 - x*1) \delta(x'*2 - x*2) \delta (s - (x'*3 - x*3)) * \; \right]
]
to ( \rho*{3\text{F}} ). Taking geometric restrictions into account (omitting terms of zero statistical weight corresponding to path-crossing, eg. ( (x*2 - x*1) (x'*2 - x'*1) < 0 ) ), we get
[
\tilde{\rho}*{3\text{F}} (s;\beta) = \frac{1}{3!} \left( 3 \tilde{\rho} (s; \beta) \tilde{\rho} (0; \beta)^2 - 4 \tilde{\rho} (s; 2\beta) \tilde{\rho} (0; \beta) + 2 \tilde{\rho} (s; 3\beta) \; \right) \; .
]

### Grand canonical, large-\(N\) approximation

In the large and fluctuating (N) and low-(T) limit, we can use the discrete Fermi function (f*n \equiv \frac{1}{\mathrm{e}^{\beta((n+\frac{1}{2})\hbar \omega - \mu)}+1}) (since we know the eigenstates of a single-particle in a harmonic trap) to construct
[
\rho*{N\text{F}} (x, x' ; \beta ) = \sum*n \cfrac{H*n (\sqrt{\frac{m\omega}{\hbar}} x) H*n (\sqrt{\frac{m\omega}{\hbar}} x') \mathrm{e}^{- \frac{m\omega(x^2 + x' ^2)}{2\hbar}}}{2^n n!} \sqrt{\frac{m\omega}{\pi \hbar}} f*n
]
where (\mu) was calculated numerically for each (N = \sum*n f*n), and (H_n (x)) are the Hermite polynomials.

As above, the distribution of (x'-x) is calculated as
[
\tilde{\rho}*{N\text{F}} (s;\beta) = \int*\infty ^\infty \, \mathrm{d}x \mathrm{d}x' \, \delta(s-(x'-x)) \rho_{N\text{F}} (x, x' ; \beta ) \; .
]

### Checking with Path integral Monte Carlo

Set (T=0.5). For (N=1), PIMC lands right on top of the diagrammatics-based (\tilde{\rho}_{N\text{F}} (s;\beta)). The grand canonical (GC) approximation, expected to work for larger systems, is understandably off. For (N=2), all 3 calculations are slightly off from each other. For (N=3), PIMC and GC seem to be converging but diagrammatics is off. For (N=20), PIMC and GC seem to have converged. (didn't attempt diagrammatics) So it looks like PIMC and GC have the correct large-(N) limit but something seems wrong with the diagrammatics in the few-particle case. Any ideas why?