Module debye 

Debye model of thermal vibration: RMS displacement and its sampling.

Model

A monatomic harmonic crystal with the Debye density of states, cut off at the Debye frequency omega_D = k_B Theta_D / hbar, has a mean-square displacement of an atom along one Cartesian axis

u1^2 = (3 hbar^2 / (m k_B Theta_D)) * [ (T/Theta_D) D1(Theta_D/T) + 1/4 ]
D1(x) = (1/x) * integral_0^x  t / (exp(t) - 1)  dt

where m is the atomic mass and D1 is the Debye function of order 1 (called Phi in the channeling literature). The 1/4 is the zero-point term; the D1 term is the thermal part.

Source. The expression is the Debye-model evaluation of the harmonic mean-square displacement, 2W(q) = q^2 <u_x^2> with <u_x^2> = (hbar / (6 m)) * int g(w) coth(hbar w / 2 k_B T) / w dw (w is angular frequency; g(w) is the per-atom density of states normalised to three modes, int_0^{w_D} g dw = 3), as in N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston, 1976), Appendix L. We could not open the printed book; the equation was read in the form given in the “Mean squared displacement” section of the Wikipedia article “Debye function” (https://en.wikipedia.org/wiki/Debye_function, retrieved 2026-10-07), which cites that appendix and writes 2W(q) = (3/2) (hbar^2 q^2 / (m hbar omega_D)) [2 (k_B T/hbar omega_D) D1(hbar omega_D / k_B T) + 1/2]; dividing by q^2 gives the form above. The same quantity in the channeling literature is D. S. Gemmell, Rev. Mod. Phys. 46, 129 (1974), doi:10.1103/RevModPhys.46.129; we could not open that paper either, so its equation number is not cited here and its text was not consulted. As an independent check of the transcription, the unit tests integrate the coth form directly with the Debye density of states g(w) = 9 w^2 / w_D^3 (so int_0^{w_D} g dw = 3) and compare it with the closed form.

Limits. T -> 0: u1^2 -> 3 hbar^2 / (4 m k_B Theta_D) (zero-point). T >> Theta_D: D1(x) -> 1 - x/4 + x^2/36 - ..., so u1^2 -> 3 hbar^2 T / (m k_B Theta_D^2) = 3 k_B T / (m omega_D^2) (the classical equipartition value). The -x/4 term of D1 cancels the zero-point 1/4 exactly, so the relative correction is +x^2/36 with x = Theta_D/T, i.e. about 0.03 % at T = 10 Theta_D.

One-dimensional versus three-dimensional amplitude

ThermalVibration::rms_1d is u1 = sqrt(u1^2), the RMS displacement along one Cartesian axis (the quantity named u_1 in the channeling literature, which sets the Gaussian width of the transverse displacement distribution). ThermalVibration::rms_3d is the RMS length of the full displacement vector, u = sqrt(3) u1 for an isotropic (cubic) crystal. Sampling draws each Cartesian component independently from a Gaussian of standard deviation u1.

Debye temperature

A Debye temperature is not unique: it depends on the property it was fitted to (low-temperature specific heat, elastic constants, or the temperature dependence of X-ray or neutron Bragg intensities), and the values differ by tens of kelvin for the same crystal. The constants THETA_D_SI and friends are the room-temperature “Debye temperature” entries of the Ioffe Institute New Semiconductor Materials handbook; that source does not state the method for Si, Ge and GaAs, and for 3C-SiC it gives Goldberg et al. (2001) as its reference (also without a method). We therefore do not claim a specific-heat or X-ray determination for them. They are the only openable, citable values we found, and they are defaults only: the model takes Theta_D as an argument so a value matched to a given measurement can be used instead. A lattice-vibration amplitude fitted from X-ray Debye-Waller factors would be the more direct choice for channeling work; no openable source for such values was found (see docs/data-provenance.md, “Debye temperatures”).

For a compound (GaAs, SiC) one Theta_D is applied to every species with that species’ own mass, which is the usual single-parameter approximation and not a property of the Debye model itself.

Sampling and determinism

ThermalVibration::sample_displacement draws from the caller’s per-particle stream (crate::rng::ParticleRng) and always consumes exactly six u64 words (three Box-Muller pairs, of which only the cosine member is used), so a history’s stream position does not depend on the sampled values and results do not depend on the thread count. Box-Muller: G. E. P. Box and M. E. Muller, Ann. Math. Statist. 29, 610 (1958).

Structs

ThermalVibration

Debye-model vibration of one atomic species at one temperature.

Enums

DebyeError

Errors from building a ThermalVibration.

Constants

THETA_D_3C_SIC

Debye temperature of 3C-SiC in kelvin: 1200 K, “Debye temperature” row of https://www.ioffe.ru/SVA/NSM/Semicond/SiC/basic.html (retrieved 2026-10-07), which credits Goldberg et al. (2001). Method not stated.

THETA_D_GAAS

Debye temperature of GaAs in kelvin: 360 K, “Debye temperature” row of https://www.ioffe.ru/SVA/NSM/Semicond/GaAs/basic.html (retrieved 2026-10-07). Method not stated by the source.

THETA_D_GE

Debye temperature of Ge in kelvin: 374 K, “Debye temperature” row of https://www.ioffe.ru/SVA/NSM/Semicond/Ge/basic.html (retrieved 2026-10-07). Method not stated by the source.

THETA_D_SI

Debye temperature of Si in kelvin: 640 K, “Debye temperature” row of https://www.ioffe.ru/SVA/NSM/Semicond/Si/basic.html (retrieved 2026-10-07). Determination method not stated by the source; see the module docs.

Functions

debye_d1

Debye function of order 1, D1(x) = (1/x) int_0^x t/(e^t - 1) dt, for x > 0. D1(0) = 1, D1(x) -> pi^2/(6x) as x -> infinity.

sample_gaussian_3d

Three independent Gaussian components of standard deviation sigma (m), by Box-Muller (module docs). Consumes exactly six u64 words from rng whatever sigma is; sigma = 0 gives components equal to zero (possibly -0.0).