"""
Debye functions.
Debye functions appear in solid-state physics for computing
thermodynamic properties of solids (heat capacity, entropy).
Performance
-----------
This module uses Numba JIT compilation with rapidly convergent series
expansions (Abramowitz & Stegun 27.1.1-27.1.3), providing ~10-50x speedup
for batch computations compared to scipy.integrate.quad.
Accuracy
--------
Typical relative error is ~1e-16 (machine precision), measured against a
30-40-digit mpmath oracle for n in {1, 4, 6, 8, 9, 10}, x in
{0.5, 0.99, 1.0, 1.01, 1.9, 1.99, 1.999, 2.0, 10.0}. The worst points on
that grid are ~2e-11 just below the x0=2.0 boundary (n=9/10 near
x=1.99-1.999, where the small-x series is nearing the edge of its
useful precision) and ~4e-12 exactly at the boundary (n=9/10 at x=2.0,
where the large-x branch's own alternating-sum cancellation happens to
be smallest). The small-x/large-x branch switch is at x0=2.0: before
this was tightened from x0=1.0, the same grid measured up to 8.6e-9 at
x=1.0-1.01 for n>=8 (catastrophic cancellation in the large-x branch
right at the old boundary).
"""
from typing import Any
import numpy as np
from numba import njit, prange
from numpy.typing import ArrayLike, NDArray
from scipy.special import zeta
# Pre-compute zeta values for common orders (n=1 to 10)
_ZETA_VALUES = np.array([zeta(k + 1) for k in range(11)])
# B_{2k} / (2k)! for k = 1..10 (Bernoulli numbers over factorials), used in
# the small-x expansion of t/(e^t - 1) = 1 - t/2 + sum B_{2k} t^{2k}/(2k)!
_BERNOULLI_COEF = np.array(
[
1.0 / 12.0,
-1.0 / 720.0,
1.0 / 30240.0,
-1.0 / 1209600.0,
1.0 / 47900160.0,
-691.0 / 1307674368000.0,
1.0 / 74724249600.0,
-3617.0 / 10670622842880000.0,
43867.0 / 5109094217170944000.0,
-174611.0 / 802857662698291200000.0,
]
)
@njit(cache=True)
def _debye_small_x(x: float, n: int, coef: np.ndarray[Any, Any]) -> float:
"""
Bernoulli series expansion for x < 1 (converges for |x| < 2*pi).
D_n(x) = 1 - n*x/(2*(n+1)) + n * sum_k B_{2k}/(2k)! * x^{2k}/(2k+n)
"""
result = 1.0 - n * x / (2.0 * (n + 1))
x2 = x * x
xp = 1.0
for k in range(len(coef)):
xp *= x2
result += n * coef[k] * xp / (2 * (k + 1) + n)
return result
@njit(cache=True)
def _debye_large_x(x: float, n: int, n_fact: float, zeta_n_plus_1: float) -> float:
"""
Complement series for x >= 1 (A&S 27.1.2-27.1.3).
D_n(x) = (n/x^n) * [n! * zeta(n+1)
- sum_{j>=1} e^{-jx} * (n!/j^{n+1}) * sum_{i=0}^{n} (jx)^i/i!]
"""
total = n_fact * zeta_n_plus_1
for j in range(1, 500):
jx = j * x
if jx > 700.0:
break
# Partial exponential sum: sum_{i=0}^{n} (jx)^i / i!
s = 1.0
term = 1.0
for i in range(1, n + 1):
term *= jx / i
s += term
contrib = np.exp(-jx) * n_fact / float(j) ** (n + 1) * s
total -= contrib
if contrib < 1e-17 * total:
break
return n / x**n * total
@njit(cache=True, parallel=True)
def _debye_batch(
n: int,
x_arr: np.ndarray[Any, Any],
zeta_n_plus_1: float,
coef: np.ndarray[Any, Any],
) -> np.ndarray[Any, Any]:
"""
Batch computation of Debye function for array input.
Parameters
----------
n : int
Order of the Debye function.
x_arr : ndarray
Array of x values.
zeta_n_plus_1 : float
Pre-computed zeta(n+1) value.
coef : ndarray
Bernoulli coefficients B_{2k}/(2k)! for the small-x series.
Returns
-------
ndarray
Debye function values.
"""
result = np.empty(len(x_arr), dtype=np.float64)
n_fact = 1.0
for k in range(1, n + 1):
n_fact *= k
for i in prange(len(x_arr)): # ty: ignore[not-iterable]
xi = x_arr[i]
if xi == 0.0:
result[i] = 1.0
elif xi < 2.0:
result[i] = _debye_small_x(xi, n, coef)
else:
result[i] = _debye_large_x(xi, n, n_fact, zeta_n_plus_1)
return result
[docs]
def debye(
n: int,
x: ArrayLike,
) -> NDArray[np.floating]:
"""
Debye function D_n(x).
The Debye function of order n is defined as:
D_n(x) = (n/x^n) * integral from 0 to x of t^n / (exp(t) - 1) dt
Parameters
----------
n : int
Order of the Debye function (positive integer).
x : array_like
Argument of the function, x >= 0.
Returns
-------
D : ndarray
Values of D_n(x).
Notes
-----
Special cases:
- D_n(0) = 1
- D_n(inf) = n! * zeta(n+1) / x^n -> 0
The Debye function D_3(x) appears in the heat capacity
of solids at low temperatures.
This implementation uses Numba JIT compilation for performance,
achieving ~10-50x speedup compared to scipy.integrate.quad for
batch computations.
Examples
--------
>>> float(debye(3, 0)[0]) # D_3(0) = 1
1.0
>>> round(float(debye(3, 1)[0]), 6)
0.674416
>>> round(float(debye(3, 10)[0]), 6)
0.019296
References
----------
- Debye, P. (1912). "Zur Theorie der spezifischen Wärmen".
Annalen der Physik, 344(14), 789-839.
"""
if n < 1:
raise ValueError(f"Order n must be >= 1, got {n}")
x = np.atleast_1d(np.asarray(x, dtype=np.float64))
# Get pre-computed zeta value if available, otherwise compute
if n < len(_ZETA_VALUES):
zeta_n_plus_1 = _ZETA_VALUES[n]
else:
zeta_n_plus_1 = zeta(n + 1)
return _debye_batch(n, x, zeta_n_plus_1, _BERNOULLI_COEF)
[docs]
def debye_1(x: ArrayLike) -> NDArray[np.floating]:
"""
First-order Debye function D_1(x).
Parameters
----------
x : array_like
Argument of the function, x >= 0.
Returns
-------
D : ndarray
Values of D_1(x).
Notes
-----
D_1(x) = (1/x) * integral from 0 to x of t / (exp(t) - 1) dt
"""
return debye(1, x)
[docs]
def debye_2(x: ArrayLike) -> NDArray[np.floating]:
"""
Second-order Debye function D_2(x).
Parameters
----------
x : array_like
Argument of the function, x >= 0.
Returns
-------
D : ndarray
Values of D_2(x).
Notes
-----
D_2(x) = (2/x^2) * integral from 0 to x of t^2 / (exp(t) - 1) dt
"""
return debye(2, x)
[docs]
def debye_3(x: ArrayLike) -> NDArray[np.floating]:
"""
Third-order Debye function D_3(x).
This is the most commonly used Debye function, appearing in
the heat capacity of solids.
Parameters
----------
x : array_like
Argument of the function, x >= 0.
Returns
-------
D : ndarray
Values of D_3(x).
Notes
-----
D_3(x) = (3/x^3) * integral from 0 to x of t^3 / (exp(t) - 1) dt
The heat capacity of a solid in the Debye model is:
C_V = 9 * N * k_B * (T/Θ_D)^3 * D_3(Θ_D/T)
where Θ_D is the Debye temperature.
"""
return debye(3, x)
[docs]
def debye_4(x: ArrayLike) -> NDArray[np.floating]:
"""
Fourth-order Debye function D_4(x).
Parameters
----------
x : array_like
Argument of the function, x >= 0.
Returns
-------
D : ndarray
Values of D_4(x).
Notes
-----
D_4(x) = (4/x^4) * integral from 0 to x of t^4 / (exp(t) - 1) dt
This appears in computing the entropy of solids.
"""
return debye(4, x)
[docs]
def debye_heat_capacity(
temperature: ArrayLike,
debye_temperature: float,
) -> NDArray[np.floating]:
"""
Debye model heat capacity (normalized).
Computes C_V / (3*N*k_B) using the Debye model.
Parameters
----------
temperature : array_like
Temperature in Kelvin.
debye_temperature : float
Debye temperature Θ_D in Kelvin.
Returns
-------
cv_normalized : ndarray
Normalized heat capacity C_V / (3*N*k_B).
Multiply by 3*N*k_B for actual heat capacity.
Notes
-----
The Debye model heat capacity is:
C_V / (3*N*k_B) = 4*D_3(x) - 3*x/(e^x - 1), with x = Θ_D/T
Limits:
- High T (T >> Θ_D): C_V -> 3*N*k_B (classical)
- Low T (T << Θ_D): C_V ~ (4*π^4/5) * (T/Θ_D)^3 (quantum)
Examples
--------
>>> # Aluminum at room temperature (Θ_D ≈ 428 K)
>>> cv = debye_heat_capacity(300, 428) # ~0.91
"""
T = np.asarray(temperature, dtype=np.float64)
theta_D = float(debye_temperature)
if np.any(T <= 0):
raise ValueError("Temperature must be positive")
if theta_D <= 0:
raise ValueError("Debye temperature must be positive")
x = np.atleast_1d(theta_D / T)
# C_V / (3*N*k_B) = 4*D_3(x) - 3*x/(e^x - 1)
# (obtained by integrating the Debye phonon spectrum by parts)
with np.errstate(over="ignore"):
boltzmann_term = np.where(x > 500, 0.0, 3.0 * x / np.expm1(np.minimum(x, 700)))
return 4.0 * debye(3, x) - boltzmann_term
[docs]
def debye_entropy(
temperature: ArrayLike,
debye_temperature: float,
) -> NDArray[np.floating]:
"""
Debye model entropy (normalized).
Computes S / (3*N*k_B) using the Debye model.
Parameters
----------
temperature : array_like
Temperature in Kelvin.
debye_temperature : float
Debye temperature Θ_D in Kelvin.
Returns
-------
s_normalized : ndarray
Normalized entropy S / (3*N*k_B).
Notes
-----
The entropy in the Debye model is:
S / (3*N*k_B) = (4/3)*D_3(Θ_D/T) - ln(1 - exp(-Θ_D/T))
"""
T = np.asarray(temperature, dtype=np.float64)
theta_D = float(debye_temperature)
if np.any(T <= 0):
raise ValueError("Temperature must be positive")
if theta_D <= 0:
raise ValueError("Debye temperature must be positive")
x = np.atleast_1d(theta_D / T)
# log1p(-e^{-x}) is accurate for all x > 0, including large x where
# e^{-x} underflows harmlessly to 0
log_term = np.log1p(-np.exp(-np.minimum(x, 700)))
return (4.0 / 3.0) * debye(3, x) - log_term
__all__ = [
"debye",
"debye_1",
"debye_2",
"debye_3",
"debye_4",
"debye_heat_capacity",
"debye_entropy",
]