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.