Maxwell-Boltzmann Distribution in Solids

The velocity distribution functions of particles in oneand three-dimensional harmonic solids are investigated through molecular dynamics simulations. It is shown that, as in the case of dense fluids, these distribution functions still obey the Maxwell-Boltzmann law and the assumption of molecular chaos remains valid even at low temperatures.


Introduction
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 [1].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 [2] [3], 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 [4].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 [4].
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 avail-able 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 [5]- [11].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.

Theoretical Background
The normalized Maxwell-Boltzmann velocity distribution function (or probability density function) of a perfect gas in one dimension is given by [12] ( ) Assuming molecular chaos, the velocity distribution function in three dimensions can be written as the product of the one-dimensional distribution func- , , exp 2π 2 where Then the probability of finding a particle with velocity in the interval ( ) , , Transforming this probability to spherical coordinates in velocity space with 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 r 0 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 ( 0 r r − ), 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 [13].
The interactions between atoms in many simple solids are described by Lennard-Jones (6-12) potential energy function [14] [15] [16], ( ) 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 [14].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 Differentiating Equation ( 7) with respect to r and setting it equal to zero gives A second differentiation of Equation (7) and its evaluation at Therefore according to Equation ( 6), for small deviations from equilibrium, the Lennard-Jones potential reduces to a harmonic potential given by ( ) ( )

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 r 0 is the unit of length.For argon, the atomic mass and the Lennard-Jones potential parameters are [15] [17], Therefore, in our system of units, the unit of mass is 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 [18] [19], ( ) ( ) Table 1.The system of units used in the molecular dynamics simulations in this work.reduced unit, where ξ 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.

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 r 0 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 K is the mean kinetic energy of the particles.The velocity of each par- ticle was then re-scaled according to the following equation to adjust the temperature of the system to the desired value [20], where a T is the actual temperature calculated from ( 16) and d T 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 10 5 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    1.0 T = , ( ) 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.

Three-Dimensional Solid
In this part of the investigation, we studies a system of 512 argon-like particles ( 8 8 8 × × ) 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 ( ) 1,1 − 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.
v v θ θ φ = and integrating over θ and φ, we obtain the Maxwell-Boltzmann speed distribution function in three dimensions, ( )

P
. Mohazzabi, S. P. Shankar DOI: 10.4236/jamp.2018.63052604 Journal of Applied Mathematics and Physics , 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 0 r .
Figure 2 shows the temperature of the system as a function of time step for 0.5 T = reduced unit.As can be seen from this figure, the temperature of the system remained very steady throughout the simulation.The results for 1.0 T = were similar.

Figure 2 .
Figure 2. Evolution of temperature as a function of time step in the one-dimensional simulation.

Figure 3 .
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
Figure 3 shows the velocity distribution functions obtained from our simulations for 0.5 T = and 1.0 T = 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 The system was allowed to thermalize for 5000 time steps.Then data were collected during the next 10 5 time steps on temperature and speed distribution of the particles for

Figure 5
Figure 5 shows the temperature of the system as a function of time step for 0.5 T = 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 1.0 T = was similar and very uniform.

Figure 6 Figure 4 .
Figure 6 shows the speed distribution functions obtained in the simulation of

Figure 5 .
Figure 5. Evolution of temperature as a function of time step in the three-dimensional simulation.

Figure 6 .
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.