It is claimed in the comments here that the Debye functions can be implemented using built-in special functions. This is clearly true for some Debye functions, e.g., $D_n^{(1)}(x)$ for $n = 1, 2, 3$ (see linked article, which gives these explicitly). However, I'm not sure whether the comment holds in general.

For example, I'm interested in the following integral:

$$D_{-1}^{(2)}(x) = \int_x^\infty \frac{dt}t \frac1{e^t - 1}.$$

Can this integral be represented exactly, in closed form, using built-in special functions? More generally, can the functions $D_n^{(1)}$ and $D_n^{(2)}$ be represented exactly for all integer $n$?

Edit: By "represented exactly, in closed form, using built-in special functions" I mean using built-in functions that are quickly evaluated numerically, without relying on NIntegrate or NSolve or other explicit numerical functions.

  • $\begingroup$ I am not entirely sure what you mean and what your end-goal is. For example, $D_{-1}^{(2)}(x)$ is respresented in Mathematica exactly with Integrate[1/(t (Exp[t] - 1)), {t, x, ∞}]. Perhaps you may want to elaborate what you want to do with these functions. $\endgroup$
    – Domen
    Oct 3, 2023 at 19:29
  • $\begingroup$ I mean represented exactly using built-in special functions that Mathematica can numerically evaluate "quickly" and to arbitrary precision. Compare, e.g., Timing@Integrate[x^2 E^-x, {x, 1., Infinity}] to Timing@Gamma[3, 1.]. In this example, Gamma is the "built-in special function" that computes the given integral. $\endgroup$
    – WillG
    Oct 3, 2023 at 20:03
  • $\begingroup$ Well, then perhaps you should make it more clear in your question that you are not really looking for an "exact" symbolic representation but for an efficient numerical evaluation. $\endgroup$
    – Domen
    Oct 3, 2023 at 20:11
  • $\begingroup$ @Domen I edited the post. $\endgroup$
    – WillG
    Oct 3, 2023 at 20:27

1 Answer 1


This only partially answers your question so I'll delete if it is not what you're looking for.

The first type of functions, $D^{(1)}_n (x)$ satisfy the differential equation $$D^{'(1)}_n (x) = \frac{x^n}{e^x-1} \\ \mathrm{with} \\ D^{(1)}_n (0) = 0$$

(I got the initial condition from looking at the integral definition on the MathWorld site at x = 0)

Edit: I missed on the MathWorld site it also mentions that $D^{(1)}_n (x)$ is only valid for $|x| < 2 \pi$, so just keep that in mind.

diff[n_] = D[d[x], x] == x^n/(E^x - 1);
soln1 = Table[ DSolveValue[{diff[n], d[0] == 0}, d[x], x], {n, 5}];
Plot[soln1, {x, -6, 6}, 
 PlotLegends -> Table[Subsuperscript[D, n, 1], {n, 5}]]

Mathematica graphics

So we can write:

debyeD1[n_, x_] := DSolveValue[{diff[n], d[0] == 0}, d[x], x]

or alternatively using the integral definition:

  x_] := (Integrate[x^n/(Exp[x] - 1), x] + (n! + Zeta[n + 1]))

The second type of functions, $D^{(2)}_n (x)$ satisfy the differential equation $$D^{'(2)}_n (x) = -\frac{x^n}{e^x-1} \\ \mathrm{with} \\ D^{(2)}_n (0) = n! \zeta (n+1)$$

This initial condition is also from the MathWorld site, as it says: $$D^{(1)}_n (x) + D^{(2)}_n (x) = n! \zeta (n+1)$$ And we know $D^{(1)}_n (0)$ is 0

diff2[n_] = D[d[x], x] == -(x^n/(E^x - 1));
soln2 = Table[ 
   DSolveValue[{diff2[n], d[0] == n! Zeta[n + 1]}, d[x], x], {n, 5}];
Plot[soln2, {x, -6, 6}, 
 PlotLegends -> Table[Subsuperscript[D, n, 2], {n, 5}]]

Mathematica graphics

And confirming that the condition $$D^{(1)}_n (x) + D^{(2)}_n (x) = n! \zeta (n+1)$$ is met:

sums = soln1 + soln2;
sums/Table[n! Zeta[n + 1], {n, 5}] // Simplify
(*{1, 1, 1, 1, 1}*)

So we can write:

debyeD2[n_, x_] := 
 DSolveValue[{diff2[n], d[0] == n! Zeta[n + 1]}, d[x], x]

or alternatively using the integral definition:

debyeD2[n_, x_] := -Integrate[x^n/(Exp[x] - 1), x]

The only problem is I don't know how to do $D^{(2)}_{-1} (x)$

neg1 = {diff2[n], d[0] == n! Zeta[n + 1]}/.n->-1
(*{Derivative[1][d][x] == -(1/((-1 + E^x) x)), d[0] == ComplexInfinity}*)

because n! is not defined I believe at negative integers. Maybe someone else can fill in on this part. The integral definition doesn't evaluate either so I can't figure out how to show an exact form for $D^{(2)}_{-1} (x)$.

  • $\begingroup$ This is very nice, but my intention was to avoid calling on numerical computation functions. E.g., in the way that Gamma[3, 1] represents (and quickly computes) Integrate[x^2 E^(-x), {x, 1, Infinity]. $\endgroup$
    – WillG
    Oct 3, 2023 at 20:08
  • $\begingroup$ These are not numerical computation functions, they are all exact solutions. But I think you are saying you want to find a general solution for Integrate[x^n/(-1 + E^x), x, Assumptions -> n \[Element] Integers] ? This seems challenging, but perhaps it is possible. $\endgroup$
    – ydd
    Oct 3, 2023 at 20:17
  • $\begingroup$ Yes, but if Mathematica had a built-in function Debye[n, z] that calculated this automatically, I would consider that a solution. I guess one could argue that there's nothing special about such built-in functions because they still need to execute algorithms to compute values, but they do tend to be very fast and easy to use. $\endgroup$
    – WillG
    Oct 3, 2023 at 20:23
  • $\begingroup$ @WillG Ah I see $\endgroup$
    – ydd
    Oct 3, 2023 at 20:37
  • $\begingroup$ There is a publication over approximations of the Debye function at researchgate.net/publication/… $\endgroup$
    – Roland F
    Oct 6, 2023 at 6:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.