Boundary integral formulations of the heat equation involve time convolutions in addition to surface potentials. If M is the number of time steps and N is the number of degrees of freedom of the spatial discretization then the direct computation of a heat potential involves order N²M² operations. This article describes a fast method to compute three-dimensional heat potentials which is based on Chebyshev interpolation of the heat kernel in both space and time. The computational complexity is order p⁴q²NM operations, where p and q are the orders of the polynomial approximation in space and time.