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?

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

URL:

Code:

<h3>Diagrammatic expansion</h3>
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 (<span style="font-style:italic;">D</span>istinguishable 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) \; .
\]

<h3>Grand canonical, large-$N$ approximation</h3>
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 ) \; .
\]

<h3>Checking with Path integral Monte Carlo</h3>
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.
<a href="http://imgur.com/I4Don"><img src="http://i.imgur.com/I4Don.png" title="Hosted by imgur.com" alt="" width="100" /></a>
For $N=2$, all 3 calculations are slightly off from each other.
<a href="http://imgur.com/NJxfJ"><img src="http://i.imgur.com/NJxfJ.png" title="Hosted by imgur.com" alt="" width="100" /></a>
For $N=3$, PIMC and GC seem to be converging but diagrammatics is off.
<a href="http://imgur.com/roBcf"><img src="http://i.imgur.com/roBcf.png" title="Hosted by imgur.com" alt="" width="100" /></a>
For $N=20$, PIMC and GC seem to have converged. (didn't attempt diagrammatics)
<a href="http://imgur.com/6ef0c"><img src="http://i.imgur.com/6ef0c.png" title="Hosted by imgur.com" alt="" width="100" /></a>
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?

Title: Name:

Generate secret URL (Learn more)

MathB.in α

Share mathematics on the web