Approximate solution of the heat equation

Question

I'm reading the article "Singularity formation for the two dimensional harmonic map flow into $S^2$" from J.Davila, M.del Pino and J.Wei and at some point, there are some computations I don't totally get (p.365). We have the following equation $$u_t - \Delta u - \frac{2r}{r^2 + \lambda^2} \dot p(t) = 0, \quad \text{on }\mathbb R^2 \times [-T, +\infty)$$ for $r = |x|$ in $\mathbb R^2$, $r \gg \lambda$ and $p$ some function depending on $t$ (which is unknown but satisfies all the good properties needed). Now, the authors say that we have a good approximate solution of this problem by considering the function $$u(r, t) = -\int_{-T}^t 2\dot p(s) r k(z(r), t - s)ds,$$ where $$z(r) = \sqrt{r^2 + \lambda^2}, \quad \text{and} \quad k(z, t) = \frac{1 - \exp(-z^2/4t)}{z^2}.$$ I do not totally get how do they find such a function. It looks like the Duhamel formula for the heat equation, but I do not understand how we can get this $1 - \exp(-z^2/4t)$ because the heat kernel is given by $\exp(-z^2/4t)/4t$. I thought maybe expand the heat kernel to the first order then use polar coordinate, but I couldn't find anything near the above expression. Any idea ?