## density matrices for few identical fermions, harmonically trapped

### 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?

Wednesday, December 12, 2012 22:56:58 GMT