Energy Equilibration between a Lattice and Electrons

Introduction

            We have already seen an animation of energy equilibration between two gases of different mass and different average initial energy. In the present animation we will see how energy is equilibrated between thermal energy of relatively free electrons and a lattice which is made up of massive nuclei (we'll call that particle an ion for brevity) that are harmonically bound to periodic lattice points and have zero initial energy.

Figures

Figure 1: Showing the electrons as blue randomly distributed and the ions as red situated at lattice points.  Also shows the histograms as well as the least squares fits of the electron and ion and sum distribution.

Calculations

Note: In the following, we will use lower case symbols for the free mass's parameters and upper case symbols for the bound mass's parameters.

1. One dimensional calculations:

As a warm-up, we will first compute the conservation of momentum and energy for a one dimensional free electron and anchored nucleus.  The equation for momentum are:

$\begin{array}{l}mv\text{'}+MV\text{'}=mv+MV\\ p\text{'}+P\text{'}=p+P\end{array}$

Since the total momentum after the collision is the same as that before the collision, we can simplify the calculations by letting δp represent the change of momentum for each mass.             

$\begin{array}{l}p\text{'}=p+\delta p\\ P\text{'}=P-\delta p\end{array}$      

and those for energy are

$\begin{array}{l}\frac{1}{2}\left(mv{\text{'}}^2+MV{\text{'}}^2\right)=\frac{1}{2}\left(m{v}^2+M{V}^2\right)\\ \frac{1}{2}\left(\frac{p{\text{'}}^2}{m}+\frac{P{\text{'}}^2}{M}\right)=\frac{1}{2}\left(\frac{p^2}{m}+\frac{P^2}{M}\right)\end{array}$

$\frac{1}{2}\left(\frac{{(p+\delta p)}^2}{m}+\frac{{(P-\delta p)}^2}{M}\right)=\frac{1}{2}\left(\frac{p^2}{m}+\frac{P^2}{M}\right)$

$\frac{1}{2}\left(\frac{(p^2+2p\delta p+\delta p^2)}{m}+\frac{(P^2-2P\delta p+\delta p^2)}{M}\right)=\frac{1}{2}\left(\frac{p^2}{m}+\frac{P^2}{M}\right)$

$\frac{2p\delta p}{m}-\frac{2P\delta p}{M}=-\delta p^2\left(\frac{1}{m}+\frac{1}{M}\right)$

$\delta p=\left(\frac{2P}{M}-\frac{2p}{m}\right)\frac{mM}{m+M}=2\frac{mM}{m+M}(V-v)$

The displacement of M with respect to its rigid anchor point is expressed by

$X(t)=A\cos (\omega t+{\theta }_0)$                                                 

where A is the amplitude just prior to collision.

The momentum of M is also expressed as a function of time in the following manner:

$P(t)=-M\omega A\sin (\omega t+{\theta }_0)$

Immediately following the collision, the displacement of M has not changed but the momentum and energy have changed.  We can write:

$P\text{'}(t)=-M\omega A\text{'}\sin (\omega t+{\theta }_0+\delta \theta )=MV\text{'}=P-\delta p$                                (0)

where A and δθ are presently unknown.  The new amplitude A' can be found from the energy conservation equation:

$\frac{1}{2}kA{\text{'}}^2=\frac{1}{2}k{A}^2+\frac{1}{2M}(P{\text{'}}^2-P^2)$

where k is the constant of the linear restoring force, F(x):

$F(x)=-kX$

To obtain δθ, we use the fact that, immediately after the collision the displacement must be the same as just before the collision:

$X\text{'}(t_{c+})=A\text{'}\cos (\omega t_{c+}+{\theta }_0+\delta \theta )=A\cos (\omega t_{c-}+{\theta }_0)=X(t_{c-})$                           (1)

Getting both the sin and cosine terms in ωt+θ₀+δθ we can have a fail-safe equation for δθ:

$-M\omega A\text{'}\sin (\omega t_c+{\theta }_0+\delta \theta )=P-\delta p$

$A\text{'}\cos (\omega t_{c+}+{\theta }_0+\delta \theta )=A\cos (\omega t_{c-}+{\theta }_0)$

$\omega t_{c+}+{\theta }_0+\delta \theta =\text{atan2}\left(-\frac{P-\delta p}{M\omega A\text{'}},\frac{A}{A\text{'}}\cos (\omega t+{\theta }_0)\right)$

 

Using equation 1 to express A' in terms of A we have:

$A\text{'}\cos (\omega t_{c+}+{\theta }_0+\delta \theta )=A\cos (\omega t_{c-}+{\theta }_0)$

$A{\text{'}}^2[1-{\cos }^2(\omega t_{c+}+{\theta }_0+\delta \theta )]=A{\text{'}}^2-A^2{\cos }^2(\omega t_{c-}+{\theta }_0)=A{\text{'}}^2{\sin }^2(\omega t_{c+}+{\theta }_0+\delta \theta )$

Let's look at the difference P'²-P² at the time of the collision:

$P\text{'}{(t)}^2-P{(t)}^2=M^2{\omega }^2[A{\text{'}}^2{\sin }^2(\omega t_c+{\theta }_0+\delta \theta )-A^2{\sin }^2(\omega t_c+{\theta }_0)]$

$\begin{array}{l}P\text{'}{(}^t-P{(}^t=M^2{\omega }^2[A{\text{'}}^2-A^2{\cos }^2(\omega t_c+{\theta }_0)-A^2(1-{\cos }^2(\omega t_c+{\theta }_0))]\\ =M^2{\omega }^2(A{\text{'}}^2-A^2)\end{array}$

which should have been expected.

This is another expression of the energy change due to collision.  Notice that the difference is independent of the time of the collision.  Let's set this equal to the kinematic energy difference:

$M^2{\omega }^2(A{\text{'}}^2-A^2)=[{(P-\delta p)}^2-P^2]=[\delta p^2-2P\delta p]$

which gives us another expression for A'.  This completes the analysis for 1 dimension.

 

 

2. Three dimensional calculations:

            Here we will consider spherical particles which have the different masses, m and M, and diameters, D₁ and D₂.  The centers of the free/anchored spheres will be labeled (x,y,z) and (X,Y,Z).  Upon collision, the momentum transferred between the spheres will always be along the unit vector between their centers:

$u=\frac{\left[(x-X)\widehat{x}+(y-Y)\widehat{y}+(z-Z)\widehat{z}\right]}{r_{12}}$                                                   (2)

where

$r_{12}=\sqrt{{(x-X)}^2+{(y-Y)}^2+{(z-Z)}^2}$

is the distance between centers.  Since the animation is illustrated in only 2 dimensions, the collision analysis will assume a containing box that is large in the x and y dimensions but very thin in the z dimension. The following vector mathematics is correct for either 2 or three dimensions.

The expression for the final momenta in terms of the initial momenta is:

$mv\text{'}+MV\text{'}=mv+MV$                                                                          (3)

where the apostrophe on the left side of the equations indicates the final velocities.  We know that the energies are conserved so

$\frac{m{v}_{}^{\text{'}2}+M{V}_{}^{\text{'}2}}{2}=\frac{m{v}_{}^2+M{V}_{}^2}{2}$                                                                            (4)

The directions of the change in momenta are along the vector of centers, u, and the values of the changes of momenta must be equal and opposite.

$m_1\Delta v_1\equiv \mu \delta vu=-m_2\Delta v_2$                                                                     (5)

where μ has units of mass and is still to be determined.  Using equation 5 in equation 3:

$\begin{array}{l}mv\text{'}=mv+\mu \delta vu\\ MV\text{'}=MV-\mu \delta vu\end{array}$                                                                               (6)

Now we can use equation 6 in equation 4 to solve for the value of Mδv.

 $\frac{(mv+\mu \delta vu)•(mv+\mu \delta vu)}{2m}+\frac{(MV-\mu \delta vu)•(MV-\mu \delta vu)}{2M}=\frac{m{v}_{}^2+M{V}^2}{2}$ (7)

where the large dot stands for the dot product and equation 7 simplifies to:

$\begin{array}{l}\frac{(2m\mu \delta vu•v+{\mu }^2\delta v^2)}{2m}+\frac{(-2M\mu \delta vu•V+{\mu }^2\delta v^2)}{2M}=0\\ {\mu }^2\delta v^2\left(\frac{1}{m}+\frac{1}{M}\right)+2\delta v\mu (u•v-u•V)=0\\ \mu \delta v=\frac{2mMu•(V-v)}{m+M}\end{array}$                                    (8)

We can now make the identification:

$\mu =\frac{2m_1{m}_2}{m_1+m_2}=2{\mu }_r$                                                                               (9)

where μ_r=μ/2 is known as the "reduced mass".

Equations 6 and 8 are a complete solution for the final momenta. The final velocities are computed by dividing both sides of equations 6 by their respective masses:

 $\begin{array}{l}v\text{'}=v+\frac{2M}{m+M}u•(V-v)u\\ V\text{'}=V-\frac{2m}{m+M}u•(V-v)u\end{array}$                                      (10)

Viewed as a harmonic oscillator, we have the following equations for displacement of M at the time of the collision:

$\begin{array}{l}X=A_x\cos (\omega t+{\theta }_{0x})=A{\text{'}}_x\cos (\omega t+{\theta }_{0x}+\delta {\theta }_x)\\ Y=A_y\cos (\omega t+{\theta }_{0y})=A{\text{'}}_y\cos (\omega t+{\theta }_{0y}+\delta {\theta }_y)\\ Z=A_z\cos (\omega t+{\theta }_{0z})=A{\text{'}}_z\cos (\omega t+{\theta }_{0z}+\delta {\theta }_z)\end{array}$           (11)

where the result on the right is for just after the collision.  For convenience I will assume that the restoring force is the same in all three directions so there is only one ω.  Since A_x and θ₀ are known, Equations 11 can be solved for either A' or δθ.

 

The velocity components in all three directions can be written from equation 11:

$\begin{array}{l}{V}_x=-\omega A_x\sin (\omega t+{\theta }_{0x})V{\text{'}}_x=-\omega A{\text{'}}_x\sin (\omega t+{\theta }_{0x}+\delta {\theta }_x)\\ V_y=-\omega A_y\sin (\omega t+{\theta }_{0y})V{\text{'}}_y=-\omega A{\text{'}}_y\sin (\omega t+{\theta }_{0y}+\delta {\theta }_y)\\ V_z=-\omega A_z\sin (\omega t+{\theta }_{0z})V{\text{'}}_z=-\omega A{\text{'}}_z\sin (\omega t+{\theta }_{0z}+\delta {\theta }_z)\end{array}$

Just as for the one dimensional case, we will equate the energy changes in the harmonic oscillator case with the energy changes in the hard sphere collision case as given below:

$MV\text{'}=MV-\mu \delta vu$

$\frac{M^{2}V{\text{'}}^2-M^2{V}^2}{2M}=\frac{1}{2M}\left({(\mu \delta v)}^2-2\mu M\delta v(V•u)\right)=\frac{1}{2}\left(k_x(A{\text{'}}_x^2-A_x^2)+k_y(A{\text{'}}_y^2-A_y^2)+k_z(A{\text{'}}_z^2-A_z^2)\right)$ Similarly we can write the components of V'-V in terms of the harmonic oscillator velocity components:

$\begin{array}{l}{V}_x\text{'}-V_x=-\frac{2m}{m+M}u•(V-v)u_x=-\omega \left[A_x\text{'}\sin (\omega t+{\theta }_{0x}+\delta {\theta }_x)-A_x\sin (\omega t+{\theta }_{0x})\right]\\ V_y\text{'}-V_y=-\frac{2m}{m+M}u•(V-v)u_y=-\omega \left[A_y\text{'}\sin (\omega t+{\theta }_{0y}+\delta {\theta }_y)-A_y\sin (\omega t+{\theta }_{0y})\right]\\ V_z\text{'}-V_z=-\frac{2m}{m+M}u•(V-v)u_z=-\omega \left[A_z\text{'}\sin (\omega t+{\theta }_{0z}+\delta {\theta }_z)-A_z\sin (\omega t+{\theta }_{0z})\right]\end{array}$       (12)

Equations 12 can be solved for either A' or δθ, which ever has not been obtained from solution of equations 11.

$\begin{array}{l}{A}_x\text{'}\sin (\omega t+{\theta }_{0x}+\delta {\theta }_x)=-\frac{V_x\text{'}-V_x}{\omega }+A_x\sin (\omega t+{\theta }_{0x})\equiv s_x\\ A_y\text{'}\sin (\omega t+{\theta }_{0y}+\delta {\theta }_y)=-\frac{V_y\text{'}-V_y}{\omega }+A_y\sin (\omega t+{\theta }_{0y})\equiv s_y\\ A_z\text{'}\sin (\omega t+{\theta }_{0z}+\delta {\theta }_z)=-\frac{V_z\text{'}-V_z}{\omega }+A_z\sin (\omega t+{\theta }_{0z})\equiv s_z\\ A_x\text{'}\cos (\omega t+{\theta }_{0x}+\delta {\theta }_x)=A_x\cos (\omega t+{\theta }_{0x})\equiv c_x\\ A_y\text{'}\cos (\omega t+{\theta }_{0y}+\delta {\theta }_y)=A_y\cos (\omega t+{\theta }_{0y})\equiv c_y\\ A_z\text{'}\cos (\omega t+{\theta }_{0z}+\delta {\theta }_z)=A_z\cos (\omega t+{\theta }_{0z})\equiv c_z\end{array}$      (12)

Using the definitions s_i and c_i and the atan2(y,x) function, we have fail-safe solutions for δθ:

 

$\begin{array}{l}\delta {\theta }_x=a\tan 2\left(s_x,c_x\right)-(\omega t+{\theta }_{0x})\\ \delta {\theta }_y=a\tan 2\left(s_y,c_y\right)-(\omega t+{\theta }_{0y})\\ \delta {\theta }_z=a\tan 2\left(s_z,c_z\right)-(\omega t+{\theta }_{0z})\end{array}$

Similarly we obtain the solutions for A':

$\begin{array}{l}{A}_x\text{'}=\frac{A_x\cos (\omega t+{\theta }_{0x})}{\cos (\omega t+{\theta }_{0x}+\delta {\theta }_x)}\\ A_y\text{'}=\frac{A_y\cos (\omega t+{\theta }_{0y})}{\cos (\omega t+{\theta }_{0y}+\delta {\theta }_y)}\\ A_z\text{'}=\frac{A_x\cos (\omega t+{\theta }_{0z})}{\cos (\omega t+{\theta }_{0z}+\delta {\theta }_z)}\end{array}$

Summary

            The result of the animation is that the energy per particle obtained by the ions is greater than that obtained by the electrons.  A reasonable fit to the ion energy density distribution is:

$N_I(E)=a_{I}E\exp (-b_{I}E)$

where a_I and b_I  are parameters found by a 2-dimensional Newton-Raphson technique that minimizes the residual errors between N_I(E) and the ions per energy bin found by iteration of collisions.  Similarly the electron energy density must be fit by a distribution like the following:

$N_E(E)=\exp (-b_{E}E)(1-a_{E}E)$

It is presently believe that the reason that the ions end up with more energy per particle is that they can store any energy that they obtain in their restoring potentials.