Based on your calculations, you have reached to the following stage
$$ U(\omega,t)= \frac{1}{2\sqrt\pi} e^{-\frac{bt^3}{3}} e^{-\frac{\omega^2}{4}(4kt+1)}. $$
Now, all you need to do is to find the inverse Fourier transform w.r.t. $\omega$ to get $u(x,t)$,
$$ u(x,t)= \int_{-\infty}^\infty U(\omega, t)e^{-i\omega x}d\omega= \frac{1}{2\sqrt\pi} e^{-\frac{bt^3}{3}}\int_{-\infty}^\infty e^{-\frac{\omega^2}{4}(4kt+1)}e^{-i\omega x}d\omega $$
$$\implies u(x,t) = \frac{1}{2\sqrt\pi} e^{-\frac{bt^3}{3}}\int_{-\infty}^\infty e^{-\frac{a\,\omega^2}{4}}e^{-i\omega x}d\omega,$$
where $ a = 4kt+1. $
Can you find the last integral now?
Consider the Cauchy problem for the heat equation:
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial u}{\partial t} - \Delta u = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
u(x,0) = u_{0} & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
where $u_{0} : \mathbb{R}^{d} \to \mathbb{R}$ is a "nice" function. (A convenient choice is $u_{0} \in \mathscr{S}(\mathbb{R}^{d})$.)
The beauty of the Fourier transform in this context is $\Delta$ becomes an algebraic "multiplier." Indeed, if I multiply all sides by $e^{-i 2 \pi \langle \xi, x \rangle}$ and integrate with respect to $x$, then I obtain
\begin{align*}
\int_{\mathbb{R}^{d}} \frac{\partial u}{\partial t}(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx &= \frac{\partial}{\partial t} \left(\int_{\mathbb{R}^{d}} u(x,t) e^{- i 2 \pi \langle \xi, x \rangle} \, dx \right) = \frac{\partial \hat{u}}{\partial t}(\xi,t) \\
\int_{\mathbb{R}^{d}} \frac{\partial^{2} u}{\partial x_{j}^{2}}(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx &= (- i2 \pi \xi_{j})^{2} \int_{\mathbb{R}^{d}} u(x,t) e^{-i 2 \pi \langle \xi, x \rangle} \, dx = - 4 \pi^{2} \xi_{j}^{2} \hat{u}(\xi,t),
\end{align*}
where in the second line I integrated by parts twice and assumed that $\lim_{|x| \to \infty} u(x,t) = 0$ uniformly (for each fixed $t$). (Note: Here $\langle \xi, x \rangle = \sum_{j = 1}^{d} \xi_{j} x_{j}$, that is, it's a fancy way of writing the dot product.) Thus, summing over $j$, we obtain
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial u}{\partial t}(\xi,t) + 4 \pi^{2} |\xi|^{2} \hat{u}(\xi,t) = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\hat{u}(\xi,0) = \hat{u}_{0}(\xi) & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
For each fixed $\xi \in \mathbb{R}^{d}$, this is an ODE. In particular,
$$\hat{u}(\xi,t) = \hat{u}_{0}(\xi) e^{- 4 \pi^{2} |\xi|^{2} t}.$$
Therefore, by Fourier inversion,
\begin{align*}
u(x,t) = \int_{\mathbb{R}^{d}} \hat{u}_{0}(\xi) e^{-4 \pi^{2} |\xi|^{2} t + i 2 \pi \langle \xi, x \rangle} \, d\xi
\end{align*}
This isn't entirely what we want: we would like to express the answer in terms of $u_{0}$, not $\hat{u}_{0}$.
Recall the following fact about the Fourier transform: if $f,g$ are "nice" functions (e.g. $f,g \in \mathscr{S}(\mathbb{R}^{d})$), then
$$(f * g)(x) = \int_{\mathbb{R}^{d}} \hat{f}(\xi) \hat{g}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d\xi.$$
In other words, "convolution in space corresponds to multiplication in frequency." Thus, let $G_{t}$ denote the function satisfying
$$\hat{G}_{t}(\xi) = e^{-4 \pi^{2} |\xi|^{2} t}.$$ By the last remark, we are led to guess that
$$\int_{\mathbb{R}^{d}} u_{0}(y) G_{t}(x - y) \, dy = \int_{\mathbb{R}^{d}} \hat{u}_{0}(\xi) \hat{G}_{t}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d\xi.$$
This leads naturally to the question: what is $G_{t}$? By Fourier inversion,
$$G_{t}(x) = \int_{\mathbb{R}^{d}} e^{- 4 \pi^{2} |\xi|^{2} t + i 2 \pi \langle \xi ,x \rangle} \, d\xi.$$
The RHS can be evaluated using complex analysis (or found in a Fourier transform table). The results is:
$$G_{t}(x) = (4 \pi t)^{-\frac{d}{2}} e^{- \frac{|x|^{2}}{4t}}.$$
Since $G_{t}$ is a very nice function (i.e. $G_{t} \in \mathscr{S}(\mathbb{R}^{d})$), our guess is justified and we obtain
$$u(x,t) = (4 \pi t)^{-\frac{d}{2}} \int_{\mathbb{R}^{d}} u_{0}(y) e^{- \frac{|y - x|^{2}}{4t}} \, dy.$$
In the case when $u_{0} = \delta_{0}$, we find $u(x,t) = G_{t}(x)$. Therefore, we showed $\tilde{G}(x,t) = G_{t}(x)$ satisfies
\begin{align*}
\left\{\begin{array}{r l}
\frac{\partial \tilde{G}}{\partial t} - \Delta \tilde{G} = 0 & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\tilde{G}(dx,0) = \delta_{0}(dx) & \text{on} \, \, \mathbb{R}^{d}.
\end{array} \right.
\end{align*}
The function $\tilde{G}$ is called the heat kernel.
A classical question is now: how can we solve inhomogeneous problem, i.e.
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial v}{\partial t} - \Delta v = f(x,t) & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
v(x,0) = u_{0} & \text{on} \, \, \mathbb{R}^{d} \end{array} \right.
\end{align*}
By linearity, it suffices to consider the case when $u_{0} = 0$. (Otherwise, add the solution with $0$ initial data to the solution of the homogeneous equation with initial data $u_{0}$.) If we take the Fourier transform in space again, we find
\begin{align*}
\left\{ \begin{array}{r l}
\frac{\partial \hat{v}}{\partial t}(\xi,t) + 4 \pi^{2} |\xi|^{2} \hat{v}(\xi,t) = \hat{f}(\xi,t) & \text{in} \, \, \mathbb{R}^{d} \times (0,\infty) \\
\hat{v}(\xi,0) = 0 & \text{on} \, \, \mathbb{R}^{d}.
\end{array} \right.
\end{align*}
Recall that this linear ODE can be solved using Duhamel's formula:
\begin{align*}
\hat{v}(\xi,t) = \int_{0}^{t} \hat{f}(\xi,s) e^{- 4 \pi^{2} |\xi|^{2} (t - s)} \, ds.
\end{align*}
If we apply the inverse Fourier transform and Fubini's Theorem, we find
\begin{align*}
v(x,t) &= \int_{0}^{t} \left( \int_{\mathbb{R}^{d}} \hat{f}(\xi,s) e^{- 4 \pi^{2} |\xi|^{2} (t - s) + i 2 \pi \langle \xi, x \rangle} \, d \xi \right) \, ds \\
&= \int_{0}^{t} \left( \int_{\mathbb{R}^{d}} \hat{f}(\xi,s) \hat{G}_{t - s}(\xi) e^{i 2 \pi \langle \xi, x \rangle} \, d \xi \right) \, ds \\
&= \int_{0}^{t} \left(\int_{\mathbb{R}^{d}} f(y,s) G_{t - s}(x - y) \, dy \right) \, ds.
\end{align*}
In the case when $f(dx,dt) = \delta_{0}(dx) \delta_{0}(dt)$, the solution is
$$v(x,t) = G_{t}(x),$$
which shows that the heat kernel is Green's function for the heat equation.
Best Answer
As $t \to 0$ the kernel $K(x,t)$ tends to the Dirac delta function $\delta(x)$ (which is not actually a function, but a so-called distribution, or generalized function). You can recognize this at least informally because the area under $K(x,t)$ stays constant, equal to $1$, but as $t \to 0$ it becomes more and more concentrated at the point $x =0$, with a higher and higher peak there.
Formally, what this convergence statement means is that the equation that you wrote down holds, namely that $$a(x) = \lim_{t \to 0} \int_{-\infty}^{\infty} a(x') K(x-x',t)\, dx'.$$
(Note also that in your first convoluton integral, there is a typo, as pointed out by Pacciu; the $dt$ there should be $dx'$.)
Now think about $K(x,t)$ as $t$ tends from $0$ to $\infty$: it flows from being a peak concentrated at $x = 0$ to a more-and-more shallow graph that spreads out over the whole $x$-axis. You should think of this as heat flowing from a point source at the origin and slowly spreading out over the whole axis. (Imagine blasting a point on a long steel beam with a blow torch for an instant, and then think about how the heat will diffuse along the beam.)
Now when you have initial conditions $a(x)$, this describes heat being applied at $t = 0$ not just at the point $x = 0$, but along the whole $x$-axis, according to the density $a(x)$. The convolution $a(x)*K(x,t)$ then describes how this heat has diffused through the line at time $t$; as sos440 writes in their answer, it is the superposition of the diffusion of the heat from each point $x$ that was present at time $t = 0$.
(If we took $a(x)$ to be $\delta(x)$, then we would be back at the situation of all the heat being initially concentrated at the single point $x = 0$; mathematically this corresponds to the formula $\delta(x)*K(x,t) = K(x,t)$ --- i.e. the $\delta$ function is the identity for convolution.)
Added in response to a question in the comments below:
Imagine for a moment that we had a certain amount of heat $A_i$ initially applied at the points $x_i$, for $i = 1, \ldots, n$. When one unit of heat is placed at $x =0$, it diffuses according to $K(x,t)$. So the amount $A_i$ of heat at $x_i$ diffuses according to $A_i K(x-x_i,t)$. (I am just changing the variable in $K(x,t)$ to shifts its centre from $x = 0$ to $x = x_i$, and scaling it by the amount $A_i$.)
So the total heat at a point $x$ and time $t$ will be $\sum_{i = 1}^n A_i K(x-x_i,t)$. (I am just adding up the heat which has arrived at the point $x$ at time $t$ from each of the points $x_1, \ldots, x_n$.)
Now imagine that instead of just having heat concentrated at $n$ point sources at time $t$, we instead have heat distributed throughout the line with density $a(x)$, so that the amount of heat in the (infinitesimally) small interval $[x',x' + dx']$ is $a(x') dx'$. Then the above sum becomes the integral $\int_{-\infty}^{\infty} a(x') K(x-x',t) dx'$, i.e. $a(x) * K(x,t)$.
Hence the amount of heat at a point $x$ at time $t$ is exactly given by $a(x) * K(x,t)$, as your professor explained.