We present a method for the fast and accurate computation of distributed heat potentials in two dimensions. The distributed source is assumed to be given in terms of piecewise space-time Chebyshev polynomials. We discretize uniformly in time, whereas in space the polynomials are defined on the leaf nodes of a quadtree. The quadtree can vary at each time step. We combine a product integration rule with fast algorithms (fast heat potentials, nonuniform FFT, fast Gauss transform) to obtain a high-order accurate method with optimal complexity. If N is the number of time steps, M is the maximum number of leaf nodes over all the time steps and the input contains a qth-order polynomial representation of f, then, our method requires OΓ°q 3 MN log MΓ work to evaluate the heat potential at arbitrary MN space-time target locations. The overall convergence rate of the method is of order q. We present numerical experiments for q = 4, 8, and 16, and we verify the theoretical convergence rate of the method. When the solution is sufficiently smooth, the 16th-order variant results in significant computational savings, even in the case in which we require only a few digits of accuracy.