### Ideal Quantum Gases

Introduction

The program implements a Monte Carlo simulation of the ideal gas in 1, 2, or 3 dimensions in equilibrium with a heat bath at temperature T. We assume that the single particle energies have the form ε = kp where k is the wavevector. (We will choose units such that the momentum and the wavevector k are the same.) The user can choose Fermi-Dirac (FD), Bose-Einstein (BE), or Maxwell-Boltzmann (MB) statistics. The output includes plots of the following quantities:

- The density of states (the number of states per energy level).
- The particle distribution function f(ε), the mean number of particles in a particular microstate with energy ε.
- The energy occupancy diagram which shows the energy of a state versus particle label.
- A k-space plot that shows which sites in phase space are occupied.

The energy of a particle in three dimensions is given by

We will choose units such that the prefactor h2/(2mL2) is equal to unity. We adopt periodic boundary conditions so that the quantum numbers are the positive and negative integers. For simplicity, we will take the ground state to be (0,0,0) in order that there be only one ground state. Note that the ground state has zero energy. Then the six first excited states are (1,0,0), (-1,0,0), (0,1,0), (0,-1,0),(0,0,1), and (0,0,-1).

Algorithm

According to quantum mechanics, particles of the same type are indistinguishable, and thus swapping two particles does not change the state of the system. To enforce this requirement the particles are labeled and the states are ordered according to their energy. Trial moves are made such that the order of the particles on this list cannot change. For example, particle number 3 is always in an energy state that is the same (possible for BE) or lower in the list than particle number 4. The Pauli exclusion principle for FD statistics is enforced by allowing only one particle per state (spin is ignored). [xx stupid Q: why can't this approach be used for interacting electrons to avoid the sign problem? xx] The only MC moves are for a particle to move up or down one state in the list of states. The algorithm is given by the following steps:

1. Set up a discrete k-space in 1, 2, or 3 dimensions with maximum wavenumber kmax. Order each state in this space by its energy.
2. Initially place the N particles in the lowest energy configuration.
3. Choose a particle at random and randomly move it up or down in the list. If the move is allowable, that is, if the state is unoccupied for FD statistics and if the order of the particles is unchanged, then compute the trial change in the energy ΔE that would occur if the particle were moved.
4. Use the Metropolis algorithm to accept or reject a move. If ΔE < 0, accept the move. If not, then if exp(-ΔE/T) > r, where r is an uniform random number between 0 and 1, accept the trial move; otherwise the move is rejected.
5. Repeat steps 3 and 4 for N attempted moves (one Monte Carlo step per particle (mcs).
6. Accumulate data for the mean number of particles in each state, the total energy, and the square of the total energy (the latter two quantities will be used to compute the specific heat).

Questions

1. Initialize the program using the default parameters. First focus on the density of states window. If kmax is large, we expect the density of states, g(ε), to be proportional to ε(d/p - 1). Try various values of the spatial dimension d and exponent p, and determine whether your results are consistent with this relation. If your results are not smooth, you can change the energy interval Δε and choose a larger value of kmax if your computer is sufficiently fast. Also you can change the steps per display to 100 or 1000 to speed up the simulation.
2. Does g(ε) depend on the statistics? Does it depend on the temperature? Does g(ε) change with the number of Monte Carlo steps per particle? Explain your answers.
3. Now focus on the energy occupancy diagram. Note that the width of each horizontal line gives the number of particles with a particular energy. Use the default values of the parameters and run the simulation for FD, BE, and then MB statistics. How does the energy occupancy diagram change with the number of Monte Carlo steps per particle? Explain the differences for the three cases. Why is it possible for more than one particle to have the same energy? As discussed, there is only one ground state and 6 first excited states with the same energy in three dimensions. How many states are there in the next energy level (ε = 2)? Repeat for p = 1 which is relevant for photons or highly relativistic particles. Do your conclusions change?
4. The program outputs the phase space diagram whenever the user presses stop. It is not continuously updated during the simulation because it is computationally slow. What is the shape of the region in phase space that is occupied for FD statistics in 1, 2, and 3 dimensions? What happens as the temperature is increased? How does the occupied region vary between the three statistics at low temperatures and at high temperatures?
5. Now consider the distribution function plot of n(ε). The results for n(ε) are averaged over all microstates with the same energy. In principle each data point at the same energy should be on top of each other, but in practice to do so would take a very long simulation. [xx why? if you are showing the average, there should be only a point. that would be less confusing xx] Use the default parameters with FD, but start with temperature T = 1. Describe the distribution function n(ε) and explain why it has the form it does. Run for thousands of mcs (displayed on the phase space diagram), and reset the data after a few thousand mcs to equilibrate the system. What happens if the temperature is increased? Now take more accurate data for the energy and specific heat at T = 0.5 and T = 1.0. Note the width of n(ε). Is it roughly proportional to T for low temperatures? If so, the specific heat should be proportional to T because this width is a measure of the change in energy of particle during a trial move. Compare ΔE/ΔT as an estimate for the specific heat with the two values you obtained at T = 0.5 and 1.0. They should be comparable. The program estimates the specific heat from the fluctuations in the energy.
6. Now consider BE statistics. Start at T = 1. What do you notice? At very low temperatures the ground state contains most of the particles. At what temperature does the ground state cease to be have a significant fraction of the particles (say less than 10%)? In the limit of an infinite system this temperature would correspond to a phase transition. Plot the specific heat as a function of temperature at low temperatures. The peak in the specific heat occurs at the phase transition in the thermodynamic limit.
7. Repeat the simulations for BE statistics for d = 1 and 2. In the thermodynamic limit there is no phase transition in one and two dimensions.
8. Consider BE statistics with p = 1. This value of p is equivalent to a gas of photons. In such a system the number of photons is not fixed as it is in our simulation. However, we can consider photons that reach the ground state as being lost to the system. The energy distribution function is what would be the blackbody spectrum. In the thermodynamic limit the location of the peak in this spectrum is proportional to the temperature. How well does the simulation produce this result for a finite systems?
9. The energy distribution function is the number of particles with energy ε between ε and ε + Δε. [xx correct? xx] Explain why this function is proportional to the density of states times the distribution function n(ε). This function is not plotted. [xx I think it is is. what is $\Delta \epsilon? xx] Is there any difference in the energy distribution function between the different statistics at low temperatures or at high temperatures?
10. The chemical potential μ can be extracted from the (equilibrium) particle distribution function n(ε). For FD statistics μ is the energy where the distribution function n(ε) equals 1/2. At T = 0 this energy is called the Fermi energy εF, and is where the step function drops from unity to zero. For low temperatures T << TF, μ - εF is proportional to T2. Collect data for μ between T = 0 and 4, and determine how your results compare.
11. For BE statistics the chemical potential μ can be found from the value of the ground state occupancy n(ε = 0), and is given by μ/T = -ln[1 + 1/n(0)] for temperatures above the Bose-Einstein condensation temperatures Tc. Note that μ is negative. How does μ behave near and above Tc? Would it be negative for FD or MB statistics? Is there a difference [xx difference in what? xx] at high temperatures (where quantum effects are negligible) and low temperatures?

References

- Harvey Gould, Lynna Spornick, and Jan Tobochnik, Thermal and Statistical Physics Simulations, John Wiley & Sons (1995), Chapter 5.

Java Classes Used

- QMC
- QMCApp

Updated 10 February 2007.