Debye function

In mathematics, the family of Debye functions is defined by $$D_n(x) = \frac{n}{x^n} \int_0^x \frac{t^n}{e^t - 1}\,dt.$$

The functions are named in honor of Peter Debye, who came across this function (with n = 3) in 1912 when he analytically computed the heat capacity of what is now called the Debye model.

Relation to other functions
The Debye functions are closely related to the polylogarithm.

Series expansion
They have the series expansion $$D_n(x) = 1 - \frac{n}{2(n+1)} x + n \sum_{k=1}^\infty \frac{B_{2k}}{(2k+n)(2k)!} x^{2k}, \quad |x| < 2\pi,\ n \ge 1,$$ where $$B_n$$ is the $n$-th Bernoulli number.

Limiting values
$$\lim_{x \to 0} D_n(x) = 1.$$ If $$\Gamma$$ is the gamma function and $$\zeta$$ is the Riemann zeta function, then, for $$x \gg 0$$, $$D_n(x) = \frac{n}{x^n} \int_0^x \frac{t^n\,dt}{e^t-1} \sim \frac{n}{x^n}\Gamma(n + 1) \zeta(n + 1), \qquad \operatorname{Re} n > 0,$$

Derivative
The derivative obeys the relation $$x D^{\prime}_n(x) = n \left(B(x) - D_n(x)\right),$$ where $$B(x) = x/(e^x-1)$$ is the Bernoulli function.

The Debye model
The Debye model has a density of vibrational states $$g_\text{D}(\omega) = \frac{9\omega^2}{\omega_\text{D}^3} \,, \qquad 0\le\omega\le\omega_\text{D}$$ with the $ω_{D}$.

Internal energy and heat capacity
Inserting $g$ into the internal energy $$U = \int_0^\infty d\omega\,g(\omega)\,\hbar\omega\,n(\omega)$$ with the Bose–Einstein distribution $$n(\omega) = \frac{1}{\exp(\hbar\omega / k_\text{B} T)-1}.$$ one obtains $$U = 3 k_\text{B}T \, D_3(\hbar\omega_\text{D} / k_\text{B}T).$$ The heat capacity is the derivative thereof.

Mean squared displacement
The intensity of X-ray diffraction or neutron diffraction at wavenumber q is given by the Debye-Waller factor or the Lamb-Mössbauer factor. For isotropic systems it takes the form $$\exp(-2W(q)) = \exp\left(-q^2\langle u_x^2\rangle\right).$$ In this expression, the mean squared displacement refers to just once Cartesian component $u_{x}$ of the vector $u$ that describes the displacement of atoms from their equilibrium positions. Assuming harmonicity and developing into normal modes, one obtains $$2W(q) = \frac{\hbar^2 q^2}{6M k_\text{B}T} \int_0^\infty d\omega \frac{k_\text{B}T}{\hbar\omega}g(\omega) \coth\frac{\hbar\omega}{2k_\text{B}T}=\frac{\hbar^2 q^2}{6M k_\text{B}T} \int_0^\infty d\omega \frac{k_\text{B}T}{\hbar\omega} g(\omega) \left[\frac{2}{\exp(\hbar\omega/k_\text{B}T)-1}+1\right].$$ Inserting the density of states from the Debye model, one obtains $$2W(q) = \frac{3}{2} \frac{\hbar^2 q^2}{M\hbar\omega_\text{D}} \left[2\left(\frac{k_\text{B}T}{\hbar\omega_\text{D}}\right) D_1{\left(\frac{\hbar\omega_\text{D}}{k_\text{B}T}\right)} + \frac{1}{2}\right].$$ From the above power series expansion of $$D_1$$ follows that the mean square displacement at high temperatures is linear in temperature $$2W(q) = \frac{3 k_\text{B}T q^2}{M\omega_\text{D}^2}.$$ The absence of $$\hbar$$ indicates that this is a classical result. Because $$D_1(x)$$ goes to zero for $$x \to \infty$$ it follows that for $$T = 0$$ $$2W(q)=\frac{3}{4}\frac{\hbar^2 q^2}{M\hbar\omega_\text{D}}$$ (zero-point motion).

Implementations

 * Fortran 77 code
 * Fortran 90 version
 * C version of the GNU Scientific Library
 * Fortran 90 version
 * C version of the GNU Scientific Library
 * C version of the GNU Scientific Library
 * C version of the GNU Scientific Library