Summary: We present boundary integral representations of several initial boundary value problems related to the heat equation. A Galerkin discretization with piecewise constant functions in time and piecewise linear functions in space leads to optimal a priori error estimates, provided that the meshwidths in space and time satisfy \(h_t=\mathcal{O}(h_x^2)\). Each time step involves the solution of a linear system, whose spectral condition number is independent of the refinement under the same assumption on the mesh. We show that if the parabolic multipole method is used to apply parabolic boundary integral operators, the overall complexity of the scheme is log-linear while preserving the convergence of the Galerkin discretization method. The theoretical estimates are confirmed numerically at the end of the paper.