Let $u$ be a classical solution for the nonhomogeneous heat equation in $\mathbb R_+ \times\mathbb R^d$: $$ \begin{cases} \partial_tu(t,x)-\Delta u(t,x) = f(t,x), \\ u(0,x)=u_0(x). \end{cases} $$ Suppose $f(t,\cdot) \in L^q(\mathbb R^n), 1 \le q \le \infty$, for each $t > 0$, and suppose $f$ has the decay behavior $$ \|f(t, \cdot)\|_{L^q} \le \frac{C}{1+t^\gamma}, \quad \gamma >0. $$ Now my question is that, given $\gamma$ and $q$, what is the decay rate of $\|u(t, \cdot)\|_{L^p(\mathbb R^d)}$ for each $1\le p \le \infty$? Probably there is no decay for all values of $p$ in this range, but I am interested in the optimal results. In particular, I would be interested in a reference (book or research paper), where this has been studied in detail.
I tried to derive the results from the representation formula for $u$: $$ u(t,x)=\int_{\mathbb R^d} H(t, x-y)u_0(y)\, dy + \int_0^t\int_{\mathbb R^d} H(t-s, x-y)f(s,y) \, dy \, ds, $$ where $ H(t,x)=(4\pi t)^{-d/2}e^{-\frac{|x|^2}{4t}} $ is the heat kernel. Here the problem comes with the latter term for large dimensions $d$. Simple arguments by using Minkowski's (or Hölder's) inequalities give some decay estimates, but then one ends up with some strange looking integrability conditions in the time integral. In particular, in this manner one only obtains decay for some values of $p$ and I doubt the estimates obtained are optimal in this sense. So I am wondering what is the correct approach and what would be the reference for this.