# Planck distribution integration

Does anyone know how to compute analytically or numerically the following integral (for $$T=10^4$$K)?:

$$n_\gamma=\frac{1}{\hbar^3\pi^2c^3}\int\limits_{2.1789\cdot 10^{-18}}^{+\infty}\dfrac{E^2\mathrm{d}E}{e^{\frac{E}{kT}}-1}$$

I tried with R, MATLAB, Maxima, Maple and Wolfram but I failed. I also search an analytical solution during a least a whole week....

Thx in advance for your help.

• What are the units of $2.1789\cdot 10^{-18}$? Nov 11, 2020 at 20:47
• What does “I failed” mean? What prevented you from obtaining a numerical integration? Nov 11, 2020 at 20:48
• Why is the lower limit not zero? Photons of arbitrarily low energy contribute to the number density. Nov 11, 2020 at 20:54
• @G.Smith The units are [J] ( the equivalent of 13.6 [eV]) Nov 11, 2020 at 21:47
• @G.Smith "i failed" mean that with all my softwares the convergence seems to fail. Nov 11, 2020 at 21:47

## 2 Answers

Sometimes computers have issues dealing with precision. Also try a change of variables like u=E/kT. Then factor out all the constants. So you get something like u^2du/(e^u-1) in your integral. This will simplify things in your numerical calculation. Once you calculate this integral, then you can just multiply back in the constants at the end.

• This 2.1789⋅10−18 is essentially zero, so you can make that assumption. You cannot assume the lower limit is zero. If you change variables to $x$, you will find out how very different from zero the lower limit is. Nov 12, 2020 at 4:22
• True, although I feel that won't be the biggest issue. Rescaling some of the variables may help though.
– Ali
Nov 12, 2020 at 4:34
• I have done the calculation, both by rescaling and not. If you set the lower limit to zero you will be wrong by about five orders of magnitude. Your statement “you can make that assumption” is false. Nov 12, 2020 at 4:42
• I just checked the lower limit is of similar order to kT, so you are correct we cannot assume it is zero. I will edit my response accordingly. Thanks!
– Ali
Nov 12, 2020 at 4:47
• Thanks for correcting. I’ve upvoted. Nov 12, 2020 at 4:49

Sometimes, one just need to do a little work before the analytic software can handle the integration. Here is a suggestion. First write it in the form: $$\frac{A}{\exp(Bx)-1}=\frac{A\exp(-Bx)}{1-\exp(-Bx)} .$$ Then expand the denominator $$\frac{A\exp(-Bx)}{1-\exp(-Bx)} = A\exp(-Bx)\sum_n \exp(-nBx) .$$ Now you (or Maple, etc.) should be able to integrate every term after which you can resum the series.

• Although one can expand the integrand in a series, it isn’t necessary to do that. Nov 12, 2020 at 4:45