def debye_integral(y):
        """
        Debye integral. Eq(5) in  doi.org/10.1016/j.comphy.2003.12.001

        Args:
            y (float): debye temperature/T, upper limit

        Returns:
            float: unitless
        """
        # floating point limit is reached around y=155, so values beyond that
        # are set to the limiting value(T-->0, y --> \infty) of
        # 6.4939394 (from wolfram alpha).
        factor = 3. / y ** 3
        if y < 155:
            integral = quadrature(lambda x: x ** 3 / (np.exp(x) - 1.), 0, y)
            return list(integral)[0] * factor
        else:
            return 6.493939 * factor