Suppose that $\rho \in L^1(\mathbb{R}^n \times (0,T))$ for every $T < \infty$ is a weak solution of the PDE
\begin{align}
\partial_t\rho &= \Delta \rho + \text{div}(\rho\nabla\Psi(x))\\
\rho(t = 0) &= \rho^0
\end{align}
where $\Psi \in C^\infty(\mathbb{R}^n)$ is given and $\Delta$ denotes the Laplacian; that is,
\begin{equation}
\int_{0}^{\infty} \int_{\mathbb{R}^n} \rho\left(\partial\xi - \nabla\Psi \cdot\nabla\xi + \Delta\xi\right)d x d t + \int_{\mathbb{R}^n}\xi(0)\rho^0d x = 0
\end{equation}
for all $\xi \in C_c^\infty(\mathbb{R}^n\times\mathbb{R})$.

I am trying to understand how to prove that this weak solution $\rho$ is in fact a classical smooth solution. 

The paper "The Variational Formulation of the Fokker-Planck Equation" http://www.imati.cnr.it/savare/Ravello2010/JKO.pdf  sketches this argument, but I am having trouble understanding one part. In particular, they obtain the following expression on page 15, Equation (55):
\begin{align}
(\rho\eta)(t_1) &= \int_{t_0}^{t_1}\left[\rho(t)(\Delta\eta - \nabla\Psi\cdot\nabla\eta)\right]*G(t_1-t)dt\\
&\quad + \int_{t_0}^{t_1}\left[\rho(t)(2\nabla\eta - \eta\nabla\Psi\right]*\nabla G(t_1 - 1)dt\\
&\quad + (\rho\eta)*G(t_1 - t_0)
\end{align}
for all $\eta \in C_c^\infty(\mathbb{R}^n)$ and for a.e. $0 \leq t_0 < t_1$, where $G(x,t) = \frac{1}{(2\pi t)^{n/2}}e^{-|x|^2/2t}$ is the heat kernel. 

After a few straightforward computations, they show that $\rho \in L^p_{loc}(\mathbb{R}^n\times(0,\infty))$, $p < \frac{n}{n-1}$.

Then, in the following line, all they say is "We now appeal to the $L^p$-estimates [18, section 3, (3.1), and (3.2)] for the potentials
in (55) - [*the above integral equality*] - to conclude by the usual bootstrap arguments that any derivative of $\rho$ is in
$L^p_{loc}(\mathbb{R}^n\times(0,\infty))$, from which we obtain the stated regularity condition ($\rho \in C^\infty(\mathbb{R}^n\times(0,\infty)$).


Could someone provide a good explanation generally for how a bootstrapping argument would work with $L^p$ estimates, and specifically in this particular case. Also how do $L^p$ estimates from reference [18], (which is O. A. Ladyzenskaja, V. A. Solonnikov, and N. N. Uralceva, *Linear and Quasi–Linear
Equations of Parabolic Type*, American Mathematical Society, Providence, RI, 1968) come into play.

Any elucidation on this or suggestions for reading would be greatly appreciated!

NOTE: This question was originally posted on Math StackExcahnge a few months ago (https://math.stackexchange.com/questions/843901/regularity-of-a-weak-solution) and a few days ago by me (https://math.stackexchange.com/questions/941505/regularity-of-the-fokker-planck-equation?lq=1) which I have deleted since I came to know of the earlier post. Since there have been no replies  I have slightly modified the earlier question and posted it here.