  In this paper, we show that efficient separated sum-of-exponentials
approximations can be constructed for the heat kernel in any dimension. In one
space dimension, the heat kernel admits an approximation involving a number of
terms that is of the order $O(\log(\frac{T}{\delta})
(\log(\frac{1}{\epsilon})+\log\log(\frac{T}{\delta})))$ for any $x\in\bbR$ and
$\delta \leq t \leq T$, where $\epsilon$ is the desired precision. In all
higher dimensions, the corresponding heat kernel admits an approximation
involving only $O(\log^2(\frac{T}{\delta}))$ terms for fixed accuracy
$\epsilon$. These approximations can be used to accelerate integral
equation-based methods for boundary value problems governed by the heat
equation in complex geometry. The resulting algorithms are nearly optimal. For
$N_S$ points in the spatial discretization and $N_T$ time steps, the cost is
$O(N_S N_T \log^2 \frac{T}{\delta})$ in terms of both memory and CPU time for
fixed accuracy $\epsilon$. The algorithms can be parallelized in a
straightforward manner. Several numerical examples are presented to illustrate
the accuracy and stability of these approximations.
