# Can the Debye functions be implemented using built-in functions?

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.

• 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. Oct 3 at 19:29
• 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. Oct 3 at 20:03
• 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. Oct 3 at 20:11
• @Domen I edited the post. Oct 3 at 20:27

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}, d[x], x], {n, 5}];
Plot[soln1, {x, -6, 6},
PlotLegends -> Table[Subsuperscript[D, n, 1], {n, 5}]] So we can write:

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


or alternatively using the integral definition:

debyeD1[n_,
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 == n! Zeta[n + 1]}, d[x], x], {n, 5}];
Plot[soln2, {x, -6, 6},
PlotLegends -> Table[Subsuperscript[D, n, 2], {n, 5}]] 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 == 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 == n! Zeta[n + 1]}/.n->-1
(*{Derivative[d][x] == -(1/((-1 + E^x) x)), d == 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)$$.

• 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]. Oct 3 at 20:08
• 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.
– ydd
Oct 3 at 20:17
• 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. Oct 3 at 20:23
• @WillG Ah I see
– ydd
Oct 3 at 20:37
• There is a publication over approximations of the Debye function at researchgate.net/publication/… Oct 6 at 6:34