The second law of thermodynamics states that, in a closed system, the motional (kinetic) energy will not have net spreading from a less energetic region to a more energetic region. The second law is a result of diffusion (spreading) and all diffusion is due to particle interaction. When the regions start out in separate locations, both low energy regions and high energy regions spread but, the higher energy regions spread faster so the net movement is from high energy regions to low energy regions. In this animation, the particle interaction is by way of hard sphere collisions but, in reality, there is no such thing as a hard sphere collision: In the real world, all collisions are mediated by potentials between the particles. Potential-mediated collisions are handled under other menu items.
We start with red spheres on the left and blue spheres on the right and observe how these mix as a function of x and time. Even if their radius can be different their mass and initial energy is the same. This mixing is a manifestation of position diffusion. The spheres can be chosen to be of different radii since this results in different self collision cross sections.

Continuity Equation for Position Diffusion

The continuity equation is fundamental to all of physics and is appropriate for changes in particle density which is the same as position diffusion:

`(dn)/dt=nabla cdot n bb v (1)`

where `n` is the number density of particles in a small volume located at position (`x,y,z`) and `bb v` is the average velocity of the particles in that small volume. Equation 1 is called the continuity equation. Please note that a bold character like ` bb u` signifies a vector that has (`x,y,z`) components (`u_x,u_y,u_z`). Some might recognize the term `n bb v` as the particle current, which is the particle flow per unit area per second. This differential equation just says that the rate `(dn)/dt` of flow of particles into the small volume is related to the slope of the density times the velocity of groups of particles in the neighboring vicinity of that small volume. The divergence operator `nabla cdot` has the followng effect on its vector argument:

`nabla cdot n bb v=(del (n v_x))/(del x)+(del (n v_y))/(del y)+(del (n v_z))/(del z) (1a)`

Of course both `n` and `bb v` can be variables of (`x,y,z`) so we must expand each of the 3 terms e.g. for the x component:

`(del n v_x)/(del x)=(del n )/(del x)v_x+n(del v_x )/(del x)`

If we could have enough particles so that each small volume is well populated we could solve equation 1 numerically but with the number of particles we can actually achieve, the results would be very coarse-grained and volatile.

Converting the Continuity Equation to the Diffusion Equation

Recall that we used the expression

`nabla cdot n bb v=(dn)/dt (1)`

to determine the rate of change of particle density in a small volume. Let us rename the `n bb v` and call it `bb C` which is similar to the current density, `bb J`, used in electromagnetics. Thus we have

`bb C=n bb v (2)`

For our case, the sign of any velocity component, say `v_x`, of any particle is random. The particle velocity is constant from one collision to the next and then, on average, the particles have a completely new velocity. So the change in position of the particles is limited by the mean distance they travel between collisions. This mean distance is called the mean free path `l`. The mean free path is defined by

`l=1/(n sigma)`

where `n` is the local average density and

`sigma=pi(r_{red}+r_{blue})^2`

and where `r` are the radii of the particles. The average time between collisions is then

`tau_c=l/|bar v|`

where

`bar v=sqrt(bar (v^2))=sqrt((2 bar E)/m)`

where `bar E` is the average kinetic energy. Thus the only asymmetry that results in the flow of groups of particles is the spatial variation of `n` (its gradient) so the rate of flow of groups of particles is

`bb C= bar v l nabla n`

Another way of saying this is to define the x component the group velocity as:

`v_x^(group)= 1/n bar v l (del n)/(del x)`

and that group velocity component is the what goes into equation (1a).

`nabla cdot n bb v=(del (n v_x))/(del x)+(del (n v_y))/(del y)+(del (n v_z))/(del z)= bar v l [del ((del n)/(del x))/(del x)+(del ((del n)/(del y)))/(del y)+(del ((del n)/(del z)))/(del z)]`

Then, in the x direction, the net rate of `C_x` is equal to

`C_x= bar v l(n(x+dx/2)-n(x-dx/2))/dx=bar v l(del n)/(del x)`

Then taking the divergence of `bb C` we have the the expression

`(dn)/dt=bar v l nabla cdot nabla n=bar v l(del^2 n )/(del x^2)`

Expanding this equation we have:

`(dn)/dt=D[( del^2 n)/(del x^2)+( del^2 n)/(del y^2)+( del^2 n)/(del z^2)] (3)`

where D is called the difusivity . D is made up of the product .

`D=l sqrt(bar (v^2))`

and has units `(meters^2)/(sec)`.
Equation 3 can be solved numerically using the Finite Difference Method. To simplify solving it, the initial n(x) will have no discontinous derivative. Even though we specify that one particle type will be on the left and the other on the right, the actual reduction in n(x) at the boundary will not be sharp because of the finite number of particles that we use. Thus, for analysis, the step in particle density at the center boundaries will be swaged and the slope of the particle density at the left and right edges will be set to 0. The finite difference method is not stable for values of average speed that cause diffusivity (D) to become large. For smaller speed values and smaller D values, the rate of mixing is not sufficient to keep the learner's attention. Therefore it was decided to show the FD solution with a greatly reduced value of D. This has the same features as the shape of the density but is much slower.

Alternative Method of Solving Equation for the Time Evolution of Density

An alternative solution to equation 3 is the error function solution:

`n(x,t)=(n_0)/2[erfc(x/(2sqrt(Dt)))] (4)`

where `n_0` is the starting density of particles on each side of the container and `erfc(u)` is the complement of the error function and is defined as:

`erfc(u)=1-2/sqrt(pi)int_0^u exp(-v^2)dv`

I will now differentiate both sides of equation 4 to show that it is a solution of equation 3:

`(del n)/(del t)=D(n_0)/4 x/t [exp(-(x^2)/(4 sqrt(D t)))]/sqrt(pi D t)`

`(del^2 n)/(del t^2)=(n_0)/4 x/t [exp(-(x^2)/(4 sqrt(D t)))]/sqrt(pi D t)`

where I have taken into account that the derivative of the integral with respect to its upper limit is just the integrand

`( d (erfc(u)))/(du)=-2/sqrt(pi) exp(-u^2)=-2/sqrt(pi) exp(-(x^2)/(4Dt))`

Equation 4 is the one that is plotted as the theory for the diffusion of particles.
As always, the learner is encouraged to press F12 (Windows) and select "Sources" to examine the HTML and javascript code that was used to create this page.