Consider the heat equation
$$\dfrac{\partial u}{\partial t}=\dfrac{\partial^2 u}{\partial x^2} +G(x,t).$$
with the condition $u(x,0)=f(x)$. When $G(x,t)=0$ it is quite easy to solve it using Fourier transform. Taking the Fourier transform in the $x$ variable we have
$$\dfrac{\partial \hat{u}}{\partial t}=-4\pi^2 \xi^2 \hat{u}(\xi,t)$$
With the obvious solution
$$\hat{u}(\xi,t)=A(\xi)e^{-4\pi^2\xi^2 t}.$$
If we then impose $u(x,0)=f(x)$ we get $\hat{u}(\xi,0)=\hat{f}(\xi)$ and hence
$$\hat{u}(\xi,t)=\hat{f}(\xi)e^{-4\pi^2\xi^2 t}$$
Which upon introducing the heat kernel is the same as $u(x,t)=f\ast \mathcal{H}_t(x)$. So this is quite straightforward.
Now, I'm trying to deal with the case $G\neq 0$, with $G\in \mathcal{S}(\mathbb{R}^2)$. I've found that the solution should be
$$u(x,t)=f\ast \mathcal{H}_t(x)+\int_0^t \int_{-\infty}^{\infty} G(y,s)\mathcal{H}_{t-s}(x-y)dyds.$$
Now I really have no idea where this comes from.
How does one arrive at this solution using Fourier transform? And how does one actually prove rigorously it is the desired solution?
Denote by $\hat G$ the Fourier transform of $G$ with respect to $x$. You write (I omit all factors coming from the scaling of FT): $$\partial_t\hat u(\xi, t) = - \xi^2 \hat u(\xi, t) + \hat G(\xi,t).$$ This is an inhomogenous linear differential equation of the first order, the solution writes by a well-known formula proven by integrating factors: $$\hat u(\xi, t) = \hat f(\xi) e^{-\xi^2t} + \int_0^te^{-\xi^2(t-s)} \hat G(\xi,s) ds$$ or $$\hat u(\xi, t) = \hat f(\xi) \hat H(\xi, t) + \int_0^t \hat H(\xi, t-s) \hat G(\xi,s) ds,$$ which then leads to: $$ u = f * H_t + \int_0^t H_{t-s}* G(s) ds.$$