# 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 [ \rho1 (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) ) \rho1 (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 ) = \sumn \mathrm{e}^{-\beta En} \Psin ^* (\mathbf{R}) \Psin (\mathbf{R}') ] where ( \mathbf{R} \equiv \left{ x1 , x2 , \cdots , xN \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 \rho1 (xj, 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} \sumn \mathrm{e}^{-\beta En} \Psin ^* (\mathbf{R}) \Psin (\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 - x1)) \delta(x'2 - x2) * \; + \delta (x'1 - x1) \delta(s - (x'2 - x2)) * \; \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 - x1)) \delta(x'2 - x2) \delta(x'3 - x3) * \right. \; \ + \delta (x'1 - x1) \delta(s - (x'2 - x2)) \delta(x'3 - x3) * \; \ \left. + \delta (x'1 - x1) \delta(x'2 - x2) \delta (s - (x'3 - x3)) * \; \right] ] to ( \rho{3\text{F}} ). Taking geometric restrictions into account (omitting terms of zero statistical weight corresponding to path-crossing, eg. ( (x2 - x1) (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 (fn \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 ) = \sumn \cfrac{Hn (\sqrt{\frac{m\omega}{\hbar}} x) Hn (\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}} fn ] where (\mu) was calculated numerically for each (N = \sumn fn), 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, 12 December 2012 22:56 GMT