In the classical theory of ideal gases, the velocity distribution function is derived from the Boltzmann transport equation based on the assumption of molecular chaos and a dilute enough gas so that ternary and higher collisions can be ignored  . In an ideal gas trajectories of the particles are straight lines, meaning that the particles interact through short range potentials, such as hard-sphere or Lennard-Jones potentials. Under these conditions, the velocity distribution of the particles is described by the Maxwell-Boltzmann distribution function.
However, it can be shown that the Maxwell-Boltzmann distribution function describes molecular velocities even in non-ideal gases   , even in the presence of long-range interactions between the particles of the system such as Coulomb potential, as long as the assumption of molecular chaos remains valid  . Not long ago, extremely dense fluid systems with strong interactions between its particles were investigated through molecular dynamics simulations. It was shown that even in these systems the velocity distribution of the particles is Maxwellian and the assumption of molecular chaos remains valid  .
In what follows, we extend this investigation to include solids with strong harmonic interactions between its particles. Solids are interesting because contrary to fluids where the particles can move throughout the entire volume available to the fluid, in solids the particles can only oscillate about their mean equilibrium positions, and for strong interactions these oscillations are restricted to a very limited region of space.
The choice of harmonic interactions between the particles is due to the fact that any interaction between particles of a solid can be approximated by a harmonic potential as long as the displacements of the particles from their equilibrium positions remain small. Of course, solids with harmonic interatomic interactions have long been studied extensively  -  . However, to the best of the author’s knowledge, velocity distribution of atoms in these solids, which is the subject of this investigation, has not been discussed or reported in the literature.
2. Theoretical Background
The normalized Maxwell-Boltzmann velocity distribution function (or probability density function) of a perfect gas in one dimension is given by 
Assuming molecular chaos, the velocity distribution function in three dimensions can be written as the product of the one-dimensional distribution functions   ,
where . Then the probability of finding a particle with velocity in the interval to is . Transforming this probability to spherical coordinates in velocity space with and integrating over q and f, we obtain the Maxwell-Boltzmann speed distribution function in three dimensions,
As stated above, in this work we assume a harmonic potential for interaction of the atoms. This is because any general interatomic potential energy function of two atoms can be expanded is a Taylor series about the equilibrium position of the atoms,
where r is the distance between the two atoms and r0 is the equilibrium distance. The first term on the right hand side is a constant, which defines the reference of the potential energy. This constant can arbitrarily be set equal to zero. The second term is zero because the force on the particle is given by
and at equilibrium the force on the particle is zero. The coefficient of the third term is positive because the potential energy is a minimum at a stable equilibrium point. Finally, for small displacements of the particles from their equilibrium ( ), we can ignore the higher-order terms in the expansion (4), and write
where is a positive constant. Therefore, for small deviations from equilibrium, a particle under any interatomic potential energy function in a solid experiences a harmonic potential. However, it should be pointed out that this analysis is valid as long as the potential energy function is not intrinsically nonlinear  .
The interactions between atoms in many simple solids are described by Lennard-Jones (6-12) potential energy function    ,
where and are constants with units of energy and length, respectively. The first term in this equation corresponds to attraction between the two atoms resulting from van der Waals interaction between them  . The second term corresponds to repulsion arising from overlap of the electron distributions of the atoms and the Pauli Exclusion Principle. This term becomes dominant at small distances between the atoms. The meaning of the constants and in Equation (7) can be inferred from Figure 1, which schematically shows the graph of Lennard-Jones potential energy function. The minimum of the potential energy corresponds to the equilibrium distance between two atoms .
Differentiating Equation (7) with respect to r and setting it equal to zero gives
A second differentiation of Equation (7) and its evaluation at gives
Figure 1. Lennard-Jones (6-12) potential energy function.
Therefore according to Equation (6), for small deviations from equilibrium, the Lennard-Jones potential reduces to a harmonic potential given by
3. Molecular Dynamics Simulations
We investigate both a one-dimensional and a three-dimensional solid. This is because in some cases one-dimensional systems behave differently from those of higher dimensions. In addition, correlating the results of the one-dimensional system with those of the three-dimensional system can further verify the validity of the assumption of molecular chaos. Furthermore, in order to mimic an actual atomic solid in our simulations, we consider argon-like particles interacting according to a harmonic potential energy function given by Equation (10).
Since the actual atomic parameters are very small in SI units, we use a reduced system of units. In this system the mass of the atoms m is the unit of mass and the equilibrium distance between the atoms r0 is the unit of length. For argon, the atomic mass and the Lennard-Jones potential parameters are   ,
Therefore, in our system of units, the unit of mass is and the unit of length (r0) is . Finally, we take and to be our units of energy and temperature, respectively, where is the Boltzmann’s constant. This completes our system of units. Any other physical quantity in our simulations can be written as a combination of these units, as shown in Table 1.
All of the simulations, in one and in three dimensions, were carried out using the velocity form of the Verlet algorithm   ,
Table 1. The system of units used in the molecular dynamics simulations in this work.
with a time step of reduced unit, where x stands for x, y, or z, and v and a are the corresponding velocity and acceleration, respectively. Furthermore, throughout the investigation periodic boundary conditions were used to reduce the size effect of the systems.
3.1. One-Dimensional Solid
We studied a linear chain of 200 argon-like particles, each of mass m, with spring-like interaction between the adjacent particles,
where r is the instantaneous distance between two adjacent atoms and r0 is the relaxed length of the spring. Since
according to Equation (10) the constant c in Equation (13) is given by
which, in our reduced system of units, has a value of 72.
The particles were initially arranged on a straight line at one unit of length intervals. Each particle was then randomly given a velocity of +1 or −1 reduced unit and then the center of mass of the system was brought to rest. The temperature of the system was calculated from
where is the mean kinetic energy of the particles. The velocity of each particle was then re-scaled according to the following equation to adjust the temperature of the system to the desired value  ,
where is the actual temperature calculated from (16) and is the desired temperature.
The simulation was allowed to run for 1000 time steps to thermalize the system before any data were collected. After thermalization, data were collected on the temperature of the system as a function of time step and the velocity distribution of the particles over the next 105 time steps. We also monitored the velocity of the center of mass of the system to ensure the center of mass did not drift during the simulation, or there was no “wind” in the system.
We simulated the system with temperatures of and reduced unit. Figure 2 shows the temperature of the system as a function of time step for reduced unit. As can be seen from this figure, the temperature of the system remained very steady throughout the simulation. The results for were similar.
Figure 2. Evolution of temperature as a function of time step in the one-dimensional simulation.
Figure 3. Velocity distribution functions for a one-dimensional chain of particles interacting according to a harmonic potential at two different temperatures. The markers are the simulation results and the solid curves are predictions from the Maxwell-Boltzmann distribution law.
Figure 3 shows the velocity distribution functions obtained from our simulations for and reduced unit. This figure also shows the one-dimensional Maxwell-Boltzmann velocity distribution functions for the same temperatures. These are obtained from Equation (1) by substituting the corresponding quantities in reduced units. Thus, for ,
and for ,
These are shown as solid curves in Figure 3. As can be seen, the simulation results are in excellent agreement with the Maxwell-Boltzmann distribution function for molecular velocities.
3.2. Three-Dimensional Solid
In this part of the investigation, we studies a system of 512 argon-like particles ( ) on a simple cubic lattice, the unit cell of which is shown in Figure 4. Each particle interacted harmonically with its 6 nearest neighbors according to Equation (13) with the same potential parameters described above. As in the one-dimensional case, periodic boundary conditions were used to reduce the size effect of the system.
Each particle was initially given a random velocity in the interval in each coordinate direction. These velocities were then reduced to bring the center of mass of the system to rest. The temperature of the system was then calculated from
as before, then the velocity components of each particle were re-scaled according to Equation (17) to adjust the temperature of the system to the desired value. The system was allowed to thermalize for 5000 time steps. Then data were collected during the next 105 time steps on temperature and speed distribution of the particles for and reduced unit.
Figure 5 shows the temperature of the system as a function of time step for reduced unit. As can be seen from the figure, similar to the one-dimensional case, the temperature of the system remained very uniform throughout the simulation. The result for was similar and very uniform.
Figure 6 shows the speed distribution functions obtained in the simulation of
Figure 4. The unit cell of the three-dimensional solid investigated. Each particle interacts harmonically with its 6 nearest neighbors.
Figure 5. Evolution of temperature as a function of time step in the three-dimensional simulation.
Figure 6. Speed distribution functions for a three-dimensional system of 512 particles interacting according to a harmonic potential at two different temperatures. The markers are the molecular dynamics simulation results and the solid curves are predictions from the Maxwell-Boltzmann distribution law.
the three-dimensional system at and reduced unit. The figure also shows the Maxwell-Boltzmann speed distribution functions at these two temperatures. These are obtained from Equation (3) by substituting the parameters of the system in reduced units. The results for and are, respectively
These are shown by solid curves in Figure 6. Again as in the one-dimensional case, the agreement between the molecular dynamics simulations and the Maxwell-Boltzmann distribution functions are excellent.
We also point out that in the three-dimensional case as well, each of the velocity distribution functions , , were similar to that of the one-dimensional system shown in Figure 3, and they were all in near perfect agreement with the one-dimensional Maxwell-Boltzmann distribution function (1).
4. Discussion and Conclusions
Figure 2 and Figure 5 show that during the entire simulation in one- and three-dimensional systems the temperature remained very stable. The large fluctuation at the beginning of the run in each case is due to small number of time steps used in the averaging process.
As pointed out earlier, the Maxwell-Boltzmann velocity distribution function, which was originally obtained for an ideal gas, remains valid even in non-ideal gases and in the presence of long-range interactions between the particles of the system as long as the assumption of molecular chaos remains valid. Later, molecular dynamics simulations showed that even in extremely dense fluid systems with strong interactions between its particles velocity distribution of the particles is Maxwellian and the assumption of molecular chaos remains valid.
In this study we have taken the investigation one step further to include solid systems where, contrary to fluid systems, motion of the particles is very limited in space, especially at low temperatures. For the argon-like particles studied, according to Table 1, the temperatures and reduced unit translate into 60 K and 120 K, respectively, which are quite low. At these low temperatures, the amplitude of the oscillations of the particles is very small and the particles remain in a very limited region of space. In fact, in all of our simulations, the mean deviation of the interatomic distance from its equilibrium value was only about 9% or less across the board. Nevertheless, the results clearly show that even in this case the velocity distribution of the particles remains Maxwellian.
Furthermore, since Equation (3) is obtained from Equation (1) based on the assumption of molecular chaos, the fact that the velocity and the speed distribution functions in one and three dimensions are both in excellent agreement with these equations is indicative of the validity of the assumption of molecular chaos in solids, even at low temperatures.
In addition, the Maxwell-Boltzmann distribution function is obtained from the Boltzmann transport equation which, in turn, is based on the assumption of a dilute gas  . According to this derivation, therefore, the system must be in the classical regime,
where n is the particle density and
is the quantum concentration of the system   . In this equation is the Plank’s constant, which in our reduced system of units has a value of 0.164 and unit of . For our three-dimensional system with simple cubic structure, there is one atom per unit cell and each side of the unit cell has a length r0, therefore, . Then with reduced units, we have
Therefore, even though the system that we have investigated is a solid at low temperatures, it is still in the classical regime.
 Talu, O. and Myers, A.L. (2001) Reference Potentials for Adsorption of Helium, Argon, Methane, and Krypton in High-Silica Zeolites. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 187-188, 83-92.