mb_dist

mdtools.statistics.mb_dist(v, temp=None, mass=None, var=None, drift=0, d=3)

Maxwell-Boltzmann speed distribution in \(d\)-dimensional space.[1]

\[p(v) = S_d \cdot (v - u)^{d-1} \left(\frac{1}{2\pi\sigma^2}\right)^{\frac{d}{2}} \exp{\left[-\frac{(v - u)^2}{2\sigma^2} \right]}\]

With the speed

\[v = \sqrt{\sum_{i=1}^d v_i^2},\]

a drift speed \(u\), the variance of the gaussian part (note that this is not the variance of \(p(v)\))

\[\sigma^2 = \frac{k_B T}{m}\]

and the surface area of a \(d\)-dimensional unit sphere[2]

\[\begin{split}S_d = \frac{2\pi^{\frac{d}{2}}}{\Gamma\left(\frac{d}{2}\right)} = \begin{cases} 2 , & d = 1\\ 2\pi , & d = 2\\ 4\pi , & d = 3\\ 2\pi^2 , & d = 4\\ \vdots & \end{cases}\end{split}\]

where \(\Gamma\) is the gamma function.

Parameters:

• v (array_like) – Array of speed values for which to evaluate \(p(v)\).

• temp (scalar or None, optional) – Temperature of the system. Must be provided if var is None.

• mass (scalar or None, optional) – Mass of the particles. Must be provided if var is None.

• var (scalar or None, optional) – Variance \(\sigma^2\) of the Maxwell-Boltzmann distribution. If None, it is calculated as \(\sigma^2 = k_B T / m\), where \(k_B\) is the Boltzmann constant, \(T\) is the temperature and \(m\) is the mass of the particles. If both, temp and mass, are provided, var is ignored.

• drift (scalar, optional) – Drift speed of the particles.

• d (int, optional) – Number of spatial dimensions in which the particles can move.

Returns:

p (scalar or array_like) – The function values \(p(v)\) at the given \(v\) values.

See also

mdtools.statistics.ekin_dist()

Distribution of kinetic energies in \(d\)-dimensional space

Notes

Be aware of the units of your input values! The unit of the Boltzmann constant \(k_B\) is J/K. Thus, temperatures must be given in K, masses in kg and velocities in m/s. If you parse var directly instead of temp and mass, v**2 and var must have the same units.

Note that all arguments should be positive to get a physically meaningful output. However, it is also possible to parse negative values to all input arguments.

Note that a (positive) drift velocity only leads to a shift of the distribution along the x axis (to higher values). Therefore, parts of the distribution that, without drift, would lie in the region of (unphysical) negative speeds now appear at positive speed values. The user is responsible to take care of this unphysical artifact. Generally, it’s easiest to calculate the Maxwell-Boltzmann distribution without drift at positive speed values and afterwards shift the outcome by the drift speed.

References

[1]

Wikipedia Maxwell-Boltzmann distribution § In n-dimensional space

[2]

Wikipedia n-sphere § Volume and area

Examples

>>> x = np.linspace(0, 4, 5)
>>> mdt.stats.mb_dist(x, var=1)
array([0.        , 0.48394145, 0.43192773, 0.07977327, 0.00428257])
>>> mdt.stats.mb_dist(x, var=2)
array([0.        , 0.21969564, 0.4151075 , 0.26759315, 0.08266794])
>>> mdt.stats.mb_dist(x, var=1, d=1)
array([7.97884561e-01, 4.83941449e-01, 1.07981933e-01, 8.86369682e-03,
       2.67660452e-04])
>>> mdt.stats.mb_dist(x, var=1, d=2)
array([0.        , 0.60653066, 0.27067057, 0.03332699, 0.00134185])
>>> mdt.stats.mb_dist(x, var=1, d=4)
array([0.        , 0.30326533, 0.54134113, 0.14997145, 0.0107348 ])
>>> mdt.stats.mb_dist(x, var=1, drift=1)
array([0.48394145, 0.        , 0.48394145, 0.43192773, 0.07977327])
>>> np.allclose(
...     mdt.stats.mb_dist(x, var=1)[:-1],
...     mdt.stats.mb_dist(x, var=1, drift=1)[1:],
...     rtol=0,
... )
True

>>> from scipy import constants
>>> temp, mass = 273.15, 1.660e-27
>>> var = constants.k * temp / mass
>>> np.allclose(
...     mdt.stats.mb_dist(x, temp=temp, mass=mass),
...     mdt.stats.mb_dist(x, var=var),
...     rtol=0,
... )
True
>>> np.allclose(
...     mdt.stats.mb_dist(x, temp=temp, mass=mass, d=2),
...     mdt.stats.mb_dist(x, var=var, d=2),
...     rtol=0,
... )
True
>>> np.allclose(
...     mdt.stats.mb_dist(x, temp=temp, mass=mass, drift=2),
...     mdt.stats.mb_dist(x, var=var, drift=2),
...     rtol=0,
... )
True

The area below the Maxwell-Boltzmann distribution is one:

>>> x = np.linspace(0, 20, 200)
>>> for d in range(1, 6):
...     np.isclose(
...         np.trapz(mdt.stats.mb_dist(x, var=1, d=d), x),
...         1,
...         rtol=0,
...         atol=1e-3,
...     )
True
True
True
True
True