Energy Distribution of 3D Rotating and Vibrating Diatomic Molecular Gas Using A Lennard-Jones Type Potential for Molecular Binding
Here we use attractive forces that are mediated by a Lennard-Jones (LJ) type potential
as the binding force between the atoms of the diatomic molecule. The form of the LJ potential is
VLJ=4*epsilon*((sigma/r)^2b-(sigma/r)^b)
where epsilon is the minimum (negative) potential, r is the distance between atom centers and sigma is a scaling
length which will usually be approximately equal to r. The exponent b is chosen to be 6 for the usual LJ potential but can be
modified for either tighter binding or more vibrational motion. The LJ potential is strongly repelling at close range and weakly attracting at longer ranges.
For this simulation the power used was b=3.
From the LJ potential a force, FLJ, can by computed by taking the derivative with respect to r:
FLJ=-dVLJ/dr
and this force is used iteratively in Newton's laws of motion to update the momentum and speed of the mobile atoms.
When the atoms have initial rotation speed, v, about their common center point, the equation for a stable rotation radius becomes:
FLJ(r)+2mv*v/r=0
where r is the separation and r/2 is the rotation radius of the atoms. If we multiply equation 3 by r,
we can solve it for both the stable r as well as the maximum rotation speed.
If we set the initial separation to be less than or greater than this stable radius, then the motion of the atoms will have sinusoids
(i.e. the separation of the atoms will vibrate) superimposed on the stable radius.
For this animation we set sigma equal to twice the hard disc radius of the atoms.
For r>sigma we use the forces from the above equation to compute the speed changes.
The form of the kinetic energy distribution for simple atoms is, with 2 dimensions, a simple exponential:
fa(E)=a*exp(-E/Eavg)
where a is a normalization constant and Eavg is the average kinetic energy of the atoms.
The form of the kinetic energy distribution for the diatomic molecules (often called rotors), with
3 dimensional motion, turns out to be the following exponential form:
fd(E)=(5E/3Eavg)3/2*exp(-5E/(2Eavg))
where Eavg is the average kinetic energy of the rotors.
Conventional thermal physics would say that the difference in these two distribution is that we have add two additional degrees of freedom.
Actually we have added a single variable constraint, the binding energy between the two atoms of the rotor.
By adding this constraint we have increased the number of modes of energy exchange by two, vibration and rotation.
For the atoms, as well as the center of mass motion of the rotors, we have three modes of energy exchange, from motion along the x axis to motion along the y and z axes, etc.
For the rotors in 3D we now have 5 total modes of energy exchange. and that gives rise to the 3/2 factor,
(E/Eavg)^(3/2), multiplying their exponential energy distribution. The actual distribution depends somewhat on the initial kinetic
energy of the rotors. If the initial energy is high enough, the vibrational modes are excited and that leads to one more mode of
energy exchange. The Debye model (https://en.wikipedia.org/wiki/Debye_model) addresses the excitation of the vibrational modes.
If they are fully excited then the expression for the distibution is
fd(E)=E2*exp(-3E/Eavg)
which is plotted in blue.